Chemical potential in active systems: predicting phase equilibrium from bulk equations of state?

We derive a microscopic expression for a quantity $\mu$ that plays the role of chemical potential of Active Brownian Particles (ABPs) in a steady state in the absence of vortices. We show that $\mu$ consists of (i) an intrinsic chemical potential similar 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 active Brownian particles show good agreement with our Fokker-Planck calculations, and confirm that $\mu(z)$ is spatially constant for several inhomogeneous active fluids in their steady states in a planar geometry. Finally, we show that phase coexistence of ABPs with a planar interface satisfies not only mechanical but also diffusive equilibrium. The coexistence can be well-described by equating the bulk chemical potential and bulk pressure obtained from bulk simulations for systems with low activity but requires explicit evaluation of the interfacial contributions at high activity.


Introduction
The non-equilibrium phase behavior of active Brownian particles (ABPs), which constantly convert energy into directed motion, has received considerable attention in recent years. The development of a thermodynamic framework to describe the clustering phenomena, the pronounced accumulation of active particles at walls, and the observed coexistence of dilute and dense phases of active matter that resemble gas-liquid and gas-solid coexistence in passive systems has been of particular interest [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17]. Even the idea of basic thermodynamic variables such as temperature and pressure of these active systems is being heavily debated. For instance, the effective temperature introduced by Loi et al. [18] and measured in experiments [10,14] was shown to depend not only on Péclet number, but also on the external potential and the particle interactions [1,4,12,15,[19][20][21][22]. Additionally, it was argued recently that the force per unit area on the wall can depend on the wall-particle interactions, which would imply that the pressure is not even a state function [9,23,24]. Similarly, a chemical potential has been introduced in the literature using phenomenological arguments [12,13,25,26], or noise approximations [11] in an approach towards a thermodynamic framework for active systems. For instance, Takatori and Brady [12] introduced a non-equilibrium chemical potential using micromechanical arguments, of similar form to the one that we will derive using the Fokker-Planck approach in this work. The authors of Ref. [12] even proceed and calculate spinodals and binodals on the basis of either a Gibbs-Duhem-like equation or 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 of state does not yield the simulated coexistence densities. Consequently, a complete and well-established thermodynamic framework to describe the phase behavior of a model as simple as ABPs is still lacking. Our Fokker-Planck approach is similar in spirit to that of Ref. [9,26], but defines an expression for the local chemical potential in terms of the new concept of a "swim potential", which is well-defined in planar geometries and curl-free particle fluxes and which may contribute, in these cases, to formulating a theoretical 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 planar geometry, for simplicity, with z the normal Cartesian direction. We confirm using Brownian Dynamics simulations that µ(z) is spatially constant for active fluids in contact with a soft planar wall, in a gravitational field, and in two-phase coexistence with a planar interface. Next, we show that the coexistence is described by diffusive and mechanical equilibrium with equal bulk pressure and bulk chemical potential of the coexisting phases, provided the swim potential that we introduce in this article, is properly taken into account. However, we conclude that the swim potential and hence the chemical potential µ(z) is not a state function of the density for a macroscopic system.

Methods and Formulation
We consider a three-dimensional dispersion of N active Brownian particles (ABPs) with positions r i = (x i , y i , z i ) and orientationsê i = (sin θ i cos φ i , sin θ i sin φ i , cos θ i ) with polar angle θ 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ê i . The motion of particle i is described by the overdamped Langevin equationṡ 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 0 is the self-propulsion speed. The collisions with the solvent are described by a stochastic force and torque characterised by random vectors Starting from (1) and (2), we average over the noise to derive the deterministic Fokker-Planck equation [9,23] ∂ψ(r,ê, t) ∂t = −∇ · j(r,ê, t) − ∇ê · jê(r,ê, t) for the time evolution of the probability distribution function ψ(r,ê, t) ≡ N i=1 δ(r − r i )δ(ê −ê i ) with . . . the averaging over the random noise. Here we defined the translational and rotational fluxes We introduced here the instantaneous full two-body correlation function ψ (2) , and hence to obtain a closed set of equations one needs a BBGKY-like hierarchy of Fokker-Planck equations for the nbody correlation functions or a mean-field approximation such as ψ (2) The zeroth moment ρ(r, t) = dê ψ(r,ê, t) defines the local particle density, and its time evolution is described by the continuity equation obtained from the zeroth moment of Eq. (3), with the particle flux J(r, t) = dê j(r,ê, t) given by Here ρ (2) (r, r ′ , t) = dê dê ′ ψ (2) (r,ê, r ′ ,ê ′ , t) is the spatial two-body correlation function and the first moment m(r, t) = dê ψ(r,ê, t)ê is the local polarization.
An equation for m(r, t) follows from the first moment of Eq. (3) which yields with d = 2, 3 the spatial dimension of interest, and with the two-rank momentum flux tensor where S = dê ψ(r,ê)(êê − I/d) is the traceless alignment tensor. We now assume that the system is only inhomogeneous in the z-direction, due to either an external potential V e (z) or due to coexistence of two phases separated by an interface parallel to the xy-plane. Without loss of generality, we consider a large, but finite system by setting V e (±∞) = ∞, such that ρ(z → ±∞) = 0. From Eq. (7), we find that the particle flux in the z-direction is given by When divided by βD t , we interpret Eq. (10) as a continuum force balance rather than at the microscopic level, which requires averaging over bins that contain enough colloids for the continuum picture to hold. In the following sections this is achieved by having bins that are very elongated in the direction(s) perpendicular to the z-direction. The term (βD t ) −1 v 0 m z (z, t) has previously been interpreted as a contribution to the divergence of the stress tensor, which has led to a debate on pressure being a state function or not in active systems [23,27,28]. Here, however, we take another point of view, and regard this term as an activity-induced body force that is exerted on the active particles by the solvent [27,29]. This allows us to define the so-called swim potential where V swim (z 0 , 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 value of this constant is determined by surfaces and interfaces, where m can be non-zero, not unlike the Donnan potential in inhomogeneous electrolyte solutions [30,31]. This is a reflection of the fact that the activity-induced body force on the active particles only averages out in the bulk, but not near interfaces.
We now combine Eqs. (10)- (12) to construct, in the spirit of the simplest dynamic density functional theory [32,33] with a density-independent diffusion coefficient, a local chemical potential-like function µ(z, t) by J z (z, t) = −D t ρ(z, t)∂ z βµ(z, t) such that 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 defined by where we have used ρ (2) (10). Eq. (13) reduces to the conventional chemical potential for a passive system, where v 0 = 0, and is constructed such that J z = 0 if µ(z, t) is a spatial constant. The local chemical potential µ(z) is therefore a prime candidate to describe diffusive equilibrium of coexisting phases in stationary states of active systems. Interestingly, all terms in Eq. (13) can be determined in Brownian Dynamics (BD) simulations of ABPs. The body-force interpretation of the polarization (11) can also be used to write the mechanical equilibrium condition of a stationary state in terms of a well-defined normal component 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 finite system, we can rewrite Eq. (10) as with the standard equilibrium-like expression for the (intrinsic) normal pressure where we used Newton's third law and the symmetry of ρ (2) (r, r ′ ) under particle exchange. The last term in Eq. (16) is the virial contribution that describes the zcomponent of the interparticle forces across a plane at z, which can be measured in a BD 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 force balance (15), which turns out to be crucial for interpreting the (osmotic) pressure as a state function [29]. However, in order to connect to existing literature, and for later reference, we do define with the zz-component of J m given by which reduces to the conventional swim pressure [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 a non-negligible role in the force balance obtained from Eq. (15) when significant spatial variations are present, e.g. in the interface of a phase coexistence. To summarize, we have introduced the concept of a swim potential here using a force balance for only the colloids. This force balance can be combined with an additional force balance for the solvent, which provides an alternative interpretation, but identical expression, for the swim 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 dP/dz = −ρ(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), and the swim potential V swim (z) and swim pressure P swim (z) have thus been constructed such that If we now invoke a Local Density Approximation (LDA), i.e. assume that the local environment behaves as a bulk such that the local pressure and chemical potential are a function of only the local density ρ(z), then Eq. (19) can be written in terms of bulk quantities as: allowing us to write dP (ρ) dρ = ρ dµ(ρ) dρ (21) with µ(ρ) = µ int (ρ) + V swim (ρ), in a zero external potential. Here, we shall take care to distinguish 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 in equilibrium (where P swim = V swim = 0) it holds true in general, we emphasize that in this case we had to make use of the LDA to derive it. This Gibbs-Duhem relation provides 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. We test the applicability of Eq. (21) in simulations and show that it works well for cases with low anisotropy (e.g. low polarization). However, Eq. (21) does not hold true in general 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 that is equal to the (incompressible) solvent volume fraction. The equilibrium analogue of Eq. (21) follows naturally if the solvent is treated grand-canonically which we implicitly assume. Both approaches are also similar in the sense that they both identify the fluxes as being proportional to the gradient of a (scalar) chemical potential.
In the next Section, we apply the formalism of Eqs. (12)- (17) to active fluids and consider four different scenarios. We perform Brownian Dynamics (BD) simulations of non-interacting as well as interacting particles in two and three dimensions by employing Eqs. (1) and (2). In Section 3.1 we study a non-interacting active fluid in contact with a short-ranged planar soft wall. We compare and verify that the stationary state is indeed described by constant µ(z) in both the Fokker-Planck calculations and particle based simulations. Next we present the results of BD simulations of an active fluid with Lennard-Jones (LJ) interactions subject to a gravitational field in Section 3.2. In Section 3.3 we consider an active Lennard-Jones fluid exhibiting gas-liquid coexistence with a planar interface and confirm mechanical and diffusive equilibrium. We perform a Maxwell equal-area construction to identify phase coexistence from bulk equations of state. We then attempt to apply the same formalism to active particles which undergo Motility Induced Phase Separation at high activity in Section 3.4.

Active Ideal Gas
We first consider a three-dimensional active ideal gas (with V (r) = 0) at Péclet number Pe = v 0 /σD r = 0 (passive), 1, 3, 5, in the external potential βV e (z) = (z/σ) 2 for z < 0 and V e (z) = 0 for z > 0, where the unit of length σ = 3D t /D r is chosen to be the particle diameter so that the Stokes-Einstein relation for spheres in three dimensions is satisfied. Note that Pe can also be perceived as the ratio of the persistence length v 0 /D r and the particle diameter [5]. For large but finite z = z b 3σ, the active fluid reaches a bulk state with bulk density ρ b = ρ(z b ), and the normal pressure reduces to the bulk pressure P b = P N (z b ) = ρ b k B T . In Fig. 1(a) and (b) we show the timeaveraged density profiles ρ(z) and orientation profiles m z (z)/ρ(z), respectively. We observe that the particles penetrate deeper into the wall at higher Pe resulting into a more extended ρ(z) within the wall accompanied by a small adsorption (that was found in Ref. [37] as well) close to z = 0. In Fig. 1(b) we see no average polarization outside or inside the wall for the passive case. At finite Pe, however, Fig. 1(b) shows that the average orientation is zero in the bulk where V e (z) = 0 and negative within the 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 ∼ 0.1k B T . Clearly, for µ(z) to be constant it is crucial that V swim (z), which is attractive towards the wall consistent with the polarization and extended density profile close to the wall, is included in Eq. (13); ignoring this contribution of 10 − 30k B T would not have yielded a spatially constant chemical potential in the stationary state. Although µ(z) was constructed to be spatially constant within the Fokker-Planck formalism, a confirmation from the simulations serves  as a useful validation. Additionally, we verify that the swim pressure (given by Eq. (17)) measured in the bulk reduces to P swim (z We use this bulk state at z b 3σ with V swim (z 0 = z b ) as the reference point for the profiles of V swim (z) and µ(z) in Fig. 1(c) and (d), respectively.

Sedimenting weakly active LJ-particles
We now consider simulations of weakly active Lennard-Jones (LJ) particles with an isotropic pair potential, V LJ (r) = 4ǫ((σ/r) 12 − (σ/r) 6 ), at k B T /ǫ = 2.0 in the gravitational potential V e (z) = Mgz for z > 0 with a hard 'bottom' at z = 0, with M the buoyant particle mass. These systems are supercritical in the passive case, and therefore even more so in the active cases since the 'critical temperature' decreases with increasing activity [6,17]. We measure the density ρ(z)σ 3 , polarization m z (z)/ρ(z), swim  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 various values of βM gσ, and Péclet number Pe=0 (blue), 10 (green), 20 (red) as obtained from BD simulations. The height z is scaled with respect to l, where l = v 0 /D r is the persistence length for Pe=10 and 20, and l = σ is the particle diameter for Pe=0. The compressibility factor P N (ρ, Pe)/ρ in (e) and the intrinsic chemical potential µ int (ρ, Pe) 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 βMgσ = 0.5 and 1.0 for Pe=0, and βMgσ = 3 and 5 for Pe=10 and 20, all plotted in Fig. 2(a)-(d). In order to obtain a comparable length scale l over which variations are observed in the passive (where we choose l = σ) and in the active cases (where l = v 0 /D r ), we used a smaller buoyant mass of the particles in the passive case. We observe that the polarization m z (z) is positive for Pe=10 and 20, and hence the mean swimming direction is opposite to the gravitational field, consistent with the findings in Ref. [14]. Moreover, Fig. 2(b) shows that the polarization profile m z (z)/ρ(z) is surprisingly constant over a large regime of heights z. As a consequence, the swim potential profile βV swim (z) essentially decreases linearly with height z for Pe = 10 and 20 and counteracts largely the gravitational field, as shown in Fig. 2(c), leading to an enormous increase in sedimentation length (βMg) −1 [10]. The chemical potential profile µ(z) is calibrated by µ(z 0 ) = 0 at the reference point z 0 determined by the condition ρ(z 0 )σ 3 = 5 × 10 −3 . µ(z) is shown in Fig. 2(d) and is indeed spatially constant within our statistical accuracy of ∼ 0.3k B T . It is important to note here that V swim (z) decreases by a few hundred k B T and the external gravitational potential V e (z) = Mgz 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. We observe 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.

Active-LJ phase coexistence
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 with bulk 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éclet number Pe= v 0 /D r σ = 3 are used in this case. In Fig. 3(a), we show a typical configuration of a liquid slab in the center of the simulation box in coexistence with a gas phase on either side. In Fig. 3(b) we plot the corresponding density profile ρ(z) which can be fitted to a hyperbolic tangent function (Eq. (A.1)), independently for z > 0 and 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 in the bulk gas and liquid respectively. In the same figure we also plot the polarization profile m z (z)/ρ(z), showing that the swimming direction of the particles at the liquidgas interface is pointing from the liquid phase towards the gas phase, i.e., against the attractive interparticle forces from the liquid [17,38].
In Fig. 3(c) and (d) we plot the profiles P (z) and µ(z), respectively, which clearly show 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 gas and liquid phase. For completeness, in Fig. 3(c) we also plot the individual contributions to the total pressure P (z) = P N (z) + P swim (z), where P swim (z) is the swim pressure obtained from Eq. (17), and P N (z) = P N,idl (z) + P N,vir (z) is the normal pressure with 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 chemical potential µ(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 using Eq. (12).
In order to investigate if we can predict phase coexistence solely from bulk quantities, we perform BD simulations of bulk states of ABPs at several temperatures k B T /ǫ and Péclet number Pe=2.67. We measure the bulk pressure P as a function of density ρ in a simulation box small enough to prevent phase separation and plot the equations of state P (ρ) for several subcritical temperatures in Fig. 4(a). Now, within a 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 shown in Fig. 4(b). We emphasize here that we refer to µ(ρ) as the µ obtained by applying Eq. (21) which is not to be confused with µ(z). The intersection of the curve µ(P ) gives the coexistence µ g = µ l and P g = P l . In the inset of Fig. 4(b) we compare the binodals in the (scaled) temperature-density plane as obtained from the density profiles from  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 that the 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 effect on the gas-liquid binodals (shown in the inset of Fig. 4(b)) as the critical temperature shifts from k B T /ǫ ≈ 1.15 in the passive case to k B T /ǫ ≈ 0.54 in the active case for Pe=2.67 (see Ref. [17] for full comparison).

Motility Induced Phase Separation
In this section we discuss the swim potential and the chemical potential in a twodimensional system of strongly active particles exhibiting motility induced phase separation at high Pe. We choose our planar geometry in the yz plane and assume homogeneity in the y direction to be consistent with previous definitions. The particles interact with the WCA potential given by V W CA (r) = V LJ (r) + ǫ, with a cut-off beyond r ≥ r c = 2 1/6 σ to make the particles purely repulsive. The particle orientations can be described in terms of a single angle θ i asê i = (cos θ i , sin θ i ). The translational equation of motion in 2D is similar to Eq. (1) and the rotational diffusion followsθ i = √ 2D r Ξ r i , with Ξ r i a zero-mean unit-variance Gaussian random variable. As before, we fix rotational and translation diffusion coefficients to correspond to the particle interaction length scale σ = 3D t /D r and change the self-propulsion speed v 0 to vary Pe. At high Pe, we find that the system phase separates into a gas phase and a dense phase, both of well-defined densities, separated by a planar interface in an elongated simulation box [16]. For Pe=50 the typical density and polarization profiles are shown in Fig. 5(a). Notably, the polarization profiles are now reversed with (c) Total pressure 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 > 0, corresponding to the system described in (a). The inset shows the ideal contribution βµ id (z) = ln ρ(z)σ 2 and that µ(z) is constant within an accuracy of 3k B T . respect to Fig. 3(b) as the particles at the interfaces are now pointing towards the dense phase. We measure the normal component of the total pressure P (z) and the chemical potential µ(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 spatially constant, 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 gas phase 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 for which there is no global phase separation at intermediate densities. We confirm the results of the homogeneous states with larger system sizes and find that the agreement is satisfactory for our analysis. Performing a Maxwell construction on P as a function of 1/ρ gives the equal-area pressure P Maxwell shown as the dashed horizontal line in Fig. 5(b). In the same figure, we also show the coexistence pressure P coex obtained from the direct coexistence simulation of the phases coexisting at the corresponding set of parameters. From the two curves it is evident that the coexistence densities predicted by the Maxwell construction and the direct-coexistence simulations do not agree. We perform the same procedure on a set of Pe in the range 30−60 and plot the corresponding coexistence densities and the densities predicted by the Maxwell construction in the inset of Fig. 5(b). From the disagreement between the two binodals we conclude that the Maxwell equal-area construction does not correspond to the coexisting states as obtained 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 potential introduced in Ref. [12] yields the same binodals as predicted here despite the difference of the factor concerning the solvent volume fraction.

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 predict the coexisting densities ρ g and ρ l [9,26] in systems of ABPs. In other words, even though µ(z g ) = µ(z l ) in a phase-separated system (where z g and z l are locations far from interfaces such that the local densities are ρ g and ρ l , respectively) , the chemical 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 the failure of the LDA assumed in the derivation of Eq. (21), as we will show below. In particular, 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 the interface between the bulk state and the reference state at z 0 . This implies that neither V swim nor µ ex as expressed in Eqs. (12) and (14), respectively, are state functions of the density. Below we show an example for V swim (z) which demonstrates this breakdown of the LDA in the case of a 2D active ideal gas (for which µ ex (z) ≡ 0) in a particular external potential.
The setup consists of a ramp-like external potential βV e (z) = λz/σ in the region 0 < z < 5σ which separates a bulk region at the left (where βV e (z) = 0 for z < 0) from the bulk on the right (where βV e (z) = 5λ for z > 5σ). These external potential are plotted in Fig. 6(a) as dash-dot lines for λ = 0, 0.5, 1, and 2. The probability density ψ(z, θ) is obtained by solving Eq. (3) for V (r) = 0 numerically, at Pe = 1 with a fixed density boundary condition ρσ 2 = 1.0 for z 0 = −10σ 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 The V swim (z) profiles, obtained equivalently from Eq. (22) or from Eq. (12), are plotted as solid lines in Fig. 6(c) where we have taken z 0 = −10σ as the reference state where V swim (z 0 ) = 0. If we would approximate the vicinity of any point z ′ as an isotropic bulk with density ρ(z ′ ) in evaluating the swim potential V swim (z ′ ), i.e. assume in Eq. (22) S zz (z ′ ) ≈ m z (z ′ ) ≈ 0 such that the term in square brackets vanishes for every z ′ , we obtain βV LDA swim (ρ(z)) = (v 2 0 /2D t D r ) ln ρ(z)σ 2 which we refer to as the local density approximation (LDA) of Eq. (22). Note that Eq. (22) follows from the Fokker-Planck formalism, 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 from the 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, we can conclude that the values for V LDA swim obtained from the Gibbs-Duhem equation are not correct in general. This is due to the failure of LDA, i.e. due to the anisotropy in the interface that renders the integral on the right hand side in Eq. (22) non-negligible as compared to the first term. In Fig. 6(b) we see that the polarization within the interface increases with λ, consistent with this idea of increasing anistropy. For an interacting system the forces between particles would add another contribution to m z (z)/ρ(z), which could also become a source of failure for the LDA. This will be studied in more detail in a future publication.
In Section 3.3 we observed that the Maxwell construction was able to predict the coexistence densities for the active LJ case with reasonable accuracy, but was in disagreement at higher activity in Section 3.4 for MIPS. We now assert that the error made in the chemical potential by assuming the LDA translates into an error in the predicted coexisting densities that is small for the active LJ particles, but significant for MIPS. We define the error in predicted coexistence densities of the gas and the dense phase, respectively, as ∆ρ err g = ρ(z g ) − ρ g and ∆ρ err l = ρ(z l ) − ρ l , where ρ(z g ) and ρ(z l ) are the bulk coexistence densities and ρ g and ρ l denote the estimates obtained by performing a Maxwell construction. If we define the gas state as the reference state for the chemical potential, i.e. µ(z 0 ) = µ(ρ g ) in Eq. (13) with z 0 = z g , then the error made in determining the chemical potential of the dense phase by using the Gibbs-Duhem equation (21) is ∆µ err l = µ(ρ l ) − µ(z l ), where we recall that µ(ρ l ) is the chemical potential of the dense phase obtained from the Gibbs-Duhem relation, whereas µ(z l ) is the true chemical potential determined in the coexistence simulation. From ∆µ err l the relative error in the predicted density of the dense phase can be estimated as ∆ρ err l /ρ(z l ) ≈ 1/ρ(z l ) · ∆µ err l /(dµ/dρ) l . Similarly, the error in the predicted density of the gas phase can be estimated by using the dense phase as the reference state (µ(z 0 ) = µ(ρ l )). The relative density error estimated in this manner is less than 5% for the active LJ case, whereas it is of the order of 100% for the MIPS case, which agrees with our findings in Fig. 4(b) and 5(b), respectively.
We wish to make a note that the anisotropy terms identified here resemble the interfacial contributions discussed in Ref. [26] for pairwise-interacting particles. Although it requires explicit measurement of these interfacial contributions by performing phase-coexistence simulations, Solon et al. were able to suggest a modified Maxwell construction for estimating the binodals in Ref. [26].
Moreover, our elongated simulation box in Section 3.3 and 3.4 forces the system to 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. In other geometries the stationary state condition ∇ · J = 0 still allows for swirls that correspond to non-zero ∇ × J, for which our expressions for mechanical and diffusive equilibrium break down and a whole new framework is needed. Furthermore, the regime of applicability of Eq. (13) is limited by the underlying dynamic DFT relation, where a ρ-independent diffusion coefficient D t is assumed; an extension to account for a ρdependent diffusion coefficient is left for a future study.

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) defined by Eq. (12) in addition to ideal, excess, and external contributions well-known from equilibrium. Our BD simulations confirm that µ(z) is spatially constant in steady states of several inhomogeneous ideal and interacting fluids of active particles, with V swim (z) an important contribution that counteracts either the external potential V e (z) or the excess contribution µ ex (z). In the low activity regime studied for active LJ fluid, the chemical potential 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 be invalid, which provides support to the conclusions of Ref. [26] that the details of the interface are necessary to determine the coexisting bulk densities. Our formalism opens new avenues towards a Fokker-Planck and dynamic density functional description (of stationary states) of active systems, especially for planar geometries.
of the number of particles n(z) in the slabs of volume L 2 ∆z (ρ(z) = n(z) /L∆z in 2D) arranged parallel to xy plane (y-direction in 2D), where L is the length of the system in the x and/or y-direction, and where ∆z = 0.1σ is the width of the slab. In a similar manner we measure the polarization profile m z (z) by summing the particle orientations in a slab at location z. The density profiles ρ(z) can be fitted to a hyperbolic tangent function given by: ρ(z) = 1 2 (ρ(z l ) + ρ(z g )) + 1 2 (ρ(z l ) − ρ(z g )) tanh 2(z − z * 0 ) D , where ρ(z l ) and ρ(z g ) are the corresponding bulk liquid and vapour coexisting densities, z * 0 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 where we use ρ(z) and m z (z) as measured in the BD simulations, and where V swim (z 0 ) is a suitably chosen reference state. In addition, we measure the normal component of the 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: 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 and the slab of width ∆z centered at z. The integral in Eq. (A.5) denotes that the virial contribution to the pressure of particle pair i and j is due to the part of r ij that lies inside the respective slab at z within the coarse-grained Irving-Kirkwood approximation [34]. We also calculate the swim pressure with the excess chemical potential µ ex (z) defined as (A.8) Here, the excess chemical potential at z with respect to a reference at z 0 is determined by integrating the averaged force that a particle feels due to the particle interactions with all other particles in the system over the distance z 0 to z. Alternatively, if V e (z) = 0, µ(z) can also be obtained using with P (z) = P N (z) + P swim (z).