A Simulation details

# Effective equilibrium states in the colored-noise model for active matter I. Pairwise forces in the Fox and unified colored noise approximations

## Abstract

The equations of motion of active systems can be modeled in terms of Ornstein-Uhlenbeck processes (OUPs) with appropriate correlators. For further theoretical studies, these should be approximated to yield a Markovian picture for the dynamics and a simplified steady-state condition. We perform a comparative study of the Unified Colored Noise Approximation (UCNA) and the approximation scheme by Fox recently employed within this context. We review the approximations necessary to define effective interaction potentials in the low-density limit and study the conditions for which these represent the behavior observed in two-body simulations for the OUPs model and Active Brownian particles. The demonstrated limitations of the theory for potentials with a negative slope or curvature can be qualitatively corrected by a new empirical modification. In general, we find that in the presence of translational white noise the Fox approach is more accurate. Finally, we examine an alternative way to define a force-balance condition in the limit of small activity.

## 1 Introduction

Active Brownian particles (ABPs) provide a simple, minimal model system to study the collective behavior of active matter. Many-body Brownian dynamics simulations of these systems have provided considerable insight into a range of interesting nonequilibrium phenomena, such as the accumulation of particles at boundaries (1); (2); (3); (4); (5) and motility-induced phase separation (6). Much of the phenomenology of ABPs can be captured using coarse-grained, hydrodynamic theories (6); (7); (8); (9); (10), which do not contain all information about the interparticle correlations. Some progress has recently been made in the linear response regime, which allows to decouple the equations of motions of the one-body density and polarization vector (11).

Due to the inherent difficulty of dealing simultaneously with both the translational and orientational degrees of freedom in active systems, attempts to develop a first-principles theory have largely focused on a simpler, related model, in which the particle dynamics are represented by a set of coupled Ornstein-Uhlenbeck processes (OUPs). Within this model an exponentially correlated noise term, with a given correlation time, serves as proxy for the persistent trajectories of ABPs (connections between the two models were explored in Ref. (12)). While the removal of orientational degrees of freedom does indeed simplify the problem, it comes at the cost that one has to deal with the non-Markovian dynamics of the translational coordinates. Fortunately, there exist various approximation methods (17); (18); (15); (16); (13); (14); (19) which enable the OUP model to be represented using an effective Markovian, and therefore tractable, dynamics. Two different approaches to doing this, (i) the Unified Colored Noise Approximation (UCNA) of Hänggi et al. (13); (14), based on adiabatic elimination on the level of the Langevin equations, and (ii) the Fox approximation (15); (16), for which an approximate Fokker-Planck equation is developed, have recently been adopted in the context of developing simple theoretical tools to describe active particles (12); (25); (20); (21); (23); (24); (22); (26).

When applied to active matter, both the UCNA and Fox approximations are referred to as ‘effective equilibrium’ approaches. The Markovian character of the dynamics implies that they obey a Fokker-Planck equation from which an effective probability distribution can in principle be obtained. Indeed, the possibility of mimicking the behavior of nonequilibrium ABPs using an equilibrium system of passive particles, interacting via effective interactions, was suggested by several researchers (see, e.g., Ref. (27)) who observed that the phase separation induced by activity in systems of repulsive ABPs closely resembles that familiar from passive systems with an attractive interaction. Despite its appeal, several years were required before this observation could be turned into something more concrete. By starting from the simpler OUP model it became possible, via application of the UCNA (20); (21); (23); (24); (22) and Fox (12); (25); (26) approximations, to put the notion of an effective equilibrium description on a firmer footing.

In this paper, we will compare and contrast the two different approaches to effective equilibrium. We will highlight the main approximations involved and assess the validity of the effective-potential approximation (EPA), which has been employed in previous work to investigate activity-induced modifications of the microstructure and motility-induced phase separation (12). This analysis both clarifies the nature of the approximations involved and suggests ways in which the description can be improved.

The paper is laid out as follows: In Sec. 2 we first specify the model under consideration and describe the UCNA and Fox approaches to obtaining an effective equilibrium picture highlighting similarities and differences between them. In Sec. 3 we describe in detail the EPA, where the emphasis is placed on the UCNA due to its simpler structure. The resulting approximate effective potentials are compared to computer simulations using a standard soft-repulsive and a non-convex (Gaussian core) potential. In Sec. 4 we consider an alternative approach to obtain pairwise forces, i.e., the low-activity limit, and make contact to the EPA. Finally, we conclude in Sec. 5.

## 2 Theory

In this section, we introduce the common starting point of both the UCNA and Fox approach. Since particles driven by Gaussian colored noise originally were not intended as a model for an active system, the choice of parameters in the literature may depend on the dimensionality and on whether contact to ABPs is made (12) or not (20). We will also clarify some notational issues.

### 2.1 Colored-noise model

We consider the coupled stochastic (Langevin) differential equations

 ˙ri(t)=γ−1Fi(r1,…,rN)+ξi(t)+vi(t) (1)

of particles. The motion of each particle at position is determined by conservative and stochastic forces and . The friction coefficient is related to the translational Brownian diffusivity and is the inverse temperature. We assume that the total interaction force can be written as the gradient of a pairwise additive many-body potential

 U(rN)=(\varv(ri)+12N∑k≠iu(ri,rk)), (2)

consisting of the one-body external fields and the interparticle potentials .

The vector represents the translational Brownian diffusion by a Gaussian (white) noise of zero mean and with the unit matrix . Here and in the following the dyadic product of two vectors with components results in a matrix. Any contraction as in a scalar product or a matrix-vector product will be explicitly indicated by a “”. Hereafter, we shortly refer to the variable as (thermal) noise. The OUPs defined by

 ˙vi(t)=−vi(t)τa+ηi(t)τa (3)

with describe a fluctuating propulsion velocity as a non-Gaussian (colored) noise of zero mean and

 ⟨vi(t)vj(t′)⟩=v20d1δije−|t−t′|τa=Daτa1δije−|t−t′|τa. (4)

Here we introduced the active time scale at which the orientation randomizes and the active diffusion coefficient , where is the average squared self-propulsion velocity and the spatial dimension.

The colored-noise model for active particles contains two parameters describing the magnitude and persistence of the self propulsion. We now aim to clarify some notational differences in the literature. The persistence time can be explicitly related to the equations of motion of run and tumble particles (28) or ABPs (29); (12). The above definitions correspond to the latter case with , where is the rotational Brownian diffusion coefficient. Another common choice (23) amounts to consider and . In the following, we use , where is the typical diameter of a particle and the dimensionless persistence time has been introduced in Refs. (12); (25).

Since the dimensionless active diffusivity implicitly depends on the persistence time, it constitutes the most general measure for the activity (together with the persistence length ). In order to connect to a system of ABPs (12); (25), it is convenient to consider instead of a dimensionless velocity , i.e., the Peclét number. In the literature, some other definitions of a Peclét number are used, which we will not consider here. One peculiar property of a system of active OUPs is that, even at vanishing self-propulsion velocity , or , there is a contribution of to Eq. (1) arising from a finite reorientation time  (26). One thus does not recover the equation of motion of a passive (Brownian) particle, as in the ABPs model. In the long-time limit, however, the contribution to the dynamics becomes irrelevant and the same steady state is described as for a passive Brownian particle, see appendix A for more details. A proper passive system can only be recovered from Eq. (1) only in the limit , in the sense that the velocity correlation in Eq. (4) reduces to a white noise. A Brownian system is then represented by when the thermal-noise variable is removed, or, trivially, by setting which amounts to neglect the contribution of the OUPs.

### 2.2 Effective equilibrium approach

The most important step towards a theoretical study of the OUPs model is to derive from the non-Markovian stochastic process (1) an equation of motion for the -particle probability distribution . In this section, we will discuss the differences between the multidimensional generalizations of the UCNA (13); (14) and the Fox (15); (16) approaches to effective equilibrium and expound the surprising similarities between these two approximations in the (current-free) steady state.

As a central quantity emerging in both cases, we define the friction tensor with the components

 Γij(rN)=1δij−~τ∇iFj=δijΓii(rN)+(1−δij)~τ∇i∇ju(ri,rj) (5)

resulting in the Hessian of and the diagonal block

 Γii(rN):=1+~τ∇i∇i(\varv(ri)+N∑k≠iu(ri,rk)) (6)

not to be confused with for particle. In the following, we briefly denote by the th block component of the inverse tensor .

The UCNA (20); (21) amounts to explicitly inserting the OUPs (3) into the overdamped limit of the time derivative of (1), resulting in the modified Langevin equation . It is now straight-forward to obtain for this (approximate) Markovian system driven by white noise the Smoluchowski equation with the probability current (the superscript (u) denotes that the UCNA has been used)

 Missing or unrecognized delimiter for \bigg (7)

Note that the UCNA remains valid as long as the friction tensor (5) is positive definite.

The Fox approximation scheme applied to (1), on the other hand, only makes use of the correlator (4) of the OUPs, which, in turn, may also be interpreted as the correlator of corresponding to a coarse-grained equation of motion representing ABPs with a constant velocity in the direction of their instantaneous orientation that is subject to Brownian rotational diffusion (12). This method directly yields the approximate Smoluchowski equation (superscript (f)) with (26); (30)

 J(f)i=Dt(βFifN−∇ifN−Da∑j∇j⋅(Γ−1jifN)), (8)

where the regime of validity is the same as for UCNA. The major difference between Eq. (7) and (8) only affects the effective description of the dynamics as a result of the additional factor arising on the level of the Langevin equation within the UCNA. Note that in the original generalization of the Fox result (12) the tensor from Eq. (5) was incorrectly obtained as , which we will later identify as the (diagonal) Laplacian approximation. It will turn out that this (or another) approximation is necessary to obtain physical expressions for the effective interaction potentials.

### 2.3 Two versions of the steady-state condition

In contrast to the dynamical problem, the (current-free) steady-state conditions

 βFiPN−∑j∇j⋅(DjiPN)=0 (9)

for the stationary distribution can be cast in a coherent form, defining the effective diffusion tensor , such that only the components

 D(u)ij(rN) :=(1+Da)Γ−1ij(rN), (10) D(f)ij(rN) :=1δij+DaΓ−1ij(rN). (11)

differ between the UCNA (u) and Fox (f) results.

Multiplying Eq. (9) with and summing over repeated indices, the steady-state condition takes the more instructive (approximate) form (21)

 0=∑iD−1ik⋅βFiPN−∇kPN−PN∇kln|detD[N]|=:βFeffkPN−∇kPN (12)

introducing the effective force . The term is an approximation for , which becomes exact in the UCNA (21). For the Fox approach, we argue in appendix B that this is still true in some important special cases, such that Eq. (12) is accurate enough for our purpose. For high particle numbers the contribution of the off-diagonal elements to becomes increasingly irrelevant (21), which amounts to setting . Assuming this diagonal form, the determinant in Eq. (12) can be replaced according to as we have before approximating the expression in the last step (compare appendix B).

Putting aside the dynamical behavior, described in Sec. 2.2, the only difference between the UCNA or Fox approximation is manifest in the definitions, (10) and (11), of . Using UCNA the active diffusivity only appears as part of a prefactor in (10), so that the friction matrix , representing a correction due to activity, contributes to the steady-state result even in the case , that is when and . Hence, the logical parameter suggested by the UCNA to tune the activity is , with the passive system () restored only in the limit . This appears to be an artifact of the pathological contribution of to the displacement of the OUPs, whereas the connection to the experimentally more relevant system of ABPs is lost. For the latter it appears more natural to tune at constant . In the derivation of the Fox result (11), on the other hand, the explicit time evolution of the OUPs in Eq. (3) is irrelevant, suggesting a better approximate representation of ABPs (12). This reflects that we recover the (same) passive system for either or (in the presence of noise).

Ignoring the noise contribution for , the UCNA and Fox approximations practically describe the equivalent effective steady states. The major advantage of this approximation, or the UCNA result in general, is that the inverse is pairwise additive, even if . Then the effective many-body potential defined as can be written in a closed form (20); (21), admitting the explicit solution of Eq. (9). Due to the more nested form of Eq. (11) the Fox approximation does in general not admit an analytic result. As is not pairwise additive in either approach, some further approximations will become necessary to construct a predictive theory, which we discuss in the following sections.

## 3 Effective-potential approximation (EPA)

Regarding the possible applications using standard methods of equilibrium liquid-state theory a desirable strategy is to approximate in Eq. (12) in terms of pair potentials. This approach allows to describe the phase behavior of ABPs approximated as particles propelled by a set of coupled OUPs, which has been discussed in detail (12); (25) for passive soft-repulsive and Lennard-Jones interactions in three dimensions. However, it can be criticized that (I.i) a system which obeys detailed balance is used to represent the interactions in an active system, (I.ii) the validity criteria of the underlying theory might be violated so that further approximations are required and (I.iii) higher-order particle interactions are neglected, which are believed to be important for the phase separation in an active system. In the following, we define the effective pair interaction and motivate different approximations, which we compare to computer simulations of two ABPs and two particles propelled by OUPs. It is our objective to comment on the aforementioned points and illustrate the qualitative differences between the Fox and UCNA. For the sake of simplicity, we will restrict the presentation of technical aspects to the UCNA results.

### 3.1 Calculation of effective potentials

To identify an effective pair potential , we consider interacting particles, i.e., the low-density limit of Eq. (12). Ignoring the external forces for now by setting , it is easy to verify that

 ∇1βueff(r)=(D−111−D−121)⋅∇1βu(r)+∇1ln|detD[2]| (13)

and analog equation for , where we used and . Keeping in mind that we seek to employ this effective potential to approximately represent the interaction of many particles, it appears undesirable that an equal statistical weight is put to both the diagonal and the off-diagonal components of the diffusion tensor. As an alternative we propose the effective potential

 ∇kβueffdiag(r)=D−1kk⋅∇kβu(r)+∇kln|detDkk|, (14)

with , obtained for a diagonal form of with . For completeness we find in the one-particle limit a quite similar formula

 Extra open brace or missing close brace (15)

for the effective external field , since for we have . Note that a quite different expression for an effective external potential can be derived starting from the equations of motion for ABPs (25); (31).

Integration of the above equalities yields the desired formulas for the effective potentials depending only on the bare potential or and the activity parameters and  (12). Alternatively, we could have directly defined (21); (23) and from the many-body potential identified in the solution of (9), which is, however, inconvenient when the Fox approach is used. Assuming a bare potential obeying , the integrated form of (13) reads

 βueff(r) =βu(r)+~τ(∂ru(r))21+Da−ln∣∣E(d−1)1(τ,r)E2(τ,r)∣∣, (16)

where and

 En(τ,r):=1+2~τrn−2∂nru(r), n∈{1,2} (17)

are the Eigenvalues of . We can further identify in the first term of Eq. (13) as the Eigenvector of corresponding to the Eigenvalue . Note that in (16) we could equally introduce an effective energy scale to absorb the factor  (21); (23). We refrain to do so as this interpretation would not be consistent with the the way enters within the Fox approach.

### 3.2 Limitations and possible corrections

Studying Eq. (16) more carefully, we notice that the effective potentials do not always behave in a physical way. This is because, in violation of the validity condition of both the UCNA and the Fox approximation, the diffusion tensor is not positive definite for a large number of relevant potentials. In general, we easily see that the logarithm will diverge whenever one of the Eigenvalues vanishes. Given a positive and convex bare potential , the Eigenvalue is strictly positive, which also means that the effective attraction solely arises from the term including the logaritm. However, as we have in this case, the Eigenvalue will vanish at a certain value of and we require a further approximation to remedy the unphysical behavior of in dimensions. At a highly non-convex or negative region of the bare potential, the same problem occurs for . Interestingly, if we only require knowledge of an effective potential on a finite interval where the eigenvalues are positive, its overall unphysical behavior is irrelevant (26).

First note that there is a broader range of admissible bare potentials when the diagonal approximation, Eq. (14), of the effective pair potential is used, or if we are interested in the one-body external field, Eq. (15). This can be understood from the explicit formula for , which we obtain from Eq. (16) by rescaling all terms proportional to with a factor . In the following, we propose different ways to generally rid the effective potential of possible artifacts of vanishing Eigenvalues in the last term of Eq. (13). A correction of the first term is not necessary and also has no noticeable effect.

Let us first assume that is convex, i.e., it represents a soft-repulsive interaction. Then a sufficient criterion (due to the presence of the term in Eq. (5), some other potentials are allowed that are only slightly negative and slightly non-convex) for the matrix to have strictly positive Eigenvalues would be that it depends on an elliptic differential operator rather than . Therefore, a convenient approximation is to redefine Eq. (5) by an elliptical operator, the simplest example of which is the Laplacian . Upon substituting

 ∇i∇j→1∇i⋅∇j (18)

the effective potential becomes

 βueffΔ(r)=βu(r)+~τ(∂ru(r))2−2(d−1)∫∞rds~τ(∂su(s))2s1+Da−ln(1+2~τ∂2ru(r)+2(d−1)~τ∂ru(r)r) (19)

where the additional term compared to (16) cannot be integrated in general. This Laplacian approximation has been successfully employed (together with the Fox and diagonal approximation) in explicit calculations (12); (25). In dimensions both differential operators reduce to the second derivative and Eq. (19) is equal to Eq. (16), which provides a good account of active particles interacting with a soft-repulsive potential (23).

An alternative way is to empirically rectify the explicit formula for in Eq. (16). Most intuitively, one can expand the argument of the logarithm up to the first order in . In fact, this small- approximation is quite similar to the Laplacian approximation (19) (and completely equivalent in one dimension), but we do not recover the additional term involving the integral. Performing the small- approximation of the full expression (16) appears too crude, as an expansion of the logarithm does not converge for . The resulting effective potential will thus become totally uncontrolled for short separations of highly-repulsive particles.

A more elaborate correction that may be applied also to highly non-convex potentials is the inverse- approximation, an empirical strategy maintaining the leading order in while not disregarding higher-order terms. This is achieved by substituting in Eq. (16) , where

 E(i)n(τ,r):={1/(2−En(τ,r))En(τ,r)if En(τ,r)<1otherwise. (20)

The major advantage of this approximation is that it yields quite similar results to the full potential whenever the validity condition is only slightly violated and the effective potential does not diverge if the bare potential is finite. The empirical motivation behind this correction is that constitutes the two leading terms of the “resummed” Taylor series of in the case . As described in appendix B, the most convenient implementation of the inverse- approximation for the Fox result is to identify the expressions for in and use Eq. (20).

### 3.3 Comparison to computer simulations of two active particles

In Sec. 3.1 we introduced different strategies to define a suitable effective interaction potential in the effective-equilibrium approximation for the colored-noise model. Now we illustrate under which conditions an approximate treatment according to Sec. 3.2 becomes necessary and compare the theoretical results to computer simulations. The easiest way to determine an effective potential numerically is to set up a two-particle simulation, measure the radial distribution function and calculate . By doing so, we make the same approximation (I.iii) as in theory to ignore the many-particle character of the interaction. However, the simulations for ABPs and OUPs, detailed in appendix A, take into account the orientation dependence and the non-Markovian character of the dynamics, respectively.

#### The role of approximations, dimensionality and thermal noise

We first discuss some general observations in the UCNA for a soft-repulsive system with the bare potential . The behavior of the Fox result is qualitatively similar. As expected, the full expression for the effective potential in Eq. (16) is impractical as it diverges at a certain distance , determined by the condition , which is when the first Eigenvalue within the logarithm vanishes, whereas is always positive. As suggested by Fig. 1a, this behavior is most problematic at larger values of , where should be rather negative, as it is the case in dimensions. We further see in Fig. 1a that this effect becomes more severe with increasing dimension. Both the inverse- and Laplacian approximations successfully cure this unphysical divergence, which we see in Fig. 1b. As employing the diagonal form of the effective potential simply amounts to a rescaling of , we observe in Fig. 1b that it results in a smaller effective diameter of the repulsive part but a flatter potential well. Accordingly, becomes smaller.

Since the definition of the active diffusivity depends on the dimension , not all parameters , and can be kept constant upon varying the dimensionality. For a constant reorientation time and propulsion velocity the effective attraction in Fig. 1c is stronger in lower spatial dimensions for both approximations considered, which coincide with the full expression in . This behavior appears sensible, as two particles have less possibilities to avoid each other upon collision, and is in qualitative agreement with computer simulations of active OUPs. Moreover, we understand that motility-induced phase separation is harder to observe in higher dimensions (32). Keeping constant instead of (which then decreases with decreasing dimension) the same trend is observed in Fig. 1d for the inverse- approximation and computer simulations, whereas the result for the Laplacian approximation barely changes with dimensionality. Comparing the qualitative behavior in Figs. 1c and d, we recognize in all spatial dimensions that the numerical effective potential is of longer range than the theoretical predictions in any approximation. This observation confirms the criterion discussed in Ref. (14) that the UCNA is expected to become less accurate for larger separations where a typical length scale of the active motion, closely related to the effective diffusion tensor, Eq. (10), exceeds the spatial scale over which the force field varies.

As the involved approximations become cruder in higher spatial dimensions, the quantitative agreement with the simulation results in Figs. 1c and d becomes worse. For , a remarkable agreement between the UCNA and simulation results for the radial distribution of two particles has been reported in Ref. (20), where in Eq. (1) was set to zero. Doing so also in our simulations, we observe in Fig. 2a that the effective potential deepens and its repulsive barrier becomes steeper. This curve is in excellent agreement with the theoretical result for zero noise, obtained by both the UCNA and Fox approach upon dropping the first term in Eq. (10) and Eq. (11), respectively. The full UCNA result is only slightly different in the repulsive regime. Intriguingly, we also recognize in Fig. 2a that the simulation data including the noise term are excellently represented by the Fox approach. We can understand these observations by recapitulating the idea behind the two approximations. The UCNA amounts to manipulate Eq. (1) by calculating the second derivative of in order to eliminate the variable in Eq. (3). The original discussion of the accuracy of this approximation does not account for the second stochastic variable . In contrast, the Fox approach is only dedicated to determine the approximate contribution of the colored noise to the effective probability current in Eq. (8), which is independent of other terms in Eq. (1). Therefore, the Fox theory has a broader range of applicability and should be accurate in both the presence and the absence of thermal noise.

Finally, we note that the excellent agreement between theory and computer simulations in one dimension implies that the diagonal approximation is not justified when it comes to describing a two-body system. It is, however, interesting to consider a single particle in an external field of the same form as the interparticle potential considered above. In agreement with the theoretical prediction, the computer simulations in Fig. 2b show nice agreement between the two-body system and a one-body system with the double value of the persistence time. We further observe that the theoretical result for one body is even closer to the simulation data than for two bodies.

#### Soft-repulsive Brownian system in three dimensions

We also performed computer simulations of ABPs (described in appendix A) for . For a finite active diffusivity, Fig. 2c reveals that the numerical effective potentials for the two considered models with and without noise are nearly identical over the full range of separations. As for , the effective potential of active OUPs (and ABPs) in the absence of thermal noise has a deeper well and a larger repulsive diameter. Quantitatively, this difference is much more pronounced in three dimensions. In the following, we restrict ourselves to systems with thermal noise and compare in Figs. 2d and 3 to the predictions of the theory. As in Figs. 2c and d, our simulations of ABPs and OUPs, shown in the first column of Fig. 3, are in nice agreement for all sets of considered parameters. This is quite surprising since, on the many-particle level, ABPs and OUPs have different steady states (33); (34). On the basis of our data for the simplistic two-body system we could rather conclude that OUPs subjected to thermal noise are an excellent model for ABPs at moderate activity (12); (25).

We see in Fig. 2d that the depth of the attractive well of all theoretical versions of the effective potential in dimensions is significantly overestimated when compared to the simulations of both ABPs and OUPs. The inverse- approximation appears to provide the best guess of the point at which the effective potential changes its sign. To facilitate the further qualitative comparison we chose the axes in Fig. 3 according to the deviation from the simulations (first column), i.e., by a factor of 10 for the Laplacian approximation (second column) and 5 for the inverse- approximation (third column). The chosen approximations exhibit a behavior similar to the full theoretical results of Eq. (16) in the physical region for , shown in the last column. The differences arising from using the exact effective force in Eq. (12) are discussed in appendix B.

For the considered soft-repulsive bare potential, we observe in Fig. 3 some notable quantitative differences between the UCNA and Fox results, even at relatively high . The effective diameter of the repulsive part is generally smaller than in the UCNA, whereas the overall attraction is weaker in the Fox approach. This becomes most apparent in the Laplacian approximation. The first row of Fig. 3 contains the effective potentials evolving for a constant persistence time when the active diffusivity (or the Peclét number ) is increased. All approaches accordingly predict an increased effective attraction and the minimum of the potential is shifted to smaller separations (12); (25). The Fox results exhibit a stronger variation with , which is also more consistent with the numerical data. Similarly, the effective potentials in the second row deepen with increasing at constant , where the location of the minimum is nearly unaffected. The most interesting behavior is observed in the third column at constant . Again, all approaches agree that the minimum is shifted to larger separations with increasing , but the effective attraction predicted by the simulations is nearly constant, as simultaneously the magnitude of the self-propulsion is decreased. This observation is not consistent with the UCNA results.

Based on the presented simple comparison, our conclusion is that the best choice for the theoretical effective potential is the Fox approach in the inverse- approximation. Upon further increasing the activity (not shown) the quantitative discrepancy of the Laplacian approximation becomes even more pronounced. Regarding Fig. 1b one might get the impression that the additional assumption of the diagonal form of the effective potential also results in a slightly better (quantitative) agreement with the simulations. However, we stress that it is not clear in how far the effective pair potentials can accurately describe the many-body situation, as there are no higher-order interactions present in a two-body simulation.

#### Difficulties for non-convex potentials

The most compelling argument in favor of the inverse- approximation arises from considering non-convex bare potentials, a case in which the Laplacian approximation becomes useless above a certain value of . This issue has already been discussed for a Lennard-Jones potential (30) and we will consider here the potential of an active Gaussian-core fluid. Although this model has not received much attention in theories for ABPs, it is quite appealing from a theoretical perspective. Most prominently, this bare potential is a known exceptional case in which a simple mean-field theory is particularly accurate (35); (36), which might, in a way, also hold for the effective potential of the active system.

The effective potential of an active Gaussian-core fluid is discussed in Fig. 4 within the UCNA. As the absolute value of both the curvature and slope of this model potential is bounded, the Fox results (not shown) are quite similar, even for the moderate values of considered here. When the persistence time is sufficiently small, the effective potential, Eq. (16), does not diverge and the different approximations behave in a quite similar way. Interestingly, we observe in Fig. 4a that the divergence of the Laplacian approximation sets in at an even smaller value of than for the full potential. The latter diverges at two points, each related to one of the two Eigenvalues, if .

In Fig. 4b, c and d we discuss the only suitable form of the effective potential, i.e., the inverse- approximation. Intriguingly, the predicted behavior depends on in which way, i.e., by means of which parameter, the activity is modified. Increasing the active diffusivity (or the Peclét number) at a constant value of results in a less repulsive core. Upon increasing , however, the height of the maximum increases and an attractive well develops at larger separations. Quite counterintuitively, we observe that in this case the effective interaction becomes more repulsive than in the passive case, also for the inverse- approximation. Our computer simulations (also carried out for values of much larger than shown in Fig. 4) confirm that this is an artifact of the theory, related to the negative curvature of the bare potential. At constant the evolution of the theoretical results agrees qualitatively with the simulations. The simulation data are, however, not very sensitive to changes in the persistence time. At constant the theory predicts the correct trend upon increasing , whereas this is not the case at constant Peclét number.

To argue about the validity of the EPA, we consider two classes of bare interactions. Firstly, soft-repulsive and convex potentials lead to quite accurate results in one dimension, but require an empirical correction in higher dimensions. Secondly, the understanding of the behavior of particles interacting with a bare potential which has a negative curvature remains one of the most urgent open problems in our theoretical framework. At the moment, the only way to obtain a workable theory in this case is by employing the inverse- correction introduced in Eq. (20). Further numerical and theoretical analysis will be needed to fully clarify this issue. Finally, we note that potentials with attractive parts do not a priori constitute a problem for the theory, but usually have regions in which they are non convex. For a discussion of such problems see also Refs. (30); (26).

## 4 Low-activity approximation

A second strategy to simplify the steady-state condition is to perform an expansion in the activity parameter . At linear order, the effective diffusion tensor becomes pairwise additive. In this low-activity approximation, a YBG-like hierarchy can be obtained by successively integrating Eq. (9) over coordinates (21), which allows defining a mechanical pressure and interfacial tension (22). Moreover, for the active system evolving according to Eqs. (1) and (4) it has been demonstrated that there exists a regime for small values of , where the principle of detailed balance is respected (37). This suggests that, at leading order in this parameter, the approximations resulting in Eq. (9) are perfectly justified.

Knowing, however, that Eq. (9) contains the same information as Eq. (12), which depends logarithmically on the parameter , we should clarify whether (II.i) the low-activity expansion converges, (II.ii) it is sufficient to only consider the leading order and (II.iii) one can obtain similar results when employing the EPA. To do so, we demonstrate how the first member () of the YBG hierarchy can be rederived from Eq. (12) and discuss the consequences of approximating in terms of pair interactions. Again, we only discuss the UCNA, where, without any further approximation, the inverse diffusion tensor is found to be pairwise additive.

### 4.1 Alternative derivation of a local force-balance condition

By integrating a multidimensional vector equation (label ) over coordinates we understand multiplying each side by , followed by summation over all particles and integation over all spatial coordinates. We further define the average of a vector , where denotes the full canonical ensemble average and is the average of the density operator . Approximating now the inverse mobility matrix in Eq. (10) as and integrating Eq. (9) over coordinates, we find the first member

 −ρ(r)⟨∇βU⟩ =(1+Da)(∇⋅(1ρ(r)−~τρ(r)⟨∇∇U⟩)), (21)

of a YBG-like hierarchy (21); (22) for the active system, where, explicitly

 ⟨DU⟩=D\varv(r)+∫dr′ρ(2)(r,r′)ρ(r)Du(r,r′) (22)

for any (nontrivial) differential operator acting on . In the derivation of (21) it turns out that the off-diagonal components of the mobility tensor do not contribute at first order in  (22). Hence, we might as well have assumed the diagonal form at linear order in beforehand.

In order to connect to the EPA, we derive a YBG-like hierarchy from Eq. (12). Assuming the diagonal form , the integration over coordinates of the first equality is carried out i n appendix C. Making use of the equilibrium version of the YBG hierarchy and expanding the expression up to first order in the result is

 0=−D−1I(r)ρ(r)⟨∇βU⟩−∇ρ(r)+~τρ(r)⟨∇⋅∇∇U⟩+ρ(r)1+Da~τ∫dr′(∇∇u(r,r′))⋅∇ρ(2)(r,r′)ρ(r) (23)

introducing the averaged inverse diffusion tensor (compare Eq. (10))

 D−1I(r):=⟨Γii⟩1+Da=1+~τ⟨∇∇U⟩1+Da. (24)

Multiplying Eq. (23) with it is easy to verify in appendix C that at first order in it becomes equivalent to (21) up to a term proportional to the expression in the second line, which we consider as a higher-order contribution.

In order to derive Eq. (23) in the Fox approach, an additional approximation is required, as the inverse of from Eq. (11) is not proportional to . We would thus need to redefine in Eq. (24) according to Eq. (11) where takes the role of . Regarding in general the presented alternative derivation of Eq. (21), its validity appears to be in question. This is because to derive the intermediate result in Eq. (23) it is necessary to expand a logarithmic term, the Taylor series of which only has a finite radius of convergence. We further explicitly assumed the diagonal form of the mobility tensor to avoid further terms that are not present in the original result. Employing in the next step the EPA will shed more light on these issues.

### 4.2 Low-activity approximation of effective potentials

Having established a connection between (9) and (12) also at linear order in , we now turn to the case in which the second equality in (12) does not hold. This is when we assume along the lines of (2) but within the EPA using the results derived in Sec. 3.1. As detailed in appendix D, the obvious result is that all correlation functions between more than two particles vanish in the approximate integrated version

 Missing or unrecognized delimiter for \left (25)

of (12). Ignoring the interparticle interactions the approximation involving only becomes exact. This situation is the same as discussed in Ref. (21).

Considering the interacting system, we multiply Eq. (25) with as done previously for (23). According to appendix D, we can only approximately reproduce Eq. (21) by doing so. This reflects both the limitations of the EPA and an inconsistency between Eq. (9) and Eq. (12) when being subject to the same type of approximation, as we discuss in the following. We observe that (III.i) the coupling between external and internal interactions is ignored by Eq. (25) (III.ii) spurious three-body correlations appear on the left-hand-side of Eq. (21) (III.iii) the second term on the right-hand side of Eq. (21) is recovered but involves a seemingly unjustified expansion and (III.iv) if we do not explicitly assume a diagonal diffusion tensor, the last term in Eq. (21) changes by a factor two. As we are mainly interested in bulk systems, the first point is only briefly commented on in appendix D.

The term including the bare interaction force in Eq. (12) depends on the position of three bodies. Hence, the pairwise approximation, which amounts to setting

 Missing or unrecognized delimiter for \left (26)

should not be too crude. Moreover, we have discussed in Sec. 3.2 that such a contribution to the effective force is usually purely repulsive and thus plays only a minor role in characterizing a possible phase transition. However, we show in appendix D that the definition (24) of the averaged diffusion tensor is not fully compatible with the EPA, resulting in point (III.ii). This is in contrast to the clean derivation of Eq. (21), where Eq. (9) is recovered from Eq. (12) by multiplication with the many-body effective diffusion tensor before integrating over positions.

The last term in (12), although considered here for a diagonal diffusion tensor, constitutes a full -body quantity. Recall from the discussion in Sec. 3 that an approximation as a pairwise quantity might be quite poor and an expansion of the logarithm does not converge. However, we demonstrate in appendix D that successively employing the EPA and expanding for small according to

 ∇ln∣∣∣det (1+~τ∇∇∑l>1u(r,rl))∣∣∣⟶∇⋅∑l>1ln|det(1+~τ∇∇u(r,rl))|⟶~τ∑l>1∇Δu(r,rl)+O(τ2) (27)

eventually results in full consistency with the respective term in Eq. (21), stated as point (III.iii). This suggests that the expansion to first order in ”implies” making the EPA when Eq. (12) is our starting point. Despite the aforementioned crudity of this expansion, we argue that Eq. (21) is valid, as its clean derivation from Eq. (9) does not require dealing with a logarithmic term. The last step in Eq. (27) is required to recover Eq. (21) without inducing undesired higher-order terms in , as, similar to point (III.ii), the integrated version is incompatible with the chosen . However, we stress that this approximation should certainly be avoided when calculating the fluid structure.

Finally, we demonstrate in appendix D that the off-diagonal elements of the diffusion tensor entering in Eq. (12) contribute to Eq. (25). Hence, the present approach would be even more inconsistent with Eq. (21) if we did not assume the diagonal form, as noted in point (III.iv). We also note that the same problem occurs for the according generalization of Eq. (23). In principle we could define in this case an additional averaged diffusion tensor , corresponding to the off-diagonal elements, which could counteract this inconsistency. Such a calculation would, however, not be useful when an effective pair potential is employed.

## 5 Conclusions

In this paper we studied different ways to define an effective pair interaction potential between active particles. Our numerical investigation reveals that a two-particle system of ABPs and active OUPs exhibits a quite similar behavior. These results serve as a benchmark to test the approximations involved in recent effective equilibrium approaches, which have been reviewed and compared in detail. For spatial dimensions higher than one we introduced an empirical way to rid the theoretical result of possible divergences, which also appears to yield the best agreement with the simulation data, although the effective attraction is still significantly overestimated. Regarding the quite accurate one-dimensional results and the qualitative features of the effective potentials in three dimensions, the Fox approximation is superior to the UCNA when the translational Brownian noise cannot be neglected. In the absence of noise both approximation schemes admit the same steady-state solution.

Further analysis is needed to better understand the role of the neglected many-body interactions in both the two-body simulations and the theory, which are thought to be imperative for a quantitative description of active systems (30). The presented theoretical approach follows two major approximate steps to define the effective pair potential. First, we map the equation of motion (1) onto a deterministic Fokker-Planck equation (effective equilibrium picture) and then we define pair forces from the two-particle limit assuming a vanishing probability current. It could well be that the mapping in the first step breaks down parts of the many-body nature of the interactions in the active system, such that the effective attraction in the many-body system becomes accessible already on the level of pair interactions. As a logical next step, it seems worthwhile to study the effective potential extracted from computer simulations of a many-particle system, in order to clarify in how far the strong attraction of the effective potential needs to be seen as the result of a fortuitous cancellation of errors.

The low-activity limit in the effective equilibrium picture also results in pairwise forces. Under this assumption, we revealed some minor inconsistencies between the two equivalent steady-state conditions in Eqs. (9) and (12), although the latter contains a logarithmic term. Relatedly, it was recognized in Ref. (22) that different routes to define the active pressure only coincide at lowest order in the activity parameter . We suspect that further differences will occur at higher orders in and when employing further approximations, such as the EPA. We conclude that the route to follow should be carefully chosen for each problem, together with the underlying approximations.

The obvious purpose of both the low-activity approximation and the EPA is to allow for an analytically tractable theory. It appears that the condition given by Eq. (12) supported by effective pair potentials is most convenient for accessing structural properties (23); (12); (25), whereas the low-activity expansion of Eq. (9) provides a direct way to define mechanical properties (22). Moreover, our analysis suggests that the thermodynamic results obtained from Eq. (12) can be rescaled in order to obtain a workable definition of mechanical active pressure and surface tension. Arguably, the most simplistic scaling factor would be the diffusivity of an ideal gas, which can be absorbed into an effective temperature (21); (22); (25). A more general approach will be detailed in the second paper of this series.

## Acknowledgements

R. Wittmann, A. Scacchi and J. M. Brader acknowledge funding provided by the Swiss National Science Foundation. C. Maggi acknowledges support from the European Research Council under the European Unionâs Seventh Framework programme (FP7/2007-2013)/ERC Grant agreement No. 307940.

## Appendix A Simulation details

We performed Brownian Dynamics simulations of a system composed of two particles of unit diameter interacting through a soft-repulsive potential or a Gaussian soft-core potential. The potential is truncated at a distance of . In the simulations of active OUPs, evolving according to Eq. (1), each particle is subjected to Gaussian thermal noise and a non-Gaussian (colored) noise. The latter yields two distinct contributions to the displacement of each particle: one drift term, proportional to the reorientation time and one Gaussian process, proportional to  (26). For a vanishing active diffusivity , the drift term decays exponentially in time and is therefore irrelevant in the long-time limit. The integration time step is fixed to where is the time scale of translational diffusion. The total run time of the simulation is . For every , we calculate the distance between the two particles. The pair-correlation function is obtained in a standard way from the distance distribution. We have verified that, for the case of and finite , the obtained pair-correlation function is independent of , although the short-time displacement is not.

We also performed Brownian Dynamics simulation of ABPs, for which the colored-noise variable in Eq. (1) is replaced with the vector describing a constant velocity of the self-propulsion in the direction of the instantaneous orientation. The equation for the time evolution for the orientation vector of each particle is evaluated as an Ito integral, where is a white noise describing rotational diffusion. The integration time step is fixed to and the total run time is .

## Appendix B Effective many-body force in the Fox approximation

In this appendix we discuss the accuracy of Eq. (12) of the main text in the Fox approximation, i.e., choosing the effective diffusion tensor from Eq. (11). The accurate definition of the effective force is

 βFeffk=∑iD−1ik⋅βFi−∑ijD−1ik⋅∇j⋅Dji (28)

since the conversion

 ∑ijD−1ik⋅∇j⋅Dji≈∇kln|detD[N]| (29)

is only correct in the following cases:

1. for a system with no thermal noise. As stated in the main text, in this case the Fox and UCNA results are equivalent. Making use of the symmetry relation (with Greek indices labeling components and particles) and Jacobi’s formula the identity in Eq. (29) can be explicitly verified (21).

2. for a passive system, since .

3. at leading order in the activity parameter , where and the same arguments as under point 1. can be used.

4. for particles in an effectively one-dimensional symmetry, i.e., if there exists a coordinate frame in which the non-trivial contributions to reduce to an at most a tensor with identical diagonal elements. This can be easily shown by an explicit calculation

1. in dimensions

2. for a planar interaction potential

3. in the Laplacian approximation (18)

As a simple counterexample to the cases listed above, we note that Eq. (29) does not generally hold for and , since , whereas under point 4. we have . As for the approximate formulas (13) and (14), we find that the difference between the effective forces, Eq. (28), for particles with and without diagonal approximation is only a factor 2 in front of each factor .

Comparing the requirements for points 1. and 2. we can say that Eq. (12) is correct for both and , indicating that it should be a good approximation over all ranges of the parameter . Moreover, the assumption of a small persistence time is required in the derivation of the effective equilibrium approach (12); (20); (21). Considering point 3., this means that the approximation in Eq. (29) is consistent with the underlying theory. Indeed, Figs. 5a and 5b show that the approximation best for either small or large and small , respectively. In general, the difference is not significant compared to other approximations shown in Fig. 3 of the main text.

If the validity criterion for the Eigenvalues, given by Eq. (17), of is violated, neither side of Eq. (29) results in physical effective potentials. Therefore, the most important benefit of the approximate form on the right-hand side is that it allows to employ the inverse- approximation, as described in Sec. 3.2: the Eigenvalues of can be written as

 E(f)n(τ,Da,r)=1+DaEn(τ,r) (30)

so that we can substitute according to Eq. (20) of the main text. The substitution of or a more general manipulation of is inconvenient since the effective potential would still diverge for although the bare potential does not. Therefore, the third column of Fig. 3 contains the optimal implementation of the the inverse- approximation for the Fox approach. Also recall that, according to point 4.(c), the results in the Laplacian approximation are the same for both expressions in Eq. (29).

## Appendix C Integration of the first equality in Eq. (12)

The derivation of Eq. (21) by integrating Eq. (9) over coordinates is quite similar to that of the YBG hierarchy in a passive system. The first member

 0=∇ρ(r)+ρ(r)∇β\varv(r)+∫dr′ρ(2)(r,r′)∇βu(r,r′)=ρ(r)∇μ (31)

is recovered from (21) when setting . The second equality reflects the interpretation of the term on the left-hand side as the gradient of a chemical potential , which is constant in equilibrium. The second member reads

 0=∇ρ(2)(r,r′′)+ρ(2)(r,r′′)∇(β\varv(r)+βu(r,r′′))+∫dr′ρ(3)(r,r′,r′′)∇βu(r,r′)=ρ(2)