Chemical potential in active systems: predicting phase equilibrium from bulk equations of state?
Siddharth Paliwal, Jeroen Rodenburg, René van Roij, Marjolein Dijkstra
aa r X i v : . [ c ond - m a t . s o f t ] J a n Chemical potential in active systems: predictingphase equilibrium from bulk equations of state?
Siddharth Paliwal
E-mail:
Soft Condensed Matter, Debye Institute for Nanomaterials Science, UtrechtUniversity, Princetonplein 5, 3584 CC Utrecht, The Netherlands
Jeroen Rodenburg, Ren´e van Roij
Institute for Theoretical Physics, Utrecht University, Princetonplein 5, 3584 CCUtrecht, The Netherlands
Marjolein Dijkstra
E-mail:
Soft Condensed Matter, Debye Institute for Nanomaterials Science, UtrechtUniversity, Princetonplein 5, 3584 CC Utrecht, The Netherlands
Abstract.
We derive a microscopic expression for a quantity µ that plays the roleof chemical potential of Active Brownian Particles (ABPs) in a steady state in theabsence of vortices. We show that µ consists of (i) an intrinsic chemical potentialsimilar to passive systems, which depends on density and self-propulsion speed, but not on the external potential, (ii) the external potential, and (iii) a newly derived one-body swim potential due to the activity of the particles. Our simulations on activeBrownian particles show good agreement with our Fokker-Planck calculations, andconfirm that µ ( z ) is spatially constant for several inhomogeneous active fluids in theirsteady states in a planar geometry. Finally, we show that phase coexistence of ABPswith a planar interface satisfies not only mechanical but also diffusive equilibrium. Thecoexistence can be well-described by equating the bulk chemical potential and bulkpressure obtained from bulk simulations for systems with low activity but requiresexplicit evaluation of the interfacial contributions at high activity.PACS numbers: 82.70.Dd,64.75.Xc, 05.65.+b, 05.40.a, 05.70.Ce, 05.10.Gg hemical potential in active systems
1. Introduction
The non-equilibrium phase behavior of active Brownian particles (ABPs), whichconstantly convert energy into directed motion, has received considerable attention inrecent years. The development of a thermodynamic framework to describe the clusteringphenomena, the pronounced accumulation of active particles at walls, and the observedcoexistence of dilute and dense phases of active matter that resemble gas-liquid andgas-solid coexistence in passive systems has been of particular interest [1–17]. Even theidea of basic thermodynamic variables such as temperature and pressure of these activesystems is being heavily debated. For instance, the effective temperature introduced byLoi et al. [18] and measured in experiments [10, 14] was shown to depend not onlyon P´eclet number, but also on the external potential and the particle interactions[1, 4, 12, 15, 19–22]. Additionally, it was argued recently that the force per unit areaon the wall can depend on the wall-particle interactions, which would imply that thepressure is not even a state function [9, 23, 24]. Similarly, a chemical potential has beenintroduced in the literature using phenomenological arguments [12, 13, 25, 26], or noiseapproximations [11] in an approach towards a thermodynamic framework for activesystems. For instance, Takatori and Brady [12] introduced a non-equilibrium chemicalpotential using micromechanical arguments, of similar form to the one that we will deriveusing the Fokker-Planck approach in this work. The authors of Ref. [12] even proceedand calculate spinodals and binodals on the basis of either a Gibbs-Duhem-like equationor a free energy for the (realistic) case of an incompressible solvent. Later, however,it was argued in Ref. [9] that a Maxwell construction on the simulated equation ofstate does not yield the simulated coexistence densities. Consequently, a complete andwell-established thermodynamic framework to describe the phase behavior of a modelas simple as ABPs is still lacking. Our Fokker-Planck approach is similar in spirit tothat of Ref. [9, 26], but defines an expression for the local chemical potential in termsof the new concept of a “swim potential”, which is well-defined in planar geometriesand curl-free particle fluxes and which may contribute, in these cases, to formulating atheoretical framework.In this study, we derive a microscopic expression for the local chemical potential µ ( z ) of active Brownian particles in a spatially inhomogeneous steady state in a planargeometry, for simplicity, with z the normal Cartesian direction. We confirm usingBrownian Dynamics simulations that µ ( z ) is spatially constant for active fluids incontact with a soft planar wall, in a gravitational field, and in two-phase coexistencewith a planar interface. Next, we show that the coexistence is described by diffusiveand mechanical equilibrium with equal bulk pressure and bulk chemical potential ofthe coexisting phases, provided the swim potential that we introduce in this article, isproperly taken into account. However, we conclude that the swim potential and hencethe chemical potential µ ( z ) is not a state function of the density for a macroscopicsystem. hemical potential in active systems
2. Methods and Formulation
We consider a three-dimensional dispersion of N active Brownian particles (ABPs) withpositions r i = ( x i , y i , z i ) and orientations ˆ e i = (sin θ i cos φ i , sin θ i sin φ i , cos θ i ) with polarangle θ i and azimuthal angle φ i , interacting via an isotropic pair potential V ( | r i − r j | )and subject to an external field V e ( r i ) for i = 1 , . . . , N at temperature T . Particle i experiences a constant self-propulsion force along its orientation ˆ e i . The motion ofparticle i is described by the overdamped Langevin equations˙ r i = − βD t ∇ i " V e ( r i ) + X j = i V ( | r i − r j | ) + v ˆ e i + p D t Ξ ti , (1)˙ˆ e i = p D r (ˆ e i × Ξ ri ) , (2)where D t and D r are the translational and rotational diffusion coefficients, β = 1 /k B T with k B the Boltzmann constant, and v is the self-propulsion speed. The collisionswith the solvent are described by a stochastic force and torque characterised by randomvectors Ξ ti and Ξ ri with h Ξ ti i = h Ξ ri i = and h Ξ ti,α ( t )Ξ tj,β ( t ′ ) i = h Ξ ri,α ( t )Ξ rj,β ( t ′ ) i = δ αβ δ ij δ ( t − t ′ ) with α, β = x, y, z [1].Starting from (1) and (2), we average over the noise to derive the deterministicFokker-Planck equation [9, 23] ∂ψ ( r , ˆ e , t ) ∂t = −∇ · j ( r , ˆ e , t ) − ∇ ˆ e · j ˆ e ( r , ˆ e , t ) (3)for the time evolution of the probability distribution function ψ ( r , ˆ e , t ) ≡ h P Ni =1 δ ( r − r i ) δ (ˆ e − ˆ e i ) i with h . . . i the averaging over the random noise. Here we defined thetranslational and rotational fluxes j = − βD t Z d r ′ Z dˆ e ′ ψ (2) ( r , ˆ e , r ′ , ˆ e , t ) ∇ V ( | r − r ′ | )+ (cid:16) − βD t ∇ V e ( r ) + v ˆ e (cid:17) ψ ( r , ˆ e , t ) − D t ∇ ψ ( r , ˆ e , t ); (4) j ˆ e = − D r ∇ ˆ e ψ ( r , ˆ e , t ) . (5)We introduced here the instantaneous full two-body correlation function ψ (2) ( r , ˆ e , r ′ , ˆ e ′ , t ) ≡ h P Ni =1 P Nj = i δ ( r − r i ) δ (ˆ e − ˆ e i ) δ ( r ′ − r j ) δ (ˆ e ′ − ˆ e j ) i , and hence to obtain a closed setof equations one needs a BBGKY-like hierarchy of Fokker-Planck equations for the n -body correlation functions or a mean-field approximation such as ψ (2) ( r , ˆ e , r ′ , ˆ e ′ , t ) ≃ ψ ( r , ˆ e , t ) ψ ( r ′ , ˆ e ′ , t ).The zeroth moment ρ ( r , t ) = R dˆ e ψ ( r , ˆ e , t ) defines the local particle density, and itstime evolution is described by the continuity equation obtained from the zeroth momentof Eq. (3), ∂ρ ( r , t ) ∂t = −∇ · J ( r , t ) , (6)with the particle flux J ( r , t ) = R dˆ e j ( r , ˆ e , t ) given by J ( r , t ) = − βD t Z ρ (2) ( r , r ′ , t ) ∇ V ( | r − r ′ | )d r ′ − βD t ρ ( r , t ) ∇ V e ( r )+ v m ( r , t ) − D t ∇ ρ ( r , t ) . (7) hemical potential in active systems ρ (2) ( r , r ′ , t ) = R dˆ e R dˆ e ′ ψ (2) ( r , ˆ e , r ′ , ˆ e ′ , t ) is the spatial two-body correlationfunction and the first moment m ( r , t ) = R dˆ e ψ ( r , ˆ e , t )ˆ e is the local polarization.An equation for m ( r , t ) follows from the first moment of Eq. (3) which yields ∂ m ( r , t ) ∂t = − ∇·J m ( r , t ) − ( d − D r m ( r , t ) , (8)with d = 2 , J m = − βD t Z dˆ e ˆ e Z d r ′ Z dˆ e ′ ψ (2) ( r , ˆ e ′ , r ′ , ˆ e ′ , t ) ∇ V ( | r − r ′ | ) − m ( r , t ) βD t ∇ V e ( r ) + v (cid:16) ρ ( r , t ) I d + S ( r , t ) (cid:17) − D t ∇ m ( r , t ) , (9)where S = R dˆ e ψ ( r , ˆ e )(ˆ e ˆ e − I /d ) is the traceless alignment tensor.We now assume that the system is only inhomogeneous in the z -direction, due toeither an external potential V e ( z ) or due to coexistence of two phases separated by aninterface parallel to the xy -plane. Without loss of generality, we consider a large, butfinite system by setting V e ( ±∞ ) = ∞ , such that ρ ( z → ±∞ ) = 0. From Eq. (7), wefind that the particle flux in the z -direction is given by J z ( z, t ) = − βD t Z d r ′ ρ (2) ( r , r ′ , t ) ∂ z V ( | r − r ′ | ) − βD t ρ ( z, t ) ∂ z V e ( z ) + v m z ( z, t ) − D t ∂ z ρ ( z, t ) . (10)When divided by βD t , we interpret Eq. (10) as a continuum force balance rather thanat the microscopic level, which requires averaging over bins that contain enough colloidsfor the continuum picture to hold. In the following sections this is achieved by havingbins that are very elongated in the direction(s) perpendicular to the z -direction.The term ( βD t ) − v m z ( z, t ) has previously been interpreted as a contribution tothe divergence of the stress tensor, which has led to a debate on pressure being a statefunction or not in active systems [23, 27, 28]. Here, however, we take another point ofview, and regard this term as an activity-induced body force − ρ ( z, t ) ∂ z V swim ( z, t ) ≡ v βD t m z ( z, t ) , (11)that is exerted on the active particles by the solvent [27, 29]. This allows us to definethe so-called swim potential V swim ( z, t ) = V swim ( z , t ) − v βD t Z zz m z ( z ′ , t ) ρ ( z ′ , t ) d z ′ , (12)where V swim ( z , t ) is a suitably chosen reference.Clearly, for a homogeneous and isotropic bulk phase, for which the polarization m = 0 in a steady state, V swim is a spatial constant. Interestingly, however, the valueof this constant is determined by surfaces and interfaces, where m can be non-zero, notunlike the Donnan potential in inhomogeneous electrolyte solutions [30, 31]. This is areflection of the fact that the activity-induced body force on the active particles onlyaverages out in the bulk, but not near interfaces. hemical potential in active systems µ ( z, t ) by J z ( z, t ) = − D t ρ ( z, t ) ∂ z βµ ( z, t ) such that µ ( z, t ) − µ ( z , t ) = µ int ( z, t ) − µ int ( z , t )+ V e ( z ) − V e ( z ) + V swim ( z, t ) − V swim ( z , t ) . (13)The external potential V e ( z ) and the intrinsic chemical potential µ int ( z, t ) = k B T ln ρ ( z, t ) + µ ex ( z, t ), consisting of an ideal part and an excess chemical potential µ ex ( z, t ), are contributions similar to those of a passive system. Here µ ex ( z, t ) is definedby µ ex ( z, t ) = µ ex ( z , t ) + Z zz d z ′′ Z d r ′ ρ ( z ′ , t ) g ( z ′′ , z ′ , R ′′ , t ) ∂ z ′′ V ( | r ′′ − r ′ | ) , (14)where we have used ρ (2) ( r , r ′ , t ) = ρ ( z, t ) ρ ( z ′ , t ) g ( z, z ′ , R, t ), with the in-plane distance R = p ( x − x ′ ) + ( y − y ′ ) , in Eq. (10). Eq. (13) reduces to the conventional chemicalpotential for a passive system, where v = 0, and is constructed such that J z = 0 if µ ( z, t ) is a spatial constant. The local chemical potential µ ( z ) is therefore a primecandidate to describe diffusive equilibrium of coexisting phases in stationary states ofactive systems. Interestingly, all terms in Eq. (13) can be determined in BrownianDynamics (BD) simulations of ABPs.The body-force interpretation of the polarization (11) can also be used to write themechanical equilibrium condition of a stationary state in terms of a well-defined normalcomponent of the stress tensor. Since the stationary state satisfies ∂ρ ( z, t ) /∂t = 0,which from Eq. (6) is equivalent to J z ( z ) = 0 for a macroscopically large, but finitesystem, we can rewrite Eq. (10) asd P N ( z )d z + ρ ( z ) ∂ z V swim ( z ) = − ρ ( z ) ∂ z V e ( z ) (15)with the standard equilibrium-like expression for the (intrinsic) normal pressure P N ( z ) = P id ( z ) + P vir ( z ) = ρ ( z ) k B T − Z z −∞ d z ′′ Z ∞ z d z ′ Z d R ′ ρ (2) ( r ′′ , r ′ ) ∂ z ′′ V ( | r ′′ − r ′ | ) , (16)where we used Newton’s third law and the symmetry of ρ (2) ( r , r ′ ) under particleexchange. The last term in Eq. (16) is the virial contribution that describes the z -component of the interparticle forces across a plane at z , which can be measured in aBD simulation [34]. Note that we did not add a swim pressure [23, 27] to the “intrinsic” P N , but instead treated the activity at the level of a swim potential V swim in the forcebalance (15), which turns out to be crucial for interpreting the (osmotic) pressure asa state function [29]. However, in order to connect to existing literature, and for laterreference, we do define P swim ( z ) − P swim ( z ) = Z zz ρ ( z ′ ) ∂ z ′ V swim ( z ′ )d z ′ = v k B T ( d − D t D r ( J m,zz ( z ) − J m,zz ( z )) (17)with the zz -component of J m given by J m,zz ( z ) = v d ρ ( z ) − m z ( z ) βD t ∂ z V e ( z ) + v S zz ( z ) − D t ∂ z m z ( z ) hemical potential in active systems − βD t Z dˆ e Z d r ′ Z dˆ e ′ ψ (2) ( r , ˆ e , r ′ , ˆ e ′ , t ) ∂ z V ( | r − r ′ | ) cos θ, (18)which reduces to the conventional swim pressure P swim ( z b ) = ρ ( z b ) v k B T / ( d ( d − D t D r )in an ideal active bulk fluid at z = z b [5, 12]. Note that our local swim pressure (17)deviates from previous expressions [35, 36] due to the gradient term ∂ z m z , which plays anon-negligible role in the force balance obtained from Eq. (15) when significant spatialvariations are present, e.g. in the interface of a phase coexistence. To summarize, wehave introduced the concept of a swim potential here using a force balance for only thecolloids. This force balance can be combined with an additional force balance for thesolvent, which provides an alternative interpretation, but identical expression, for theswim pressure as an excess solvent pressure [29].With the definition (17) one can thus define a total pressure P ( z ) = P N ( z ) + P swim ( z ), such that Eq. (15) can be written as d P/ d z = − ρ ( z ) ∂ z V e ( z ); in the case where V e ( z ) = 0 a steady state is then characterized by a spatially constant total pressure P ( z ). The intrinsic chemical potential µ int ( z ) and intrinsic normal pressure P N ( z ), andthe swim potential V swim ( z ) and swim pressure P swim ( z ) have thus been constructed suchthat d P N ( z )d z = ρ ( z ) d µ int ( z )d z , and d P swim ( z )d z = ρ ( z ) d V swim ( z )d z . (19)If we now invoke a Local Density Approximation (LDA), i.e. assume that the localenvironment behaves as a bulk such that the local pressure and chemical potential area function of only the local density ρ ( z ), then Eq. (19) can be written in terms of bulkquantities as: d P N ( ρ )d ρ = ρ d µ int ( ρ )d ρ , and d P swim ( ρ )d ρ = ρ d V swim ( ρ )d ρ , (20)allowing us to writed P ( ρ )d ρ = ρ d µ ( ρ )d ρ (21)with µ ( ρ ) = µ int ( ρ ) + V swim ( ρ ), in a zero external potential. Here, we shall take care todistinguish the notation µ ( ρ ) for the chemical potential obtained via Eq. (21) from µ ( z )which denotes the chemical potential calculated from Eq. (13). We recognize Eq. (21)as a generalization of the Gibbs-Duhem relation for equilibrium systems. Whereas inequilibrium (where P swim = V swim = 0) it holds true in general, we emphasize that inthis case we had to make use of the LDA to derive it. This Gibbs-Duhem relationprovides a way to obtain the chemical potential µ ( ρ ) from the bulk equation of state P ( ρ ), whereas to obtain µ ( z ) from Eq. (13) we require complete spatial profiles. Wetest the applicability of Eq. (21) in simulations and show that it works well for caseswith low anisotropy (e.g. low polarization). However, Eq. (21) does not hold true ingeneral as V swim ( z ) = V LDA swim ( ρ ( z )) for high anisotropy as we discuss later.We note that Eq. (21) is akin to the one in Ref. [12], apart from a factor thatis equal to the (incompressible) solvent volume fraction. The equilibrium analogue ofEq. (21) follows naturally if the solvent is treated grand-canonically which we implicitly hemical potential in active systems µ ( z ) in both the Fokker-Planck calculations and particlebased simulations. Next we present the results of BD simulations of an active fluidwith Lennard-Jones (LJ) interactions subject to a gravitational field in Section 3.2. InSection 3.3 we consider an active Lennard-Jones fluid exhibiting gas-liquid coexistencewith a planar interface and confirm mechanical and diffusive equilibrium. We performa Maxwell equal-area construction to identify phase coexistence from bulk equations ofstate. We then attempt to apply the same formalism to active particles which undergoMotility Induced Phase Separation at high activity in Section 3.4.
3. Results
We first consider a three-dimensional active ideal gas (with V ( r ) = 0) at P´eclet numberPe = v /σD r = 0 (passive), 1, 3, 5, in the external potential βV e ( z ) = ( z/σ ) for z < V e ( z ) = 0 for z >
0, where the unit of length σ = p D t /D r is chosen to be theparticle diameter so that the Stokes-Einstein relation for spheres in three dimensionsis satisfied. Note that Pe can also be perceived as the ratio of the persistence length v /D r and the particle diameter [5]. For large but finite z = z b & σ , the active fluidreaches a bulk state with bulk density ρ b = ρ ( z b ), and the normal pressure reducesto the bulk pressure P b = P N ( z b ) = ρ b k B T . In Fig. 1(a) and (b) we show the time-averaged density profiles ρ ( z ) and orientation profiles m z ( z ) /ρ ( z ), respectively. Weobserve that the particles penetrate deeper into the wall at higher Pe resulting intoa more extended ρ ( z ) within the wall accompanied by a small adsorption (that wasfound in Ref. [37] as well) close to z = 0. In Fig. 1(b) we see no average polarizationoutside or inside the wall for the passive case. At finite Pe, however, Fig. 1(b) showsthat the average orientation is zero in the bulk where V e ( z ) = 0 and negative withinthe wall, corresponding to particles oriented towards the wall. Fig. 1(c) and (d) show V swim ( z ) and µ ( z ) as obtained from Eq. (12) and (13), respectively. We find that µ ( z )is indeed constant within our statistical accuracy of ∼ . k B T . Clearly, for µ ( z ) to beconstant it is crucial that V swim ( z ), which is attractive towards the wall consistent withthe polarization and extended density profile close to the wall, is included in Eq. (13);ignoring this contribution of 10 − k B T would not have yielded a spatially constantchemical potential in the stationary state. Although µ ( z ) was constructed to be spatiallyconstant within the Fokker-Planck formalism, a confirmation from the simulations serves hemical potential in active systems . . . . . . . ρ ( z ) / ρ b (a) P e β V s w i m ( z ) (c) βV e ( z ) − − − − − z/σ − . − . − . − . . m z ( z ) / ρ ( z ) (b) − − − − − z/σ − . . . . . . β µ ( z ) (d) Figure 1. (a) Density profile ρ ( z ), (b) polarization profile m z ( z ) /ρ ( z ), (c) swimpotential V swim ( z ) and the soft external potential V e ( z ) (see text), and (d) localchemical potential µ ( z ) (with error bars), all as a function of z for an active idealgas in contact with a planar soft wall, as obtained from BD simulations (solid lines)and Fokker-Planck calculations (dashed lines), for varying P´eclet numbers as labeled.In (d) the solid lines represent βµ ( z ) obtained by the integration of J z ( z ) /ρ ( z ), whichfluctuates about zero, whereas the square symbols show the resultant from Eq. (13).The errorbars represent the error induced in βµ ( z ) due to the statistical error in ρ ( z ).The deviation from the Fokker-Planck calculations deep into the wall for high Pe isdue to the correlation of error upon integration. as a useful validation.Additionally, we verify that the swim pressure (given by Eq. (17)) measured inthe bulk reduces to P swim ( z b ) = k B T v ρ b / (6 D r D t ). V swim can similarly be obtainedas V swim ( z b ) = ( k B T v / D r D t ) ln ρ b σ . We use this bulk state at z b & σ with V swim ( z = z b ) as the reference point for the profiles of V swim ( z ) and µ ( z ) in Fig. 1(c)and (d), respectively. We now consider simulations of weakly active Lennard-Jones (LJ) particles with anisotropic pair potential, V LJ ( r ) = 4 ǫ (( σ/r ) − ( σ/r ) ), at k B T /ǫ = 2 . V e ( z ) = M gz for z > z = 0, with M the buoyant particle mass. These systems are supercritical in the passive case, andtherefore even more so in the active cases since the ‘critical temperature’ decreases withincreasing activity [6,17]. We measure the density ρ ( z ) σ , polarization m z ( z ) /ρ ( z ), swim hemical potential in active systems . . . ρ ( z ) σ (a) P e, βM gσ (0 , . , .
0) (10 , ,
5) (20 , , z/l β µ ( z ) (d) . . . . m z ( z ) / ρ ( z ) (b) β V s w i m ( z ) (c) βV e ( z ) .
001 0 .
01 0 . ρσ . . . . . β P N ( ρ ) / ρ (e) P e, βM gσ (0 , . , . , , , , . . . . . ρσ β µ i n t ( ρ ) (f) Figure 2.
Height-dependence of (a) density ρ ( z ), (b) polarization m z ( z ) /ρ ( z ),(c) swim potential V swim ( z ), and (d) chemical potential µ ( z ) (with an offset for clarity),all for an active LJ fluid in an external gravitational potential V e ( z ) = M gz for variousvalues of βM gσ , and P´eclet number Pe=0 (blue), 10 (green), 20 (red) as obtainedfrom BD simulations. The height z is scaled with respect to l , where l = v /D r is thepersistence length for Pe=10 and 20, and l = σ is the particle diameter for Pe=0. Thecompressibility factor P N ( ρ, P e ) /ρ in (e) and the intrinsic chemical potential µ int ( ρ, P e )shown with an offset in (f) show a proper collapse in the dilute limit for different βM gσ but not for Pe. potential βV swim ( z ), and chemical potential µ ( z ) for βM gσ = 0 . βM gσ = 3 and 5 for Pe=10 and 20, all plotted in Fig. 2(a)-(d). In order to obtain acomparable length scale l over which variations are observed in the passive (where wechoose l = σ ) and in the active cases (where l = v /D r ), we used a smaller buoyantmass of the particles in the passive case. We observe that the polarization m z ( z ) ispositive for Pe=10 and 20, and hence the mean swimming direction is opposite to thegravitational field, consistent with the findings in Ref. [14]. Moreover, Fig. 2(b) showsthat the polarization profile m z ( z ) /ρ ( z ) is surprisingly constant over a large regime ofheights z . As a consequence, the swim potential profile βV swim ( z ) essentially decreaseslinearly with height z for Pe = 10 and 20 and counteracts largely the gravitationalfield, as shown in Fig. 2(c), leading to an enormous increase in sedimentation length hemical potential in active systems βM g ) − [10]. The chemical potential profile µ ( z ) is calibrated by µ ( z ) = 0 at thereference point z determined by the condition ρ ( z ) σ = 5 × − . µ ( z ) is shown inFig. 2(d) and is indeed spatially constant within our statistical accuracy of ∼ . k B T .It is important to note here that V swim ( z ) decreases by a few hundred k B T and theexternal gravitational potential V e ( z ) = M gz increases by a few hundred k B T in the z -range of interest as shown in Fig. 2(d).In addition, we show in Fig. 2(e) and (f) both P N and µ int as a function of ρ ,obtained by eliminating z from P N ( z ) and ρ ( z ), and µ int ( z ) and ρ ( z ), respectively. Weobserve that the data collapse at fixed Pe, and it is alluring to interpret that P N ( ρ, Pe)and µ int ( ρ, Pe) are state functions of the density in this regime.
We now consider a weakly active LJ fluid without any external potential ( V e ( z ) = 0),and at subcritical temperatures such that coexistence of a gas and a liquid phase withbulk densities ρ g and ρ l , respectively, is to be expected at overall intermediate densities ρ g < ρ < ρ l in an elongated simulation box with periodic boundary conditions [6, 17].A temperature k B T /ǫ = 0 .
43 and a P´eclet number Pe= v /D r σ = 3 are used in thiscase. In Fig. 3(a), we show a typical configuration of a liquid slab in the center of thesimulation box in coexistence with a gas phase on either side. In Fig. 3(b) we plot thecorresponding density profile ρ ( z ) which can be fitted to a hyperbolic tangent function(Eq. (A.1)), independently for z > z <
0, to obtain the coexistence densities ρ ( z g ) and ρ ( z l ) of the two bulk phases as fit parameters, with z g and z l a position inthe bulk gas and liquid respectively. In the same figure we also plot the polarizationprofile m z ( z ) /ρ ( z ), showing that the swimming direction of the particles at the liquid-gas interface is pointing from the liquid phase towards the gas phase, i.e., against theattractive interparticle forces from the liquid [17, 38].In Fig. 3(c) and (d) we plot the profiles P ( z ) and µ ( z ), respectively, which clearlyshow that both are spatially constant. We hence conclude that P ( z g ) = P ( z l ) and µ ( z g ) = µ ( z l ), demonstrating mechanical and diffusive equilibrium of the coexisting gasand liquid phase. For completeness, in Fig. 3(c) we also plot the individual contributionsto the total pressure P ( z ) = P N ( z ) + P swim ( z ), where P swim ( z ) is the swim pressureobtained from Eq. (17), and P N ( z ) = P N, idl ( z ) + P N, vir ( z ) is the normal pressurewith the ideal pressure P N, idl ( z ) and the virial contribution to the normal pressure P N, vir ( z ) as obtained from Eq. (16). Similarly we plot the contributions to the chemicalpotential µ ( z ) = µ int ( z ) + V swim ( z ) in Fig. 3(d), where the intrinsic chemical potential µ int ( z ) = k B T ln ρ ( z ) + µ ex ( z ) represents the sum of ideal and excess chemical potential.The swim potential V swim ( z ) is calculated from the measured polarization profiles usingEq. (12).In order to investigate if we can predict phase coexistence solely from bulkquantities, we perform BD simulations of bulk states of ABPs at several temperatures k B T /ǫ and P´eclet number Pe=2.67. We measure the bulk pressure P as a function of hemical potential in active systems (a) − . . . . . . (b) ρ ( z ) σ m z ( z ) /ρ ( z ) − . − . − . . . . (c) βP N, idl ( z ) σ βP N, vir ( z ) σ βP N ( z ) σ βP swim ( z ) σ βP ( z ) σ − −
10 0 10 20 z/σ − − (d) ln ρ ( z ) σ βµ ex ( z ) βµ int ( z ) βV swim ( z ) βµ ( z ) -20 0 20-0.40.00.4 Figure 3. (a) A typical configuration of a three-dimensional gas-liquid coexistenceof an active LJ fluid at Pe=3, and temperature k B T /ǫ = 0 .
43, along with (b)the corresponding density profile ρ ( z ) and polarization profile m z ( z ) /ρ ( z ), (c) totalpressure P ( z ) = P N ( z ) + P swim ( z ) and the individual contributions, and (d)total chemical potential µ ( z ) obtained from Eq. (13), and individual contributions,along with an inset showing a magnified view of µ ( z ). Both P ( z ) and µ ( z ) arespatially constant within numerical accuracy, demonstrating mechanical and diffusiveequilibrium of the coexisting gas and liquid phase. density ρ in a simulation box small enough to prevent phase separation and plot theequations of state P ( ρ ) for several subcritical temperatures in Fig. 4(a). Now, withina Local Density Approximation (LDA), we apply the Gibbs-Duhem relation Eq. (21)and obtain µ ( P ) by integrating the equation of state P ( ρ ) for several T ’s as shownin Fig. 4(b). We emphasize here that we refer to µ ( ρ ) as the µ obtained by applyingEq. (21) which is not to be confused with µ ( z ). The intersection of the curve µ ( P ) givesthe coexistence µ g = µ l and P g = P l . In the inset of Fig. 4(b) we compare the binodalsin the (scaled) temperature-density plane as obtained from the density profiles from hemical potential in active systems . . . . . ρσ − . − . − . . . . β P σ (a) k B T /ǫ . . .
48 0 . . − . − . − . . . βP σ β µ ( P ) (b) ρσ k B T / ǫ directMaxwell Figure 4. (a) Scaled pressure-density P - ρ , and (b) chemical potential-pressure µ - P relations of an active LJ fluid at several temperatures k B T /ǫ and Peclet Pe= v /σD r =2 .
67. The inset shows the temperature-density gas-liquid binodals as obtained fromdirect coexistence simulations ( (cid:4) ) and from equating µ and P in the coexisting phases( N ) of an active LJ fluid. direct coexistence simulations ( ρ ( z g ) and ρ ( z l )) and from the bulk µ ( P ) intersections( ρ g and ρ l ). We find good agreement between the two results and thus conclude thatthe corresponding coexistence densities ρ g and ρ l could, in this (low Pe) case at least,be determined from the bulk equations of state. Note that the activity has a huge effecton the gas-liquid binodals (shown in the inset of Fig. 4(b)) as the critical temperatureshifts from k B T /ǫ ≈ .
15 in the passive case to k B T /ǫ ≈ .
54 in the active case forPe=2.67 (see Ref. [17] for full comparison).
In this section we discuss the swim potential and the chemical potential in a two-dimensional system of strongly active particles exhibiting motility induced phaseseparation at high Pe. We choose our planar geometry in the yz plane and assumehomogeneity in the y direction to be consistent with previous definitions. The particlesinteract with the WCA potential given by V W CA ( r ) = V LJ ( r ) + ǫ , with a cut-off beyond r ≥ r c = 2 / σ to make the particles purely repulsive. The particle orientations can bedescribed in terms of a single angle θ i as ˆe i = (cos θ i , sin θ i ). The translational equationof motion in 2D is similar to Eq. (1) and the rotational diffusion follows ˙ θ i = √ D r Ξ ri ,with Ξ ri a zero-mean unit-variance Gaussian random variable.As before, we fix rotational and translation diffusion coefficients to correspondto the particle interaction length scale σ = p D t /D r and change the self-propulsionspeed v to vary Pe. At high Pe, we find that the system phase separates into a gasphase and a dense phase, both of well-defined densities, separated by a planar interfacein an elongated simulation box [16]. For Pe=50 the typical density and polarizationprofiles are shown in Fig. 5(a). Notably, the polarization profiles are now reversed with hemical potential in active systems − −
50 0 50 100 z/σ − . . . . . (a) ρ ( z ) σ m z ( z ) /ρ ( z ) (c) βP N,idl ( z ) σ βP N,vir ( z ) σ βP swim ( z ) σ βP ( z ) σ . . . . ρσ (b) βP ( ρ ) σ βP maxwell σ βP coex σ z/σ − − (d) ln ρ ( z ) σ βµ ex ( z ) βV swim ( z ) βµ ( z ) -100 0 100-303 ρσ P e directMaxwell Figure 5. (a) Density ρ ( z ) σ and polarization m z ( z ) /ρ ( z ) profiles of an active fluidwith WCA interactions exhibiting MIPS at Pe = 50, and temperature k B T /ǫ = 0 . βP ( ρ ) σ vs. density ρσ curve obtained from bulk simulations of smallsystems (solid circles) and large systems (open circles), together with Maxwell equal-area pressure (dashed line) and coexistence pressure P coex = P ( z ) (dotted line) asmeasured in (c). The inset shows a comparison of bulk densities from direct coexistencesimulations ( (cid:4) ) and the Maxwell equal-area construction ( N ) for various Pe. (c) Totalpressure P ( z ) = P N ( z )+ P swim ( z ) profile, with the ideal, virial and swim contributions,and (d) total chemical potential µ ( z ) profile with individual contributions, for z > βµ id ( z ) = ln ρ ( z ) σ and that µ ( z ) is constant within an accuracy of 3 k B T . respect to Fig. 3(b) as the particles at the interfaces are now pointing towards the densephase. We measure the normal component of the total pressure P ( z ) and the chemicalpotential µ ( z ) by summing the individual contributions, and plot them in Fig. 5(c) and(d), respectively. We clearly observe that both the quantities P ( z ) and µ ( z ) are spatiallyconstant, demonstrating mechanical and diffusive equilibrium of the coexisting phases.With the polarization profiles reversed, P swim ( z ) and V swim ( z ) are now higher in the gasphase as compared to the denser phase.Further, we perform a Maxwell equal-area construction on the equation of state.The P − ρ curves shown in Fig. 5(b) are obtained again using a small system size forwhich there is no global phase separation at intermediate densities. We confirm theresults of the homogeneous states with larger system sizes and find that the agreementis satisfactory for our analysis. Performing a Maxwell construction on P as a functionof 1 /ρ gives the equal-area pressure P Maxwell shown as the dashed horizontal line in hemical potential in active systems P coex obtained fromthe direct coexistence simulation of the phases coexisting at the corresponding set ofparameters. From the two curves it is evident that the coexistence densities predictedby the Maxwell construction and the direct-coexistence simulations do not agree. Weperform the same procedure on a set of Pe in the range 30 −
60 and plot the correspondingcoexistence densities and the densities predicted by the Maxwell construction in theinset of Fig. 5(b). From the disagreement between the two binodals we conclude thatthe Maxwell equal-area construction does not correspond to the coexisting states asobtained from the direct coexistence simulations, noted previously as well in Ref. [9,26].We have checked that using our P ( ρ ) data with the definition of the chemical potentialintroduced in Ref. [12] yields the same binodals as predicted here despite the differenceof the factor concerning the solvent volume fraction.
4. Discussion
The results from the previous section show that the Maxwell equal-area construction,and hence the Gibbs-Duhem equation (20), cannot be used in general to predictthe coexisting densities ρ g and ρ l [9, 26] in systems of ABPs. In other words, eventhough µ ( z g ) = µ ( z l ) in a phase-separated system (where z g and z l are locationsfar from interfaces such that the local densities are ρ g and ρ l , respectively) , thechemical potentials obtained from the Gibbs-Duhem equation (21) may not be equal,i.e. µ ( ρ g ) = µ ( ρ l ). The nonzero difference between µ ( ρ g ) and µ ( ρ l ) is caused by thefailure of the LDA assumed in the derivation of Eq. (21), as we will show below. Inparticular, the values of V swim ( z ) and µ ex ( z ) in a bulk state at position z and density ρ b do not only depend on ρ b (and other system parameters such as Pe) but also on theinterface between the bulk state and the reference state at z . This implies that neither V swim nor µ ex as expressed in Eqs. (12) and (14), respectively, are state functions of thedensity. Below we show an example for V swim ( z ) which demonstrates this breakdownof the LDA in the case of a 2D active ideal gas (for which µ ex ( z ) ≡
0) in a particularexternal potential.The setup consists of a ramp-like external potential βV e ( z ) = λz/σ in the region0 < z < σ which separates a bulk region at the left (where βV e ( z ) = 0 for z <
0) fromthe bulk on the right (where βV e ( z ) = 5 λ for z > σ ). These external potential areplotted in Fig. 6(a) as dash-dot lines for λ = 0 , . ,
1, and 2. The probability density ψ ( z, θ ) is obtained by solving Eq. (3) for V ( r ) = 0 numerically, at Pe = 1 with a fixeddensity boundary condition ρσ = 1 . z = − σ and with a hard wall placed at z = 15 σ . The density and polarization profiles for increasing λ are plotted in Fig. 6(a)and (b), respectively.In order to determine V swim ( z ) for this non-interacting system with V ( r ) ≡
0, Eq. (17)can be rewritten as V swim ( z ) − V swim ( z ) = v βD t D r ln (cid:18) ρ ( z ) ρ ( z ) (cid:19) hemical potential in active systems . . . . . ρ ( z ) σ (a) . . . . . m z ( z ) / ρ ( z ) (b) λ = 0 . λ = 0 . λ = 1 . λ = 2 . − − z/σ − − − − β V s w i m (c) V swim ( z ) V LDA swim ( ρ ( z )) β V e ( z ) ρ ( z ) V e ( z ) Figure 6. (a) Density profiles ρ ( z ) and (b) polarization profiles m z ( z ) /ρ ( z ) of anon-interacting active fluid at Pe=1 in a ramp-shaped external potential with slope λ = 0 , . , , V swim ( z ) obtained usingEq. (12) and βV LDA swim ( ρ ( z )) = ( v / D t D r ) ln ρ ( z ) σ obtained using LDA. + v βD r Z zz ρ ( z ′ ) dd z ′ (cid:20) v D t S zz ( z ′ ) − βm z ( z ′ ) ∂V e ( z ′ ) ∂z ′ − ∂m z ( z ′ ) ∂z ′ (cid:21) d z ′ . (22)The V swim ( z ) profiles, obtained equivalently from Eq. (22) or from Eq. (12), are plottedas solid lines in Fig. 6(c) where we have taken z = − σ as the reference state where V swim ( z ) = 0. If we would approximate the vicinity of any point z ′ as an isotropic bulkwith density ρ ( z ′ ) in evaluating the swim potential V swim ( z ′ ), i.e. assume in Eq. (22) S zz ( z ′ ) ≈ m z ( z ′ ) ≈ z ′ ,we obtain βV LDA swim ( ρ ( z )) = ( v / D t D r ) ln ρ ( z ) σ which we refer to as the local densityapproximation (LDA) of Eq. (22). Note that Eq. (22) follows from the Fokker-Planckformalism, and this LDA does not refer to an approximation of a free-energy functional.This V LDA swim ( ρ ( z )), plotted as dotted lines in Fig. 6(c), is equal to V swim ( ρ ) obtained fromthe swim component of the Gibbs-Duhem-like relation (21). We find that V swim ( z ) and V LDA swim ( ρ ( z )) start to deviate at high λ and do not coincide in the right bulk. Hence, wecan conclude that the values for V LDA swim obtained from the Gibbs-Duhem equation are notcorrect in general. This is due to the failure of LDA, i.e. due to the anisotropy in theinterface that renders the integral on the right hand side in Eq. (22) non-negligible as hemical potential in active systems λ , consistent with this idea of increasing anistropy. For an interactingsystem the forces between particles would add another contribution to m z ( z ) /ρ ( z ), whichcould also become a source of failure for the LDA. This will be studied in more detailin a future publication.In Section 3.3 we observed that the Maxwell construction was able to predictthe coexistence densities for the active LJ case with reasonable accuracy, but was indisagreement at higher activity in Section 3.4 for MIPS. We now assert that the errormade in the chemical potential by assuming the LDA translates into an error in thepredicted coexisting densities that is small for the active LJ particles, but significant forMIPS. We define the error in predicted coexistence densities of the gas and the densephase, respectively, as ∆ ρ errg = ρ ( z g ) − ρ g and ∆ ρ errl = ρ ( z l ) − ρ l , where ρ ( z g ) and ρ ( z l ) are the bulk coexistence densities and ρ g and ρ l denote the estimates obtained byperforming a Maxwell construction. If we define the gas state as the reference statefor the chemical potential, i.e. µ ( z ) = µ ( ρ g ) in Eq. (13) with z = z g , then the errormade in determining the chemical potential of the dense phase by using the Gibbs-Duhem equation (21) is ∆ µ errl = µ ( ρ l ) − µ ( z l ), where we recall that µ ( ρ l ) is the chemicalpotential of the dense phase obtained from the Gibbs-Duhem relation, whereas µ ( z l )is the true chemical potential determined in the coexistence simulation. From ∆ µ errl the relative error in the predicted density of the dense phase can be estimated as∆ ρ errl /ρ ( z l ) ≈ /ρ ( z l ) · ∆ µ errl / ( dµ/dρ ) l . Similarly, the error in the predicted densityof the gas phase can be estimated by using the dense phase as the reference state( µ ( z ) = µ ( ρ l )). The relative density error estimated in this manner is less than 5% forthe active LJ case, whereas it is of the order of 100% for the MIPS case, which agreeswith our findings in Fig. 4(b) and 5(b), respectively.We wish to make a note that the anisotropy terms identified here resemblethe interfacial contributions discussed in Ref. [26] for pairwise-interacting particles.Although it requires explicit measurement of these interfacial contributions byperforming phase-coexistence simulations, Solon et al. were able to suggest a modifiedMaxwell construction for estimating the binodals in Ref. [26].Moreover, our elongated simulation box in Section 3.3 and 3.4 forces the systemto phase separate with a planar interface. Only for such a geometry J z ( z ) = 0,allowing us to write explicit expressions for mechanical and diffusive equilibrium. Inother geometries the stationary state condition ∇ · J = 0 still allows for swirls thatcorrespond to non-zero ∇ × J , for which our expressions for mechanical and diffusiveequilibrium break down and a whole new framework is needed. Furthermore, the regimeof applicability of Eq. (13) is limited by the underlying dynamic DFT relation, wherea ρ -independent diffusion coefficient D t is assumed; an extension to account for a ρ -dependent diffusion coefficient is left for a future study. hemical potential in active systems
5. Conclusions
In conclusion, we have constructed expression (13) for the local chemical potential µ ( z )for active fluids in a planar geometry, which includes the swim potential V swim ( z ) definedby Eq. (12) in addition to ideal, excess, and external contributions well-known fromequilibrium. Our BD simulations confirm that µ ( z ) is spatially constant in steady statesof several inhomogeneous ideal and interacting fluids of active particles, with V swim ( z ) animportant contribution that counteracts either the external potential V e ( z ) or the excesscontribution µ ex ( z ). In the low activity regime studied for active LJ fluid, the chemicalpotential provides a method to predict the coexisting densities from bulk simulations.At high activity the anisotropy in the interface causes the Gibbs-Duhem relation to beinvalid, which provides support to the conclusions of Ref. [26] that the details of theinterface are necessary to determine the coexisting bulk densities. Our formalism opensnew avenues towards a Fokker-Planck and dynamic density functional description (ofstationary states) of active systems, especially for planar geometries. Acknowledgments
S.P. and M.D. acknowledge funding from the Industrial Partnership Programme“Computational Sciences for Energy Research” (Grant No.14CSER020) of theFoundation for Fundamental Research on Matter (FOM), which is part of theNetherlands Organization for Scientific Research (NWO). This research programme isco-financed by Shell Global Solutions International B.V. This work is part of the D-ITPconsortium, a program of the Netherlands Organisation for Scientific Research (NWO)that is funded by the Dutch Ministry of Education, Culture and Science (OCW). Weacknowledge funding of an NWO-VICI grant. We would like to thank V. Prymidis forproviding initial estimates and configurations for simulating MIPS in 2D, and also L.Filion and J. Tailleur for stimulating discussions.
Appendix A. Simulation details
We perform Brownian Dynamics (BD) simulations for three-dimensional and two-dimensional geometries in Sections 3.1-3.3 and Section 3.4, respectively. We use theEuler-Maruyama method to integrate the equations of motion (1) and (2) with a timestep size dt = 10 − τ where τ = 3 Dr − is the unit of time. We keep the temperature ofthe bath fixed at T by keeping the translational and rotational diffusion coefficients ( D t and D r respectively) fixed and vary v to change Pe and interaction strength ǫ to changethe temperature of the colloidal particles. We employ periodic boundary conditions inonly x - and y -direction in Section 3.1 and 3.2, in all three directions in Section 3.3, andin both y - and z - directions in Section 3.4. The system sizes are about 2500 particlesin 3D and about 6500 particles in 2D for elongated box simulations. We measure thedensity profile ρ ( z ) in the z -direction as ρ ( z ) = h n ( z ) i /L ∆ z by measuring the average hemical potential in active systems h n ( z ) i in the slabs of volume L ∆ z ( ρ ( z ) = h n ( z ) i /L ∆ z in2D) arranged parallel to xy plane ( y -direction in 2D), where L is the length of the systemin the x and/or y -direction, and where ∆ z = 0 . σ is the width of the slab. In a similarmanner we measure the polarization profile m z ( z ) by summing the particle orientationsin a slab at location z . The density profiles ρ ( z ) can be fitted to a hyperbolic tangentfunction given by: ρ ( z ) = 12 ( ρ ( z l ) + ρ ( z g )) + 12 ( ρ ( z l ) − ρ ( z g )) tanh (cid:20) z − z ∗ ) D (cid:21) , (A.1)where ρ ( z l ) and ρ ( z g ) are the corresponding bulk liquid and vapour coexisting densities, z ∗ is the location of the dividing plane and D represents the thickness of the interface.Subsequently, the swim potential profile V swim ( z ) is obtained as V swim ( z ) = V swim ( z ) − v βD t Z zz m z ( z ′ ) ρ ( z ′ ) d z ′ , (A.2)where we use ρ ( z ) and m z ( z ) as measured in the BD simulations, and where V swim ( z )is a suitably chosen reference state. In addition, we measure the normal component ofthe stress tensor using P N ( z ) = P id ( z ) + P vir ( z ) (A.3)with the ideal gas pressure P id ( z ) and the virial pressure P vir ( z ) given by: P id ( z ) = ρ ( z ) k B T (A.4) P vir ( z ) = − L ∆ z * N X i =1 N X j = i z ij r ij d V ( r ij )d r ij Z C ij ˆ z · d l + , (A.5)where r ij = | r ij | = | r j − r i | denotes the center-of-mass distance between particle i and j , z ij = z j − z i where z i is the z position of particle i , C ij is the intersection of r ij andthe slab of width ∆ z centered at z . The integral in Eq. (A.5) denotes that the virialcontribution to the pressure of particle pair i and j is due to the part of r ij that lies insidethe respective slab at z within the coarse-grained Irving-Kirkwood approximation [34].We also calculate the swim pressure P swim ( z ) = v k B Td ( d − D t D r ρ ( z ) − v m z ( z )( d − D r ∂ z V e ( z ) − v ( d − D r Z d ˆe Z d r ′ d ˆe ′ ψ (2) ( r , ˆe , r ′ , ˆe ′ )( ∇ V ( | r − r ′ | ) · ˆ z ) cos θ + v k B T ( d − D t D r S zz ( z ) − v k B T ( d − D r ∂ z m z ( z ) , (A.6)and the chemical potential profile µ ( z ) using µ ( z ) = µ ( z ) + k B T ln ρ ( z ) + µ ex ( z ) + V e ( z ) + V swim ( z ) − k B T ln ρ ( z ) − µ ex ( z ) − V e ( z ) − V swim ( z ) , (A.7)with the excess chemical potential µ ex ( z ) defined as µ ex ( z ) = µ ex ( z ) + Z zz d z ′ * n ( z ′ ) n ( z ′ ) X i =1 N X j = i z ′ ij r ′ ij d V ( r ′ ij )d r ′ ij + . (A.8) hemical potential in active systems z with respect to a reference at z is determinedby integrating the averaged force that a particle feels due to the particle interactionswith all other particles in the system over the distance z to z .Alternatively, if V e ( z ) = 0, µ ( z ) can also be obtained using µ ( z ) = µ ( z ) + Z zz d z ′ ρ ( z ′ ) d P ( z ′ )d z ′ (A.9)with P ( z ) = P N ( z ) + P swim ( z ). References [1] Fily Y and Marchetti M C 2012
Physical Review letters
Physical Review letters
EPL (Europhysics Letters)
Physical Review letters
Physical Review Letters
Physical Review E Soft Matter Nature Communications Physical Review Letters
Physical Review Letters
Soft Matter Physical Review E Soft Matter Physical Review Letters
EPL (Europhysics Letters)
Annual Review of Condensed Matter Physics The Journal of Chemical Physics
Soft Matter (8) 3726[19] Szamel G 2014 Physical Review E Physical Review X The European Physical Journal Special Topics
Current Opinion in Colloid & Interface Science Nature Physics Physical Review Letters
Soft Matter ArXiv e-prints ( Preprint )[27] Yan W and Brady J F 2015
Soft Matter Physical Review E Soft Matter (47) 8957[30] Verwey E and Niessen K 1939 The London, Edinburgh, and Dublin Philosophical Magazine andJournal of Science Physical Review Letters
The Journal of Chemical Physics
The Journal of Chemical Physics
Molecular Simulation hemical potential in active systems [35] Yang X, Manning M L and Marchetti M C 2014 Soft Matter Soft Matter Journal of Fluid Mechanics
R1[38] Paliwal S, Prymidis V, Filion L and Dijkstra M 2017
The Journal of Chemical Physics147