Lack of an equation of state for the nonequilibrium chemical potential of gases of active particles in contact

TL;DR

This paper investigates the concept of nonequilibrium chemical potential in active particle gases, revealing it lacks an equation of state and depends on potential details, challenging equilibrium thermodynamics assumptions.

Contribution

It demonstrates that the nonequilibrium chemical potential for active particles does not follow an equation of state and depends on potential barrier details, unlike equilibrium systems.

Findings

Chemical potential depends on potential barrier shape

No equation of state for the nonequilibrium chemical potential

Maxwell relation is invalid in this context

Abstract

We discuss the notion of nonequilibrium chemical potential in gases of non-interacting active particles filling two compartments separated by a potential energy barrier. Different types of active particles are considered: run-and-tumble particles, active Brownian particles, and active Brownian particles with a stochastic reorientation along an external field. After recalling some analytical results for run-and-rumble particles in one dimension, we focus on the two-dimensional case and obtain a perturbative expression of the density profile in the limit of a fast reorientation dynamics, for the three models of active particles mentioned above. Computing the chemical potentials of the non-equilibrium systems in contact from the knowledge of the stationary probability distribution of the whole system ---which agrees with a recently proposed general definition of the chemical potential in non-equilibrium systems in contact--- we generically find that the chemical potential lacks an equation of state, in the sense that it depends on the detailed shape of the potential energy barrier separating the compartments and not only on bulk properties, at odds with equilibrium. This situation is reminiscent of the properties of the mechanical pressure in active systems. We also argue that the Maxwell relation is no longer valid and cannot be used to infer the nonequilibrium chemical potential from the knowledge of the mechanical pressure.

Figures

Equations

$$\dot{\mathbf{r}}=v_0(\mathbf{r})\mathbf{e}(\theta )-\eta \nabla U,$$

$$\dot{\theta }=\xi (t)$$

$$⟨\xi (t)\xi (t^{\prime })⟩=2D_{\mathrm{r}}\delta (t-t^{\prime }),$$

$$P(N_A|N)\sim e^{-V\mathcal{I}({\rho }_A,{\rho }_B)}$$

$$\mathcal{I}({\rho }_A,{\rho }_B)={\gamma }_A{I}_A({\rho }_A)+{\gamma }_B{I}_B({\rho }_B),$$

$$I_A^{\prime }({\rho }_A)-I_B^{\prime }({\rho }_B)={\mu }_A^{\mathrm{cont}}-{\mu }_B^{\mathrm{cont}}$$

$$\frac{\partial P_R}{\partial t}(x,t)=$$

$$-\alpha \left(P_R(x,t)-P_L(x,t)\right)$$

$$\frac{\partial P_L}{\partial t}(x,t)=$$

$$+\alpha \left(P_R(x,t)-P_L(x,t)\right).$$

$$\frac{\partial P}{\partial t}(x,t)$$

$$\frac{\partial M}{\partial t}(x,t)$$

$$-2\alpha M(x,t),$$

$$P(x)\propto \frac{1}{1-{\left(\frac{\eta U^{\prime }(x)}{v_0}\right)}^2}$$

$$\times \exp \left[-{\int }^x\mathrm{d}q\left(\frac{\eta \alpha }{v_0^2}\right)\frac{U^{\prime }(q)}{1-{\left(\frac{\eta U^{\prime }(q)}{v_0(q)}\right)}^2}\right],$$

$$\rho (x)$$

$$\times \exp \left[-{\int }_{x_k^{\ast }}^x\mathrm{d}q\frac{\eta \alpha }{v_k^2}\frac{U^{\prime }(q)}{1-{\left(\frac{\eta U^{\prime }(q)}{v_k}\right)}^2}\right],$$

$$v_A{\rho }_A^{\ast }{e}^{-\Delta Q_A}=v_B{\rho }_B^{\ast }{e}^{-\Delta Q_B}$$

$$\Delta Q_k={\int }_{x_k^{\ast }}^0\mathrm{d}q\frac{\alpha \eta }{v_k^2}\frac{U^{\prime }(q)}{1-(\eta U^{\prime }(q)/v_k{)}^2}.$$

$$N={\int }_{-L_A}^{L_B}\mathrm{d}x\rho (x)=L_A{\rho }_A^{\ast }+L_B{\rho }_B^{\ast }$$

$$P(N_A|N)=(\genfrac{}{}{0ex}{}{N}{N_A})p_A^{N_A}{p}_B^{N-N_A}$$

$$p_A$$

$$p_B$$

$$P({\rho }_A|\bar{\rho })\sim e^{-L\mathcal{I}({\rho }_A,{\rho }_B)}$$

$$\mathcal{I}({\rho }_A,{\rho }_B)={\gamma }_A{\rho }_A\ln \frac{{\rho }_A}{{\rho }_A^{\ast }}+{\gamma }_B{\rho }_B\ln \frac{{\rho }_B}{{\rho }_B^{\ast }}$$

$${\gamma }_A{\rho }_A+{\gamma }_B{\rho }_B=\bar{\rho },$$

$$I_k({\rho }_k)={\rho }_k\ln \frac{{\rho }_k}{{\rho }_k^{\ast }}(k=A,B).$$

$$\frac{1}{{\gamma }_A}\frac{d\mathcal{I}}{d{\rho }_A}=I_A^{\prime }({\rho }_A)-I_B^{\prime }({\rho }_B)=\ln \frac{{\rho }_A}{{\rho }_B}-\ln \frac{{\rho }_A^{\ast }}{{\rho }_B^{\ast }}.$$

$$I_A^{\prime }({\rho }_A)-I_B^{\prime }({\rho }_B)=\ln \frac{{\rho }_A{v}_A{e}^{-\Delta Q_A}}{{\rho }_B{v}_B{e}^{-\Delta Q_B}}.$$

$$I_A^{\prime }({\rho }_A)-I_B^{\prime }({\rho }_B)={\mu }_A^{\mathrm{cont}}({\rho }_A)-{\mu }_B^{\mathrm{cont}}({\rho }_B),$$

$${\mu }_k^{\mathrm{cont}}({\rho }_k)=\ln \frac{{\rho }_k{v}_k}{\alpha }-\Delta Q_k.$$

$${\mu }_A^{\mathrm{cont}}({\rho }_A^{\ast })={\mu }_B^{\mathrm{cont}}({\rho }_B^{\ast })$$

$$\frac{\mathrm{d}P}{\mathrm{d}t}({\rho }_A|\bar{\rho })$$

$$=v_A\rho (0^{-})\left(P({\rho }_A+\frac{1}{V_A}|\bar{\rho })-P({\rho }_A|\bar{\rho })\right)$$

$$+v_B\rho (0^{+})\left(P({\rho }_A-\frac{1}{V_A}|\bar{\rho })-P({\rho }_A|\bar{\rho })\right),$$

Full text

Self-propelled particles, Active matter, Chemical potentials, Large deviations

I Introduction

Active matter systems like assemblies of macroscopic particles in which energy is continuously fed at the particle scale Marchetti et al. (2013), are a useful benchmark to test nonequilibrium extensions of thermodynamic concepts Speck (2016); Takatori and Brady (2015); Solon et al. (2015a, 2018a). These systems include both artificial particles like self-propelled colloids Palacci et al. (2010); Theurkauff et al. (2012); Buttinoni et al. (2013); Palacci et al. (2013); Bricard et al. (2013); Gomez-Solano, Blokhuis, and Bechinger (2016) or vibrated macroscopic particles Deseigne, Dauchot, and Chaté (2010); Narayan, Menon, and Ramaswamy (2006); Kudrolli et al. (2008), and biological systems like assemblies of bacteria for instance Peruani et al. (2012); Zhou et al. (2014). On the theoretical side, several simple models of non-interacting active particles have been introduced (see Ref. Bechinger et al., 2016; Fodor and Marchetti, 2018 for recent reviews). Two of these models have emerged as paradigmatic models of active particles, namely the Run-and-Tumble Particle (RTP) model Tailleur and Cates (2008); Cates and Tailleur (2013); Schnitzer (1993) modeling the motion of bacteria like E. Coli, and the Active Brownian Particle (ABP) Cates and Tailleur (2013); Stenhammar et al. (2014) model describing active colloids for instance.

In recent years, it has been realized that even an apparently well-defined notion like mechanical pressure can exhibit an unexpected behavior, as it lacks an equation of state Fily et al. (2017); Solon et al. (2018a, 2015a, 2015b); Speck and Jack (2016); Takatori, Yan, and Brady (2014); Winkler, Wysocki, and Gompper (2015); Joyeux and Bertin (2016). This means that unlike at equilibrium, the mechanical pressure exerted on a wall depends not only on the bulk density, but also on the detailed shape of the diverging potential energy profile defining the wall. This dependence results from a nonequilibrium density profile within the (soft) wall, and can be ultimately traced back to the lack of momentum conservation in the bulk Fily et al. (2017), which is a consequence of the self-propulsion mechanism.

At equilibrium, intensive parameters (like temperature, pressure, chemical potentials,…) are generally defined as derivatives of the free energy with respect to their associated extensive variable (resp. energy, volume, number of particles,…). However, out of equilibrium, the free energy is in general no longer defined and the natural thermodynamic definition can no longer be applied. There are essentially two ways to circumvent this issue. On the one hand, one can apply operational definitions. For instance, pressure can always be considered as the average mechanical force per unit area on the surface of the container. On the other hand, one can try to generalize the thermodynamic definition which relies on the equalization of intensive parameters when two systems in contact are allowed to exchange a conserved quantity, e.g., energy for temperatureBertin, Dauchot, and Droz (2004); Shokef, Shulkind, and Levine (2007) or number of particules for chemical potentialsBertin, Dauchot, and Droz (2006); Bertin et al. (2007). As energy is generically not conserved in driven systems, there is no clear way to extend the validity of the equilibrium definition along this line. In contrast, the number of particles is often conserved even for driven systems, offering a way to define a nonequilibrium chemical potential (whose operational definition in terms of work is less straightforward than for the pressure). Based on a large deviation approach, a precise definition of such a nonequilibrium chemical potential has been proposed recently for generic driven systems in weak contact.Guioth and Bertin (2018); Bertin, Dauchot, and Droz (2006); Bertin et al. (2007). This well-grounded definition notably helps clarifying some of the results reported in the context of lattice gas models in contact.Guioth and Bertin (2018); Pradhan, Ramsperger, and Seifert (2011); Dickman and Motai (2014) The main result is that the nonequilibrium chemical potential defined in this way equalizes by construction between the two systems in contact, but generically depends on the details of the contact dynamics.Guioth and Bertin (2018)

Here, we wish to explore the consequences of this definition of a nonequilibrium chemical potential in the more realistic context of systems of active particles, for which experimental realizations are availablePalacci et al. (2010); Theurkauff et al. (2012); Buttinoni et al. (2013); Palacci et al. (2013); Bricard et al. (2013); Gomez-Solano, Blokhuis, and Bechinger (2016); Deseigne, Dauchot, and Chaté (2010); Narayan, Menon, and Ramaswamy (2006); Kudrolli et al. (2008) (see also the very recent work [Chakraborti and Pradhan, ] for a related attempt in the Vicsek model). In the spirit of the seminal work of Sasa and Tasaki on the generalization of thermodynamics to nonequilibrium systems,Sasa and Tasaki (2006) a weak contact can be realized between two systems by connecting them through a high potential energy barrier. Such a setting can be easily implemented for active particles. With this goal in mind, we consider three distinct models of non-interactive active particles. The first two models are standard models, namely the Run-and-Tumble Particle (RTP) model, and the Active Brownian Particle (ABP) model Cates and Tailleur (2013). In both cases, particles have an overdamped dynamics with a velocity characterized, in two dimensions, by an angle $\theta(t)$. In the presence of an external potential $U({\bf r})$, the position ${\bf r}(t)$ of a particle evolves according to

[TABLE]

with ${\bf e}(\theta)$ the unit vector in the direction $\theta$, $v_{0}(\mathbf{r})$ the (possibly space-dependent) speed of the particles, and $\eta$ is the translational mobility coefficient. The dynamics of the angle $\theta(t)$ depends on the model. For the RTP dynamics, a new value of the angle $\theta$ is drawn from a uniform distribution over $[0,2\pi)$ with a probability rate $\alpha$. The run-and-tumble moded can also be defined in one dimension, in which case the angle $\theta$ can only take the values [math] and $\pi$. In the case of ABPs, the angle $\theta$ continuously diffuses,

[TABLE]

where $\xi(t)$ is a white noise satisfying $\langle\xi(t)\rangle=0$ and

[TABLE]

with $D_{\rm r}$ the angular diffusion coefficient. As the angle $\theta$ varies continuously, no one-dimensional version of the ABP model can be defined.

As a third model, we generalize the ABP dynamics by adding to the diffusive dynamics of $\theta$ given by Eq. (2) a stochastic jump dynamics that tends to align the velocity along an externally given direction. This model is inspired by experiments on microalgae Chlamydomonas Reinardtii Martin et al. (2016), that randomly reorient their swimming direction towards the direction $\theta_{0}$ of a light source, with some noise in the reorientation process. The model is thus defined as an ABP model of angular diffusion coefficient $D_{\rm r}$, with a random reorientation, at rate $\lambda$, to a new direction $\theta^{\prime}$ drawn from a distribution $\psi(\theta-\theta_{0})$.

The definition of a nonequilibrium chemical potential proposed in Refs. Bertin, Dauchot, and Droz, 2006; Bertin et al., 2007; Guioth and Bertin, 2018 relies on a large deviation formalism, that we briefly recall here. Splitting a given system into two subsystems A and B such that the total number of particles $N_{A}+N_{B}=N$ is conserved (while particles can be exchanged between A and B), one determines the distribution $P(N_{A}|N)$ under the asymptotic form, valid for a large volume $V=V_{A}+V_{B}$,

[TABLE]

where $\rho_{A}=N_{A}/V_{A}$ and $\rho_{B}=N_{B}/V_{B}$ are constrained by the conservation law $\gamma_{A}\rho_{A}+\gamma_{B}\rho_{B}=\overline{\rho}$, with $\gamma_{k}=V_{k}/V$ and $\overline{\rho}$ the overall density. The function $\mathcal{I}(\rho_{A},\rho_{B})$ is called a large deviation function. In the weak contact limit, when the rate of particle exchange is very small (for instance due to a high energy barrier separating the two systems), and when a single particle is exchanged at a time, it has been shownGuioth and Bertin (2018) that the large deviation function $\mathcal{I}(\rho_{A},\rho_{B})$ is additive,

[TABLE]

provided that the contact dynamics satisfies a factorization property with respect to compartments A and B. When the additivity property of the large deviation function holds, nonequilibrium chemical potentials can be defined from the relation

[TABLE]

for systems A and B in contact (the prime denotes a derivative). The goal of this paper is thus to determine the large deviation function $\mathcal{I}(\rho_{A},\rho_{B})$ in the different models of active particles mentioned above, and to deduce from it the chemical potentials $\mu_{k}(\rho_{k})$, $k=A$, $B$. Note that in the approach presented here, the evaluation of $\mathcal{I}(\rho_{A},\rho_{B})$ requires first to determine the stationary density profile $\rho(x)$. As a result, the obtained chemical potentials are not used to predict the stationary densities in systems A and B in contact, as these stationary densities are already known from the density profile (yet, equalizing the chemical potentials yields back the stationary densities). The focus of our study is rather to determine whether the nonequilibrium chemical potentials that we have defined obey an equation of state, in the case of the simplest systems of active particles.

On the technical side, the lack of time-reversibility (or detailed balance) in the dynamics —which appears in the presence of a non-uniform potential as we will see henceforth— makes it difficult to compute the stationary density profile, even for independent particles. However, perturbative computations may be feasible in certain limits of the parameters. This is the case in particular when the reorientation time scale is very small compared to the translational time scales as we will see in Sect. III for the three models defined above. Before presenting this generic perturbative method, we recall in Sect. II the exact results obtained for the one-dimensional run-and-tumble model, and discuss how the potential energy barrier unveils the otherwise hidden irreversibility of the dynamics Fodor et al. (2016); Nardini et al. (2017).

II Run-and-Tumble particles in 1D

In one dimension, there are only two possible directions of motion, namely left or right. Calling $P_{L}(x,t)$ and $P_{R}(x,t)$ the probability density distributions associated with a velocity oriented to the left or to the right respectively, the master equation reads

[TABLE]

Defining $P(x,t)=P_{R}(x,t)+P_{L}(x,t)$ and $M(x,t)=P_{R}(x,t)-P_{L}(x,t)$, equations (7) become

[TABLE]

with $J(x,t)=v_{0}(x)M(x,t)-\eta U^{\prime}(x)$.

Now, let us formally split our system into two parts A and B that are separated by a high barrier $U(x)$. In this one-dimensional context, we define compartment A as the interval $(-L_{A},0)$ and compartment B as the interval $(0,L_{B})$. One can also assume that $v_{0}(x)=v_{A}$ if $-L_{A}<x<0$ and $v_{0}(x)=v_{B}$ if $0<x<L_{B}$, i.e., that the velocities in A and B can be different. This situation may be particularly relevant for light-controlled self-catalytic propulsion as described for instance in Ref. Buttinoni et al., 2013; Palacci et al., 2013; Stenhammar et al., 2016; Gomez-Solano, Blokhuis, and Bechinger, 2016, where the light intensity could differ in compartments A and B, leading to different self-propulsion speeds (note however, that ABP are a better description of active colloids that RTPCates and Tailleur (2013)).

In the bulk of each system, where the potential $U(x)$ is flat, the current reads $J(x,t)=v_{A,B}M(x,t)$. Since the boundaries are closed, the current is vanishing at the boundaries. The question we ask is how the two densities $\rho_{A}$ and $\rho_{B}$ in the bulk of each system are related to each other when particles are allowed to cross the potential barrier $U(x)$ centered at $x=0$.

To address this question, we look for a solution in the presence of the potential barrier. For the self-propulsion to be strong enough so that a particle can overcome the barrier, the condition $v_{0}(x)>\eta|U^{\prime}(x)|$ has to be fulfilled for all $x$. If this is the case, it has been shown Tailleur and Cates (2008) that a steady-state solution of equations (8) reads, on an interval where $v_{0}(x)$ takes a constant value $v_{0}$,

[TABLE]

with $P(x)=0$ if $v_{0}(x)<\eta|U^{\prime}(x)|$. One can see that this solution is associated with a vanishing current $J(x)$. Note that an accumulation of particles occurs around the potential energy barrier (see, e.g., Ref. Fodor and Marchetti, 2018), which already suggests that out-of-equilibrium effects are present near the separating barrier.

In order to compute the stationary distribution of density in $A$ and $B$, one can proceed in two ways. The first one is straightforward: as one knows the exact stationary distribution for the whole system, one can compute exactly the probability density of a configuration $\{X_{i}\}_{i=1}^{N}$ and then find the probability to have $N_{A}$ particles in $A$. In a word, this way is the one someone would apply in an equilibrium situation where the Gibbs-Maxwell-Boltzmann distribution is known. The second one proceeds in a similar fashion as Ref. Guioth and Bertin, 2018 and focuses directly on the exchange dynamics at the contact area.

II.1 Evaluation of the chemical potentials from the density profile

In the case described above of two compartments with speed $v_{A}$ and $v_{B}$, the stationary density field $\rho(x)=NP(x)$ reads for $x\in\Lambda_{k}$ ($k=A,\,B$), using Eq. (9)

[TABLE]

with $x_{k}^{\ast}$ a point in the bulk of $k$ where the potential barrier $U(x)=0$, and $\rho_{k}^{\ast}$, the densities in the bulks to be determined. The current being continuous at contact between $A$ and $B$, the constants $\rho_{A}^{\ast}$ and $\rho_{B}^{\ast}$ are not independent from each other but are related by the equality of the stationary fluxes of particles going from $A$ to $B$ and from $B$ to $A$, yielding

[TABLE]

with

[TABLE]

At the thermodynamic limit where $L_{A}$ and $L_{B}$ tend to infinity, the normalization condition of $\rho(x)$ reads

[TABLE]

since the non-uniformity caused by the external potential $U(x)$ stays localized around the contact area. Equations (11) and (13) of course completely determine $\rho_{A}^{\ast}$ and $\rho_{B}^{\ast}$ but one would like to see how chemical potentials of systems in contact can be defined in this simple case.

Knowing the density field $\rho(x)$, one can now calculate the probability $P(N_{A}|N)$ to observe $N_{A}$ particles in compartment A (and thus $N_{B}=N-N_{A}$ particles in B). Particles being independent, the corresponding probability distribution is simply a binomial

[TABLE]

where $p_{A}$ and $p_{B}$ are the probabilities for a given particle to be in compartment $A$ or $B$ respectively. These probabilities are simply expressed in terms of the density field $\rho(x)$ as

[TABLE]

In the thermodynamic limit, the distribution $P(N_{A}|N)$ takes the large deviation form

[TABLE]

with $L=L_{A}+L_{B}$, and a large deviation function

[TABLE]

where $\rho_{A}$ and $\rho_{B}$ are constrained by

[TABLE]

with $\gamma_{k}=L_{k}/L$. Using Eq. (11), the large deviation function $\mathcal{I}(\rho_{A},\rho_{B})$ can be rewritten in the additive form Eq. (5), with $I_{k}(\rho_{k})$ given by

[TABLE]

Taking the derivative of $\mathcal{I}(\rho_{A},\rho_{B})$ with respect to $\rho_{A}$ under the contraint (19), one obtains

[TABLE]

Using Eq. (11) to reexpress the ratio $\rho_{A}^{\ast}/\rho_{B}^{\ast}$, one finds

[TABLE]

Since one has by definition

[TABLE]

the chemical potential $\mu_{k}^{\mathrm{cont}}(\rho_{k})$ can be expressed as

[TABLE]

In addition, the stationary densities $\rho_{k}^{*}$ obey

[TABLE]

thanks to the property $d\mathcal{I}/d\rho_{A}(\rho_{A}^{\ast},\rho_{B}^{\ast})=0$.

The parameter $\alpha$ as been introduced in Eq. (24) as a typical rate (inverse time scale) to get a dimensionless argument in the logarithm; yet, the choice of this inverse time scale is arbitrary. One notices that beyond the expected “perfect gas” contribution $\ln\rho_{k}$ in Eq. (24), a dependence on the particle speed appears, and more importantly a non-trivial dependence on the contact appears in the expression of a chemical potential at contact, through the term $\Delta Q_{k}$. Even when the velocities on each side are equals ($v_{A}=v_{B}$), $\Delta Q_{A}$ and $\Delta Q_{B}$ (see equation (12)) are different as soon as the potential barrier $U(x)$ is asymmetric. As a result, the chemical potential does not satisfy an equation of state, in the sense that it does not depend only on the bulk density, but also on the details of the contact dynamics.

II.2 Computation based on the Poisson dynamics over the particle number

As an alternative determination of the chemical potential, we would like to compute the large deviation function $\mathcal{I}(\rho_{A},\rho_{B})$ for the one-dimensional run-and-tumble particle model, by directly using a model of the exchange dynamics. The independence of the particles implies that the number of particles $N_{A}$ in system $A$ obeys a Poisson process. The associated master equation reads

[TABLE]

with

[TABLE]

$\rho_{k}$ being the actual densities in the bulk of system $k$ such that $\gamma_{A}\rho_{A}+\gamma_{B}\rho_{B}=\bar{\rho}$, and $\Delta Q_{k}$ being defined in equation (12). Note that $v_{k}\rho_{k}e^{-\Delta Q_{k}}$ is nothing but the average flux of particles (and also the transition rates here for this Poisson process) directed from $k\to k^{\prime}\neq k$, $k,k^{\prime}=A,B$.

Using the large deviation form (17) in Eq. (26) yields an equation on the derivative of the large deviation function associated with $P(\rho_{A}|\bar{\rho})$. Solving this equation, one finds

[TABLE]

From the relation $\mathcal{I}^{\prime}(\rho_{A},\rho_{B})=\mu_{A}^{\mathrm{cont}}-\mu_{B}^{\mathrm{cont}}$, we thus recover, up to the introduction of the (arbitrary) normalizing time scale $\alpha$, the chemical potentials $\mu_{k}^{\mathrm{cont}}$ introduced in Eq. (24).

As emphasized in Ref. Guioth and Bertin, 2018, we note that the additivity of the large deviation function is a consequence of the factorization of the macroscopic transition rates (which are here of the so-called Sasa-Tasaki typeSasa and Tasaki (2006)) as well as macroscopic detailed balance, which results from the exchange of a single particle at a time between compartments A and B.

II.3 Entropy production of 1D RTPs

Entropy production in systems of active particles is a natural thermodynamic quantity to look at, because active particles dissipate energy and are thus expected not to obey time-reversal symmetryFodor et al. (2016); Nardini et al. (2017); Pietzonka and Seifert (2018). With this idea in mind, it is interesting to note that from the point of view of trajectories, the run-and-tumble dynamics in the absence of external potential is time-reversible (however, the underlying processes that generate the self-propulsion are irreversible and do generate entropy production Pietzonka and Seifert (2018), but they are not described explicitly in the run-and-tumble model). One may thus wonder why some nonequilibrium effects would come into play. As seen above, nonequilibrium effects occur around the energy barrier, as the shape of the energy barrier matters for the definition of the chemical potentials of the systems in contact. One thus expects that the presence of an external potential may break time-reversibility, and should thus lead to a non-zero entropy production. Let us then compare trajectories with their time-reversal counterpart to examine the non-equilibrium character of this model. For one particle, a trajectory of length $T$ corresponds to the position $x(t)$ as well as the direction of its velocity $\pm v_{0}$, or equivalently the times of velocity reversals $\{t_{k}\}_{k=1}^{\mathcal{N}_{j}}$, over the interval $[0,T]$, $\mathcal{N}_{j}$ being the number of jumps during this interval. The position $x(t)$ thus obeys the equation of motion

[TABLE]

with

[TABLE]

Calling $x^{\mathrm{R}}(t)=x(T-t)$ the time-reversed trajectory, it necessarily obeys the equation

[TABLE]

The question is to know whether this time-reversed trajectory obeys the same stochastic dynamics as the original one. Clearly, the last term cannot be realised by the natural noise present in the system since the latter is only bi-valued, unless the external force vanishes ($U^{\prime}(x)=0$). Thus, the presence of an external potential breaks the time-reversal symmetry and the entropy production gets infinite, because the probability to observe the time-reversed trajectory is simply zero. This infinite value of the entropy production is an artifact of the absence of thermal diffusion. In principle, a small thermal diffusive motion would come on top of the ballistic motion described by Eq. (1). Taking into account this thermal contribution restores a finite entropy production.

III Active dynamics in 2D: perturbative evaluation

As mentioned in the introduction of Sec. II, very few situations even for independent particles allow explicit exact calculations. To see if the effect of the potential described before is not specific to the one-dimensional run-and-tumble dynamics, we now study several two-dimensional models of active particles, namely the two-dimensional run-and-tumble particle (RTP) model, the active Brownian particle (ABP) model, and the ABP with an external field. We generically call $\Lambda_{A}$ and $\Lambda_{B}$ the regions of space defining the compartments A and B.

III.1 Run-and-Tumble particles in 2D

In two dimensions, relevant variables are the position $\mathbf{r}=x\bm{e}_{x}+y\bm{e}_{y}$ and the angle $\theta\in[0,2\pi)$. The master equation which describes the evolution of the probability density function $P(\mathbf{r},\theta,t)$ reads

[TABLE]

Finding the stationary solution of Eq. (31) is difficult. However, one expects an equilibrium-like behaviour in the limit where an infinite tumbling rate would lead to a translational diffusive-like behaviour (if the velocity $v_{0}$ scales in the proper way). For the model to remain well defined, the effective diffusion coefficient $D\sim v_{0}^{2}/\alpha$ should remain finite when $\alpha\to\infty$. Our goal is then to compute the first order correction in $\alpha$ to this equilibrium-like limit. For the sake of simplicity and in order to see only the genuine out-of-equilibrium effect of an asymmetry of the potentiel barrier, we will assume here that $v_{0}(\bm{r})$ is uniform all along the whole system, $v_{0}(\bm{r})=v_{0}$.

To perform this study, we introduce several quantities and notations (see Ref. Peshkov et al., 2014 for a review of these techniques). First, we define the angular Fourier transform

[TABLE]

with $f_{k}^{*}=f_{-k}$, $z^{*}$ being the complex conjugate of $z\in\mathbb{C}$. With these notations, one gets

[TABLE]

Projecting the master equation (31) onto the Fourier basis we have just introduced, gives a hierarchy of equations on $f_{k}$:

[TABLE]

where one has introduced the complex derivatives

[TABLE]

In order to perform the perturbative expansion, let us write the first three equations:

[TABLE]

For $k\geqslant 1$, it is not difficult to see that $f_{k}$ has a fast relaxation dynamics for large $\alpha$, so that the time derivative $\partial_{t}f_{k}$ can be neglected on time scales large compared to $\alpha^{-1}$. Under these assumptions, we obtain

[TABLE]

Fixing the large scale diffusion coefficient $D=v_{0}^{2}/2\alpha$, one obtains that $v_{0}\sim\alpha^{1/2}$. A careful look at Eq. (III.1) indicates that a consistent scaling is given by $f_{k}\sim\alpha^{-|k|/2}$. Then, by inspection, one sees from Eqs. (36) that to expand the continuity equation (36a) to order $\alpha^{-1}$, one needs to expand $f_{1}$ to order $\alpha^{-3/2}$, and $f_{2}$ to order $\alpha^{-1}$. At these orders, $f_{1}$ and $f_{2}$ read

[TABLE]

where $\Delta$ is the Laplacian. One eventually gets a closed equation on the density field $\rho(\mathbf{r})$ which reads

[TABLE]

Looking for a stationary perturbative solution, we propose the following ansatz

[TABLE]

where $\phi_{0}$ and $\phi_{1}$ do not depend on $\alpha$. Eq. (39) can be rewritten in terms of the current ${\bf J}$ of particles, as $\partial_{t}\rho=\bm{\nabla}\cdot{\bf J}$. We now assume that $U(\bm{r})=U(x)$, i.e., that $U$ does not depend on $y$. In this case, the stationary current vanishes, and one finds

[TABLE]

Note that $D/\eta$ plays the role of an effective temperature. At equilibrium, this ratio would correspond to the true thermodynamic temperature, from the fluctuation-dissipation theorem.

More importantly, one notices the presence of a non-local contribution of the external potential barrier at order $\alpha^{-1}$, because the term $\int_{0}^{x}\mathrm{d}q\left(U^{\prime}(q)\right)^{3}$ does not vanish in the bulk, contrary to the other terms that involve only local derivatives of the potential. This cubic contribution is quite similar to the 1D case at order $\alpha^{-1}$ and gives different contributions to $\Delta Q_{k}$ in (12) depending on $k=A,\,B$. As a result, the stationary densities $\rho_{A}$ and $\rho_{B}$ depend on the detailed shape of the potential energy barrier.

III.2 Active Brownian Particles in 2D

For the active Brownian particles, the dynamics is very similar except for the angle that follows a Brownian motion rather than a jump process. The corresponding master equation reads

[TABLE]

Projecting this Fokker-Planck equation on the Fourier basis (32), one obtains

[TABLE]

For ABPs, the quasi-equilibrium limit is obtained by taking $D_{\rm r}\to\infty$ while keeping $D=v_{0}^{2}/2D_{\rm r}$ fixed. It turns out that up to $k=2$, the hierarchy of equations is essentially the same as for RTP, up to the replacement of $\alpha$ by $D_{\rm r}$. The only difference is in the equation for $k=2$, where a factor $k^{2}=4$ appears in the angular diffusion term. Taking into account the fast relaxation of the $f_{k}$ for $k\geq 1$, the closed equation over the density field $\rho(\bm{r})$ reads

[TABLE]

Eventually, looking for a stationary solution with a similar ansatz (40) as in the RTP case,

[TABLE]

one gets, assuming $U(\bm{r})=U(x)$,

[TABLE]

Conclusions are the same as for RTPs in 2D since the expressions are very closed to each other, except for numerical prefactors. In particular, one recovers that the stationary densities $\rho_{A}$ and $\rho_{B}$ depend on the detailed shape of the potential energy barrier.

As a check of the validity of the perturbative expansion, Fig. 1 compares the stationary density profile (averaged in the $y$ direction) measured in a numerical simulation of ABPs for a reasonably large value of $D_{r}$, with the perturbative analytical prediction. The potential barrier is uniformly smooth but asymmetric. This asymmetry generally leads to different densities in the bulks of each subsystems ($x>0$ and $x<0$) contrary to the equilibrium situation (see details in subsection III.4).

III.3 Active Brownian particles in an external field

As a last example, we consider the case of ABP with a stochastic reorientation in a given direction $\theta_{0}$, with some noise. More precisely, with a rate $\lambda$, a new direction $\theta^{\prime}$ is chosen randomly from a distribution $\psi(\theta^{\prime}-\theta_{0})$. For simplicity, the distribution $\psi(\theta)$ is assumed to be symmetric, $\psi(-\theta)=\psi(\theta)$. As mentioned in the introduction, this model is inspired by the dynamics of the microalgae Chlamydomonas Reinhardtii, which can reorient its direction of motion towards a light source.Martin et al. (2016) On a theoretical level, it can be thought of as a way to include an external field in a model of active particles. The evolution equation for the joint distribution $P(\bm{r},\theta)$ reads

[TABLE]

Expanding in angular Fourier modes, one has

[TABLE]

with $\psi_{k}=\int_{-\pi}^{\pi}d\theta\,\cos(k\theta)\psi(\theta)$. Similarly to the case of ABP, we perform an expansion at large angular diffusion coefficient $D_{r}$, to leading order in $1/D_{r}$. From now on, we assume that the potential energy is invariant along the direction $y$, and we choose $\theta_{0}=\frac{\pi}{2}$ so that the external field is parallel to the energy barrier. The stationary equation for the density profile $\rho(x)$ (invariant in the $y$ direction) reads

[TABLE]

with

[TABLE]

Assuming the form Eq. (45) for the density $\rho(x)$, we eventually obtain

[TABLE]

Again, conclusions are the same as the two previous 2D independent active particles models since expressions are very similar, except for the non-crucial term proportional to $\lambda$ in $\phi_{1}$ that results from the presence of the external field. As a result, we see that the dependence of the bulk densities $\rho_{A}$ and $\rho_{B}$ on the detailed shape of the potential energy barrier is a generic effect, that does not require any fine-tuning of parameters.

III.4 Chemical potentials for independent active particles in 2D

In the three active particles models in 2D, the stationary density profile $\rho(x)$ has been computed perturbatively and generally reads

[TABLE]

with $\phi(x)=\phi_{0}(x)+\tfrac{1}{\zeta}\phi_{1}(x)$, $\zeta=\alpha\text{ or }D_{\rm r}$ depending on the model considered.

Since particles are independent, the dynamics on the global density in each subsystems A and B is also expected to be Poissonian as for the RTP in 1D. Therefore, the large deviation function $\mathcal{I}(\rho_{A},\rho_{B})$ associated with $P(\rho_{A}|\bar{\rho})$ (see Eq. (17)) is exactly the same as eq. (18). The stationary densities in the bulks of each systems, $\rho_{A}^{\ast}$ and $\rho_{B}^{\ast}$, are obtained from the condition that the density profile is continuous in $x=0$ (a condition valid here because we have considered an equal speed $v_{0}$ in both compartments), leading to the condition

[TABLE]

with $x_{k}^{\ast}$, $k=A,\,B$, two arbitrary points in the uniform bulks of systems $A$ and $B$. Contrary to the equilibrium situation ($\zeta\to\infty$) for which $\phi(x)=\beta U(x)$, the presence of a non-local contribution proportional to $\int_{0}^{x}\mathrm{d}q\,[U^{\prime}(q)]^{3}$ in $\phi(x)$ for each model leads to an additional bias as long as $\zeta^{-1}>0$.

The associated chemical potentials are thus very similar to the RTP in 1D (24) and read

[TABLE]

with

[TABLE]

$K$ being a numerical constant which is equal to $\tfrac{1}{2}$ for RTP and $\tfrac{7}{8}$ for both ABP and ABP with an external field. Equalization of the chemical potentials $\mu_{A}^{\rm cont}$ and $\mu_{B}^{\rm cont}$ yields back the densities $\rho_{A}^{\ast}$ and $\rho_{B}^{\ast}$. Fig. 2 displays the ratio $\rho_{A}^{\ast}/\rho_{B}^{\ast}$ as a function of the asymmetry of the potential barrier. A strong effect is observed, even for relatively large values of $D_{\rm r}$.

We thus conclude that the nonequilibrium chemical potential defined in models of independent active particles does not satisfy an equation of state in general, as it does not depend only on the bulk density, but also on the shape of the energy barrier. As mentioned in Ref. Guioth and Bertin, 2018, the zeroth law of thermodynamics can be restored if one includes the (half) barrier in the definition of the system.

IV Relation between chemical potential and pressure

The lack of an equation of state for the chemical potential is reminiscent of a similar situation occuring for the pressure in models of active particles. As mentioned in the introduction, many studies on the notion of pressure in active systems have been conducted recently Solon et al. (2015b, a, 2018a); Fily et al. (2017); Solon et al. (2018b); Joyeux and Bertin (2016). The mechanical pressure can be computed as the average force exerted by the particles present in compartment $k=A$, $B$ on the wall separating the two compartments:

[TABLE]

(we recall that $x_{k}^{\ast}$ is a point in the bulk of compartment $k$). When studying the pressure, the potential energy barrier is assumed to be infinite. In constrast, in the study of the chemical potential, the barrier needs to be finite (though high to implement a weak contact) to allow for particle transfer between compartments.

For simplicity, we discuss here the connection between chemical potential and pressure in the framework of the one-dimensional RTP model, but generalization to other models is straightforward. If the system behaves as if it was at equilibrium (which can be achieved in the 1D RTP by taking the limit $\alpha\to\infty$, $D=v_{0}^{2}/\alpha$ kept fixed with $\beta=\eta/D$ the effective temperature), the density profile reads

[TABLE]

(where $\rho_{k}^{\ast}$ is the bulk pressure in compartment $k$), leading to a pressure on the side $k$ of the wall

[TABLE]

consistently with what one would expect for a perfect gas. But since $\rho(x)$ is non-local in $U(x)$ for RTPs in 1D, $P_{k\to{\rm wall}}$ has been shown to depend on the external potential barrier profile $U(x)$. It turns out that the pressure can be exactly computed here (see, e.g., appendix of Ref. Solon et al., 2015a in the case of an infinite barrier):

[TABLE]

At equilibrium, there is a well-known relationship between pressure and chemical potentials which is often called a Maxwell relation (see Refs. Sekimoto, 2010; Sasa and Tasaki, 2006 or Refs. Takatori and Brady, 2015; Takatori, Yan, and Brady, 2014 in the context of active matter). The Maxwell relation reads

[TABLE]

It is important to note that the chemical potential $\tilde{\mu}$ appearing in the Maxwell relation (60) has the dimension of an energy, and is related to the dimensionless chemical potential $\mu$ considered in this paper through $\mu=\beta\tilde{\mu}$, with $\beta$ the inverse temperature of the equilibrium system considered. Out of equilibrium, the temperature of the system is, strictly speaking, ill-defined. Yet, it is possible to circumvent the problem by defining an effective temperature $\beta_{\rm eff}^{-1}=v_{0}^{2}/(\alpha\eta)$, and a chemical potential $\tilde{\mu}=\beta_{\rm eff}^{-1}\mu$. In this case, one can wonder whether the Maxwell relation is still valid. For independent RTPs in one dimension, one obtains

[TABLE]

and

[TABLE]

Hence strictly speaking, the Maxwell relation is not valid. We note, however, that in the weak contact limit needed to define the nonequilibrium chemical potentialGuioth and Bertin (2018), the energy barrier separating compartments A and B is high, implying that $\Delta Q_{k}$ defined in Eq. (12) is large. It thus follows that the correction to the Maxwell equation is small in this regime. In addition, it might also be argued that the definition of the effective temperature $\beta_{\rm eff}^{-1}$ is phenomenological and not necessarily well grounded, so that including appropriate corrections to this effective temperature might restore the Maxwell relation (on condition, of course, that proper justifications to the modified effective temperature be given).

It turns out that more general arguments in favor of the breaking of the Maxwell relation can be given. The first argument goes as follows. As shown in Ref. Sasa and Tasaki, 2006, the Maxwell relation is intimately related to a balance between forces applied by an external potential on the particles and by the pressure exerted by the particles. Yet, since the two stationary densities $\rho_{k}^{\ast}$ for $k=A$, $B$ are related through the equality of currents (11), one finds that

[TABLE]

as long as $v_{A}\neq v_{B}$ and $\Delta Q_{A}\neq\Delta Q_{B}$. Hence, there exists a net force applied on the wallSolon et al. (2015a), which needs to be balanced to keep the potential energy barrier $U(x)$ centered at $x=0$. In reaction, the wall exerts a net force on the particles which needs in turn to be balanced for the average current to vanish. Such a balance is provided by the active force, understood as the force exerted by the environnement on the active particles.Fily et al. (2017) This means that at the end, one should experimentally observe a deformation or a flow in the substrate. As a consequence, if the force on the wall is not balanced, the wall is movingSolon et al. (2015a); Joyeux and Bertin (2016).

An alternative interpretation of the breaking of the Maxwell relation comes from the large deviation approach, which lies at the root of our definition of the chemical potential. At equilibrium, the pressure $P$ and the chemical potential $\tilde{\mu}$ are directly related to the free energy $F$ through the definitions $P=-\partial F/\partial V$ and $\tilde{\mu}=\partial F/\partial V$. The scaling relation $F(N,V)=V\,f(N/V)$ straightforwardly leads to the Maxwell relation (60). In terms of large deviation functions, this means that the large deviation function associated with an exchange of particles is $I_{N}(\rho)=f(\rho)$, and the one associated with an exchange of volume is $I_{V}(v)=vf(1/v)$, with the same function $f$ in both cases. Following Ref. Bertin, Dauchot, and Droz, 2006; Bertin et al., 2007; Guioth and Bertin, 2018, the definition of an intensive parameter conjugated to a conserved quantity can be extended out of equilibrium through the use of large deviation functions. However, as shown in Ref. Guioth and Bertin, 2018, the large deviation function generically depends on the contact dynamics, so that exchanging particles or volume no longer leads to the same scaling function $f$, at odds with equilibrium (the thermodynamic pressure could even be ill-definedGuioth and Bertin (2018)). This difference with the equilibrium case is thus another reason explaining the breaking of the Maxwell relation (60).

V Conclusion

In this paper, we have discussed the notion of nonequilibrium chemical potential in the simple framework of gases of noninteracting active particles placed in two compartments separated by a high energy barrier. Our definition of a nonequilibrium chemical potential is well-grounded in the theory of large deviation functionsGuioth and Bertin (2018), and does not rely on the assumed validity of the Maxwell relation connecting pressure and chemical potentialTakatori and Brady (2015). The main prediction of Ref. Guioth and Bertin, 2018 is that the nonequilibrium chemical potential generically depends on the contact dynamics. We have confirmed this result for several models of noninteracting active particles in one and two dimensions.

Although for two-dimensional systems our arguments are based on a perturbative expansion in the vicinity of an effective equilibrium limit, this result is enough to show qualitatively that the nonequilibrium chemical potential generically depends on the detailed shape of the potential energy barrier between the two compartments. In other words, the nonequilibrium chemical potential does not satisfy an equation of state in terms of bulk properties of the gas of active particles, but the shape of the potential barrier has to be taken into account. As already mentioned above, this result is reminiscent of a similar result obtained when evaluating the mechanical pressure exerted by active particles on a (soft) wall modeled by a diverging potential energy profile. An interesting extension of this work could be to try to incorporate interactions between particles in the evaluation of the chemical potential. Although no conclusion can be drawn at this stage, it is unlikely that the inclusion of interactions could restore a bulk equation of state. At least this is not the case when considering pressure Solon et al. (2015a).