Effective Interactions in Active Brownian Suspensions
T.~F.~F. Farage
P. Krinninger
J.~M. Brader$^1$
Active colloids exhibit persistent motion, which can lead to motility-induced phase separation (MIPS). However, there currently exists no microscopic theory to account for this phenomenon. %which occurs entirely in the absence of attractive interactions. %arising from self-propulsion. %In sufficiently dense suspensions this motion becomes hindered by %interparticle interactions and a motility-induced phase separation can occur. %It has been suggested that this nonequilibrium phase transition %, which occurs entirely in the absence of attractive interactions, %could be %accounted for by mapping the active system onto an effective passive system, %interacting via an attractive pair potential. %However, a rational basis for such a mapping has so far been lacking and %no statistical mechanical theory exists for relating the microscopic properties of %individual particles to their collective behaviour. % We report a first-principles theory, free of fit parameters, for active spherical colloids, which shows explicitly how an {effective many-body interaction potential is generated by activity and how this can rationalize MIPS. %can arise %and how this is related to the passive pair potential. For a passively repulsive system the theory predicts phase separation and pair correlations in quantitative agreement with simulation. For an attractive system the theory shows that phase separation becomes suppressed by moderate activity, consistent with recent experiments and simulations, and suggests a mechanism for reentrant cluster formation at high activity.
Farage_Krinninger_Brader_manuscript.tex
\documentclass[pre,aps,twocolumn,superscriptaddress,showpacs]{revtex4}
%\documentclass[pre,aps,preprint,superscriptaddress,showpacs]{revtex4}

%% ADDITIONAL OPTIONAL STYLE FILES
\usepackage[dvips]{graphicx}
%\usepackage{pnastwoF}
\usepackage{amssymb,amsfonts,amsmath}
\usepackage{color}
\usepackage{ulem}


%% OPTIONAL MACRO DEFINITIONS
\newcommand{\sig}{\boldsymbol{\sigma}}
\newcommand{\kapp}{\boldsymbol{\kappa}}
\newcommand{\pdf}{\Psi(\{\textbf{r}_i\},t)}
\newcommand{\om}{\hat{\Omega}}
\newcommand{\ups}{\Upsilon}
\newcommand{\vt}{\vartheta}
\newcommand{\vp}{\varphi}
\newcommand{\X}{\boldsymbol{\chi}}

\newcommand{\omadj}{{\Omega}^{\dagger}}
\newcommand{\F}{\boldsymbol{F}}
\newcommand{\sigg}{\boldsymbol{\sigma}}
\newcommand{\siggh}{\hat{\boldsymbol{\sigma}}}
\newcommand{\B}{\textbf{B}}
\newcommand{\kk}{\boldsymbol{k}}
\newcommand{\qq}{\boldsymbol{q}}
\newcommand{\pp}{\boldsymbol{p}}
\newcommand{\J}{\boldsymbol{J}}
\newcommand{\E}{\textbf{E}}
\newcommand{\gam}{\dot{\gamma}}
\newcommand{\strain}{\boldsymbol{\varepsilon}}
\newcommand{\joe}{\color{red} }
\newcommand{\pvec}{\boldsymbol{p}}
\newcommand{\rr}{\boldsymbol{r}}
\newcommand{\nab}{\boldsymbol{\nabla}}
\newcommand{\xxi}{\boldsymbol{\xi}}
\newcommand{\xchi}{\boldsymbol{\chi}}
\newcommand{\eeta}{\boldsymbol{\eta}}
\newcommand{\R}{\boldsymbol{R}}
\newcommand{\Phik}{\Phi_{\boldsymbol{k}}}
\newcommand{\Gammak}{\Gamma_{\boldsymbol{k}}}
\newcommand{\mk}{m_{\boldsymbol{k}}}


\newcommand{\rhoq}{\rho_{\boldsymbol{q}}^{\phantom{\ast}}}
\newcommand{\rhok}{\rho_{\boldsymbol{k}}^{\phantom{\ast}}}
\newcommand{\rhop}{\rho_{\boldsymbol{p}}^{\phantom{\ast}}}
\newcommand{\rhol}{\rho_{\boldsymbol{l}}^{\phantom{\ast}}}
\newcommand{\rhosq}{\rho^{\ast}_{\boldsymbol{q}}}
\newcommand{\rhosk}{\rho^{\ast}_{\boldsymbol{k}}}
\newcommand{\rhosp}{\rho^{\ast}_{\boldsymbol{p}}}
\newcommand{\rhosl}{\rho^{\ast}_{\boldsymbol{l}}}

\newcommand{\rhoqq}{\rho_{\boldsymbol{q}'}^{\phantom{\ast}}}
\newcommand{\rhokk}{\rho_{\boldsymbol{k}'}^{\phantom{\ast}}}
\newcommand{\rhopp}{\rho_{\boldsymbol{p}'}^{\phantom{\ast}}}
\newcommand{\rholl}{\rho_{\boldsymbol{l}'}^{\phantom{\ast}}}
\newcommand{\rhosqq}{\rho^{\ast}_{\boldsymbol{q}'}}
\newcommand{\rhoskk}{\rho^{\ast}_{\boldsymbol{k}'}}
\newcommand{\rhospp}{\rho^{\ast}_{\boldsymbol{p}'}}
\newcommand{\rhosll}{\rho^{\ast}_{\boldsymbol{l}'}}


\begin{document}

\title{Effective Interactions in Active Brownian Suspensions}


\author{T.~F.~F. Farage}
\affiliation{Department of Physics,
  University of Fribourg, CH-1700 Fribourg, Switzerland}
\author{P. Krinninger}
\affiliation{Theoretische Physik II, Physikalisches Institut, 
  Universit{\"a}t Bayreuth, D-95440 Bayreuth, Germany}
\author{J.~M. Brader$^1$}







\begin{abstract}
Active colloids 
exhibit persistent motion, which can lead to motility-induced phase separation (MIPS). 
However, there currently exists no microscopic theory to account for this phenomenon.  
%which occurs entirely in the absence of attractive interactions.
%arising from self-propulsion. 
%In sufficiently dense suspensions this motion becomes hindered by 
%interparticle interactions and a motility-induced phase separation can occur. 
%It has been suggested that this nonequilibrium phase transition
%, which occurs entirely in the absence of attractive interactions, 
%could be 
%accounted for by mapping the active system onto an effective passive system, 
%interacting via an attractive pair potential. 
%However, a rational basis for such a mapping has so far been lacking and 
%no statistical mechanical theory exists for relating the microscopic properties of 
%individual particles to their collective behaviour. 
%
We report a first-principles theory, free of fit parameters, for active spherical colloids, 
which shows explicitly how an {\it effective} many-body interaction potential is generated 
by activity and how this can rationalize MIPS.  
%can arise 
%and how this is related to the passive pair potential. 
For a passively repulsive system the theory predicts 
phase separation and pair correlations in quantitative agreement with simulation. 
For an attractive system the theory shows that phase 
separation becomes suppressed by moderate activity, consistent with recent 
experiments and simulations, 
and suggests a mechanism for reentrant cluster formation at high activity. 
\end{abstract}
\pacs{82.70.Dd, 64.75.Xc, 05.40.-a}

\maketitle



Active colloidal particles in suspension are currently the 
subject of considerable attention, due largely to their ability to model 
self-organisation phenomena in biological systems, but also as a new branch 
of fundamental research in nonequilibrium statistical mechanics: 
assemblies of active colloids are intrinsically out-of-equilibrium systems. 
In contrast to their passive counterparts, active colloids undergo both solvent 
induced Brownian motion and a self-propulsion which requires a continual consumption of 
energy from the local environment.  
Several idealized experimental model systems have been developed, such as catalytic 
Janus particles \cite{lyderic2010,erbe2008,howse2007}, colloids with artificial 
flagella \cite{flagella2005} and light activated particles \cite{palacci2013}. 
The understanding of active systems has been further aided by the development of 
simple theoretical models, which aim to capture the essential physical mechanisms 
and which have been used to study e.g. bacteria, cells or filaments in the cytoskeleton
\cite{collective_motion2012, ramaswamy2010, romanczuk2012, cates_bacteria2012}.

Active particles are characterised by a persistent motion, which can lead to 
`self trapping' dynamics and a rich variety of related collective phenomena 
\cite{collective_motion2012, ramaswamy2010, romanczuk2012, cates_bacteria2012, 
cates_tailleur2014}. 
Even the simplest models of active spherical particles with purely repulsive 
interactions can display the phenomenon of motility-induced phase separation (MIPS) 
\cite{cates_tailleur2014}.  
In many respects, MIPS resembles the equilibrium phase separation familiar from 
passive systems with an attractive component to the interaction potential 
(e.g.~the Lennard-Jones 
potential) 
\cite{tailleur2008, fily2012, stenhammar2013, redner2013, levis_berthier2014}. 
This apparent similarity has motivated several recent attempts to map 
an assembly of active particles onto a passive equilibrium system, 
interacting via an effective attraction  
(usually taken to be a very short range sticky-sphere potential 
\cite{schwarz-linek2012,bocquetPRX_2014}). 
Despite the intuitive appeal of mapping to an equilibrium system, there 
exists no systematic theoretical approach capable of predicting an effective 
equilibrium potential directly from the bare interactions.	

Our current understanding of MIPS has largely been gained 
through either simulation 
\cite{fily2012,redner2013,stenhammar2013,levis_berthier2014,wysocki2014} 
or phenomenological theory 
\cite{tailleur2008,cates2013,stenhammar2013,cates_tailleur2014}. 
The phenomenological theory is based on an  
equation for the coarse-grained density, featuring a local speed  
and a local orientational relaxation time. 
Although the precise relationship between these one-body fields and the 
interparticle interaction potential remains to be clarified, some progress 
in this direction has been made \cite{bialke2013}.
%because the microscopic details of the coupling between interparticle 
%interactions, external fields and activity cannot be resolved.  
%%the relationship between the observed phenomena and the interaction potential
%e.g. reentrant phase behaviour of active Lennard-Jones particles
%\cite{schwarz-linek2012}, 
%%cannot be fully elucidated. 
On a more microscopic level, it has recently been shown that a general 
system of active particles
does not have an equation of state \cite{cates_kardar}, due to the influence of the 
confining boundaries, however, one can be recovered for the special case 
of active Brownian spheres \cite{cates_kardar,takatori}.  



%A higher level of complexity is reached when active particles move 
%in an external potential field. 
%There remains much to be understood about the general nature of the interaction between 
%persistent motion and external fields, although very recent work suggests that this 
%interaction may have profound implications for even the most basic 
%thermodynamic concepts (e.g. the equation of state) for certain classes of active system 
%\cite{cates_kardar}. 
%Progress has been made in treating noninteracting particles in certain special external fields. 
%Stark and coworkers have addressed the motion and distribution of active particles confined 
%in harmonic traps and find the emergence of orientational order \cite{enculescu2011,pototsky2012}. 
%One-dimensional, asymmetric periodic potentials are found to generate 
%rectified motion \cite{angelani2011,ai2013}. 


Here we report a first-principles theory for systems of active Brownian 
spheres, which demonstrates explicitly how an effective 
many-body interaction potential is induced by activity. 
An appealing feature of this approach is that intuition gained from equilibrium can be 
used to understand the steady-state properties of active systems.  
The required input quantities are the passive (`bare') interaction potential, 
the rotational diffusion coefficient and the particle propulsion speed. 
The theory generates as output the static correlation functions and 
phase behaviour of the active system. 
%involves neither empirical nor self-consistently determined 
%fit parameters. 
%The effective potential is then employed within the 
%framework of equilibrium statistical mechanics.  
For a repulsive bare interaction, activity generates an attractive 
effective pair potential, thus providing an intuitive explanation for the 
MIPS observed in simulations \cite{fily2012,levis_berthier2014,stenhammar2014}. 
For an attractive bare potential, we find that increasing activity 
first reduces the effective attraction, consistent with the experiments of 
Schwarz-Linek {\it et al.} \cite{schwarz-linek2012}, before leading at higher activity 
to the development of a repulsive potential barrier. 
We speculate that this barrier may be related to the reentrant phase behaviour 
observed in simulation by Redner {\it et al.} \cite{redner2013}. 
%by rendering the  
%homogeneous fluid unstable with respect to a cluster phase. 

%Finally, we extend our theory to treat 
%external fields and predict: (i) wetting of a repulsive 
%substrate by repulsive spheres, (ii) capilliary condensation of repulsive spheres 
%in a slit. 

The paper will be structured as follows: 
In \S\ref{theory} we specify the microscopic dynamics and 
describe how to eliminate orientational degrees of freedom. 
From the resulting coarse-grained, non-Markovian Langevin equation we 
derive a Fokker-Planck equation for the positional degrees of freedom, 
from which we identify an effective pair potential. 
%The liquid state integral equation 
%we employ to calculate pair correlations and the Brownian dynamics simulations 
%re specified. 
In \S\ref{results} we employ the effective pair potential in an equilibrium 
integral equation theory and investigate the structure and phase behaviour of 
both repulsive and 
attractive bare potentials. In the former case we predict MIPS, whereas 
in the latter case phase separation is suppressed by activity. 
Finally, in \S\ref{discussion} we discuss our findings and provide an outlook 
for future research. 


\section{Theory}\label{theory}
\subsection{Microscopic dynamics}  
We consider a three dimensional system of $N$ active, interacting, spherical Brownian 
particles with spatial coordinate 
$\rr_i$ and orientation specified by an embedded unit vector $\pvec_i$. 
Each particle experiences a self propulsion of speed $v_0$ in its direction of orientation.
Omitting hydrodynamic interactions the particle motion can be modelled by the 
overdamped Langevin equations
%
\begin{align}
&\dot{\rr}_i = v_0\,\pvec_i  + \gamma^{-1}\F_i + \xxi_i, 
\label{full_langevin}
\\
&\dot{\pvec}_i = \eeta_i\times\pvec_i ,
\label{langevin_orientation}
\end{align}
%	
where $\gamma$ is the friction coefficient and the force on particle $i$ is generated from the 
total potential energy according to $\F_i\!=\!-\nabla_i U_N$.
The stochastic vectors $\xxi_i(t)$ and $\eeta_i(t)$ are Gaussian distributed with zero mean and 
have time correlations	
$\langle\xxi_i(t)\xxi_j(t')\rangle=2D_t\boldsymbol{1}\delta_{ij}\delta(t-t')$ and 
$\langle\eeta_i(t)\eeta_j(t')\rangle=2D_r\boldsymbol{1}\delta_{ij}\delta(t-t')$, where 
$D_t$ and $D_r$ are the translational and rotational diffusion coefficients. 

Equations \eqref{full_langevin} and \eqref{langevin_orientation} are convenient for simulation, 
but are perhaps not the most suitable starting point for developing a first-principles microscopic 
theory. 
For a homogeneous system, averaging over the angular degrees of freedom generates a 
coarse-grained equation \cite{fily2012} 
%
\begin{align}\label{effective_langevin}
\dot{\rr}_i(t) &= \gamma^{-1}\F_i(t) + \xxi_i(t) + \xchi_i(t), 
\end{align}
%
where $\xchi_i(t)$ is a Markov process with zero mean and where the
time correlation function is given by
%
\begin{align}\label{exp_correlation}
\langle\xchi_i(t)\,\xchi_j(t')\rangle = \frac{v_0^2}{3} e^{-2D_r|t-t'|}\boldsymbol{1}\delta_{ij}.   
\end{align} 
%
The average in \eqref{exp_correlation} is over both noise and initial orientation. 
The distribution of $\xchi_i(t)$ is Gaussian to a good approximation.
This point and further technical details of the coarse graining are discussed 
in Appendix \ref{appendix1}. 
Equation \eqref{effective_langevin} provides a mean-field level of description, which 
deviates from the exact equations \eqref{full_langevin} and \eqref{langevin_orientation} 
by neglecting the coupling of fluctuations in orientation and positional degrees 
of freedom. 

The Langevin equation \eqref{effective_langevin} describes a non-Markovian 
process, which approximates the stochastic time evolution of the 
positional degrees of freedom.  
%%is statistically equivalent to 
%equations \eqref{full_langevin} and
The persistent motion of active particles is here encoded by 
the exponential decay of the time correlation \eqref{exp_correlation}, with persistence 
time $\tau_p\!=\!(2D_r)^{-1}$\!. 
For small $\tau_p$ the time correlation becomes 
$\langle\xchi_i(t)\xchi_j(t')\rangle = 2D_a\boldsymbol{1}\delta_{ij}\delta(t-t')$ and 
the dynamics reduce to that of an equilibrium system 
with diffusion coefficient $D_t+D_a$, where $D_a\!=\!v_0^2/(6D_r)$. 
This limit is realized when $\tau_p$ is shorter 
than the mean free time between collisions, i.e. in a dilute suspension.  
To treat finite densities requires an approach which deals with 
persistent trajectories. 
With this aim, we adopt \eqref{effective_langevin} as the starting point for contructing  
a closed theory. 

\subsection{Fokker-Planck equation}
A stochastic process driven by colored noise, such as that described by equation 
\eqref{effective_langevin}, is always non-Markovian. 
Consequently, it is not possible to derive an {\it exact} Fokker-Planck equation 
for the time evolution of the probability distribution \cite{vankampen1976}.  
Nevertheless, an approximate Fokker-Planck description capable of making accurate predictions 
can usually be found. 
The approximate Fokker-Planck equation implicitly defines a Markov
process which best approximates the process of physical interest (although precisely what 
constitutes the `best' approximation remains a matter of debate). 
From the extensive literature on this subject (see \cite{vankampen1976,vankampen1998,faetti1988} 
and references therein) has emerged a powerful method due to Fox \cite{fox1986,fox_weak1986}, 
in which a perturbative expansion in powers of correlation time is partially resummed using 
functional calculus. 
The resulting Fokker-Planck equation is most accurate for short correlation times 
(`off white' noise \cite{vankampen1998}) and for one-dimensional models 
makes predictions in good agreement with simulation data \cite{faetti1988}.

%\subsection{Many-body forces}
%Initially restricting our attention to homogeneous systems 
We now consider applying the method of Fox \cite{fox1986,fox_weak1986} to 
equation \eqref{effective_langevin}. 
This approach consists of first formulating the configurational probability distribution as 
a path (functional) integral and then making a time-local, Markovian approximation to this quantity. 
Technical details of the method are given in Appendix \ref{appendix2}.
Fox's approach was originally developed to treat one-dimensional problems 
\cite{fox1986,fox_weak1986}, however the generalization to three dimensions is 
quite straightforward.  
This enables us to directly obtain the following Fokker-Planck equation  
%
\begin{align}
	\partial_t\Psi(\rr^N\!\!,t) = -\sum_{i=1}^{N}\nabla_i\cdot \J_i(\rr^N\!\!,t) ,
\label{smol_eq}
\end{align}
%
where $\Psi(\rr^N\!\!,t)$ is the configurational probability distribution.  
Within the generalized Fox approximation the many-body current is given by
%
\begin{align}\label{current_i}
	\J_i(\rr^N\!\!,t) = -D_i(\rr^N)\left[\,\nabla_i 
-\beta \F_i^{\rm eff}(\rr^N)\,\right] \Psi(\rr^N\!\!,t) ,
\end{align}
%
where $\beta\!\equiv\!(k_BT)^{-1}$. 
The diffusion coefficient is given by 
%
\begin{align}\label{diffusion_i}
D_i(\rr^N) = D_t + D_a\left(
1 +
\frac{\tau \nabla_i\!\cdot\!\beta\F_i(\rr^N)}{1 - \tau \nabla_i\!\cdot\!\beta\F_i(\rr^N)}
\right), 
\end{align}
%
where we have defined a dimensionless persistence time, $\tau\!=\!\tau_p D_t/d^2$.  
The effective force is given by
%
\begin{align}
\F_i^{\rm eff}(\rr^N) =& 
 \frac{1}{\mathcal{D}_i(\rr^N)}\Big( 
\F_i(\rr^N) - k_BT\,\nabla_i \mathcal{D}_i(\rr^N)
\Big),
\label{force_i}
\end{align}
%
where $\mathcal{D}_i(\rr^N)\!=\!D_i(\rr^N)/D_t$ is a dimensionless diffusion coefficient. 
Either in the absence of interactions or in limit of large $D_r$ 
the diffusivity \eqref{diffusion_i} reduces to
$D_t\!+\!D_a$ and the effective force becomes 
$D_t \F_i(\rr^N)/(D_t+D_a)$. 
In this diffusion limit the system behaves as an equilibrium system at 
effective temperature $T_{\rm eff}=T(1+D_a/D_t)$. 

For weakly persistent motion, $\tau\!\rightarrow\! 0$, equations 
\eqref{smol_eq} to \eqref{force_i} become exact and the theory provides 
the leading order correction to the diffusion approximation. 
However, the Fox approximation goes beyond this by including contributions 
to all orders in $\tau$. Indeed, detailed studies of one-dimensional systems 
have demonstrated 
good results over a large range of $\tau$ values \cite{faetti1988}. 
%The Fox approach has the property of uniform convergence \cite{fox_weak1986} 
The only caveat is that the condition $1\!-\!\tau\, \nabla_i\!\cdot\!\beta\F_i\!>\!0$  
must be satisfied \cite{fox1986,fox_weak1986}. The range of accessible $\tau$ 
values thus depends upon the specific form of the bare interaction potential.



Within our stochastic calculus approach, the effective many-body force 
\eqref{force_i} emerges in a natural way from the coarse grained Langevin equation 
\eqref{effective_langevin}.
The more standard route (adopted in all attempts made so 
far \cite{pototsky2012,bialke2013}) to approach this problem is to derive from the 
Markovian equations \eqref{full_langevin} the exact Fokker-Planck 
equation for the joint distribution of positions and orientations, $P(\rr^N\!\!,\pvec^N\!\!,t)$. 
However, coarse graining strategies based on integration of $P$ over 
orientations generate intractable integral terms.  
By starting from \eqref{effective_langevin} we are able to circumvent these difficulties. 
As we shall demonstrate below, our effective force accounts for several important 
collective phenomena in active systems. 


\subsection{Effective pair potential}
In the low density limit we need only consider isolated pairs of particles. 
In this limit \eqref{smol_eq} reduces to an equation of motion for the radial 
distribution function, $g(r,t)\!\equiv\!\Psi(r,t)/\rho_b^2$, where 
$\rho_b$ is the bulk density. 
This equation of motion, the pair Smolochowski equation, is given by
%
\begin{align}
\partial_tg(r,t)\!=\!-\nabla\cdot {\bf j}(r,t), 
\end{align}
%
where $r\!=\!|\rr_{12}|$ is the particle separation and $\nabla\!=\!\nabla_{\rr_{\!12}}$.   
The pair current is given by
%
\begin{align}\label{current_pair}
	{\bf j}(r,t) = -2 D(r) g(r,t)\left[\,\nabla\ln g(r,t)
- \beta \F^{\rm eff}(r)\,\right].
\end{align}
%
where the radial diffusivity 
%
\begin{align}\label{diffusion_pair}
D(r) = D_t + D_a\left(
1 -
\frac{\tau \nabla^2\beta u(r)}{1 + \tau \nabla^2 \beta u(r)}
\right), 
\end{align}
%
interpolates between the value $D_t$ at small separations, where $u(r)$ is strongly repulsive, 
and $D_t+D_a$ at large separations. 
The effective interparticle force is given by
%
\begin{align}
\F^{\rm eff}(r) = 
 \frac{1}{\mathcal{D}(r)}\Big( 
\F(r) - k_BT\,\nabla\mathcal{D}(r)
\Big),
\label{force_pair}
\end{align}
where the bare force is related to the pair potential by 
$\F(r)\!=\!-\nabla u(r)$. 
The symmetry of the two-body problem can be exploited to calculate from 
\eqref{force_pair} an effective interaction 
potential	 
\begin{align}
\beta u^{\rm eff}(r)=\int_r^{\infty}\!\!dr' \left( \frac{\beta F(r')}{\mathcal{D}(r')} 
- \frac{\partial}{\partial r'}\ln \mathcal{D}(r') \right),  
\label{potential_pair}
\end{align}
where $F(r)\!=\!|\F(r)|$. 
We have thus identified an effective interaction pair potential, which requires as input 
the bare potential and the activity parameters $\tau$ and $D_a$. 

\begin{figure}[t!]
\hspace*{-0.1cm}\includegraphics[width=8.65cm]{figure_one.pdf}
\vspace*{0.2cm}
\caption{\label{r12_fig}{\bf Activity induces effective attraction.}
Passive potential $\beta u(r)\!=\!r^{-12}\!$.
{\bf(a)} Increasing $Pe$ (in steps of $8$) from $0$ to $40$ 
generates an effective interparticle attraction. 
Points indicate the potential minima. 
{\bf(b)} Radial distribution function, $g(r)$, from simulation (points) and theory (lines) 
for $\rho_b\!=\!0.5$ and $Pe\!=\!0$ to $20$ (in steps of $4$). Curves are shifted vertically 
for clarity.  
{\bf(c)} As in (b), but focusing on larger separations for $Pe=4$ (squares), 
$12$ (circles) and $20$ (diamonds). Inset: position of the first peak in $g(r)$ as a function 
of $Pe$. 
{\bf(d)} Spinodals for $\tau=0.045$ (dot-dashed) to $0.065$ (long dashed) in steps of 
$0.005$.   
} 
\end{figure}


\section{Results}\label{results}
\subsection{Motility-induced phase separation (MIPS)}
To illustrate how activity can generate an effective attraction in a passively repulsive system 
we consider the non-specific potential $\beta u(r)=r^{-12}$. 
In Fig.\ref{r12_fig}a we show the evolution of the effective potential \eqref{potential_pair} 
for fixed $\tau$ as a function of the dimensionless velocity $Pe\!=\!v_0d/D_t$. 
For $Pe\!\gtrsim\!10$ the effective potential develops an attractive tail. 
As $Pe$ is increased the potential well deepens, the minimum moves to 
smaller separations and the radius of the soft repulsive core decreases. 
These trends are consistent with the intuitive picture that persistent motion drives soft 
particles into one another (the soft core radius reduces) and 
that they remain dynamically coupled (`trapped') for longer than in the corresponding 
passive system. 
Within our equilibrium picture the trapping is accounted for by the effective attraction. 
%The effective potential thus displays behaviour in line with physical intuition.

For systems at finite density the pair potential \eqref{potential_pair} is an approximation, because 
three- and higher-body interactions will play a role (see equation \eqref{force_i}). 
However, for simplicity we henceforth employ the pair potential \eqref{potential_pair} for 
all calculations, as we anticipate that this will provide the dominant contribution. 
Although corrections to this assumption can be made, they obscure the physical picture and 
come at the expense of a more complicated theory. 
The validity of the pair potential approximation is justified {\it a posteriori} by 
the comparison with simulation for the finite density pair correlations.

In Fig.\ref{r12_fig}b we show the steady-state (isotropic) radial distribution 
function for $\rho_b\!=\!0.5$ for various values of $Pe$.
We employ the effective pair potential \eqref{potential_pair} together with 
liquid state integral equation theory and compare theoretical predictions
%We assume that the average current \eqref{current_i} vanishes, such that 
%the effective potential \eqref{potential_pair} can be input to an equilibrium calculation. 
%The equilibrium theory we employ to calculate $g(r)$ is the soft-core mean spherical approximation 
%(SMSA) proposed by Madden and Rice \cite{madden1980}, which is known to provide reliable and 
%accurate results for the structure and thermodynamics of Lennard-Jones-type potentials. 
with direct Brownian dynamics simulation of equations \eqref{full_langevin} and 
\eqref{langevin_orientation}. 
The integral equation theory we employ is the soft mean-spherical approximation (SMSA) 
proposed by Madden and Rice \cite{madden1980}.   
This approximate closure of the Ornstein-Zernike equation is known to provide reliable 
results for the pair stucture of Lennard-Jones type potentials.
Given the form of the effective pair potential shown in Fig.\ref{r12_fig} the SMSA 
would seem to be a reasonable choice of closure. 
Details of the integral equation theory and the simulation procedure are given in Appendices 
\ref{integralequation} and \ref{simulations}, respectively.  

\begin{figure}[t!]
\hspace*{-0.1cm}\includegraphics[width=8.65cm]{figure_two.pdf}
\vspace*{0.2cm}
\caption{\label{figure2}{\bf From active suppression of phase separation to a cluster phase.} 
Passive potential $\beta u(r)\!=\!4\,\varepsilon(r^{-12}\!-r^{-6})$ with $\varepsilon\!=\!1.4$. 
{\bf (a)} Effective potential \eqref{potential_pair} for $\tau\!=\!0.025$ and 
$Pe\!=\!0$ (black) $4, 8, 12$ (red broken lines) and $20, 28, 36$ (blue).
Inset: zoom of the repulsive peak for $Pe\!=\!20, 28, 36$.
{\bf (b)} Spinodals for $Pe\!=\!0, 4, 8, 12$. Increasing $Pe$ increases $\varepsilon_{\rm crit}$, 
the critical value of $\varepsilon$. 
{\bf (c)} $\varepsilon_{\rm crit}$ as a function of $Pe$ (black full line) 
and the locus of points for which the repulsive peak of $\beta u^{\rm eff}$ 
takes the value $0.1$ (red dashed) and $0.05$ (blue dot-dashed). 
Open (closed) circles indicate points where BD simulation find a mixed (phase-separated) 
state (see figure \ref{sim_snapshot}).
Arrow indicates path taken in (a). 
{\bf (d)} Theory (lines) and simulation (symbols) data (shifted for clarity) 
for $g(r)$ at $\varepsilon=0.5$, $\rho_b=0.3$ for $Pe=0$ (black), $12$ (blue) and $20$ (green). 
Dotted lines indicate peak positions.   
} 
\end{figure}

We find that as $Pe$ is increased the main peak of $g(r)$ grows in height and 
shifts to smaller separations (see 
inset to Fig.\ref{r12_fig}c), reflecting the changes in the effective potential. 
%This is consistent with the simulated structure factors of Ref.\cite{ni2013}. 
%who find 
%that the main peak of the static structure factor shifts to larger wavevector with 
%increasing $Pe$. 
In the main panel of Fig.\ref{r12_fig}c we focus on the second and third peaks. 
The quantitative accuracy of the theory in describing the decay of $g(r)$ 
is quite striking, in particular the phase shift induced by increasing activity 
is very well described. 
Further comparison for other parameter values (not shown) suggests that 
\eqref{potential_pair} combined with the SMSA theory 
provides an accurate account of the asymptotic decay of pair correlations.    

In Fig.\ref{r12_fig}d we show the spinodal lines mapping the locus of points for which the 
static structure factor, $S(q)=(1-\rho_b c(q))^{-1}$, diverges at vanishing wavevector. 
Simulations have shown that MIPS 
is consistent with a spinodal instability \cite{fily2012}. 
%Systems with parameters within the spinodal region are mechanically unstable with respect to density 
%fluctuations. 
%and will undergo phase separation. 
As $\tau$ is decreased the critical point moves to higher values of $Pe$ and to slightly higher densities. 
When compared with the spinodal of a standard Lennard-Jones system (e.g. the 
black curve in Fig.\ref{figure2}b) 
the critical points in Fig.\ref{r12_fig}d lie at rather higher 
values of $\rho_b$. 
This suggests that typical coexisting liquid densities for MIPS will be larger than those found  
in equilibrium phase separated systems, as has been observed in simulation 
\cite{fily2012,redner2013}.


\subsection{Suppression of phase separation} 
We next consider the influence of activity on a Lennard-Jones system, 
$\beta u(r)\!=\!4 \varepsilon (r^{-12}\!-r^{-6})$. 
For a phase-separated passive system, recent experiments  
and simulations have demonstrated that increasing $Pe$ first suppresses 
the phase separation \cite{schwarz-linek2012} and then leads at higher $Pe$ to a reentrant MIPS 
\cite{redner2013}. 
Schwarz-Linek {\it et al.} have argued that the suppression of phase separation at lower to intermediate 
$Pe$ occurs in their system because particle pairs bound by the attractive (depletion) 
potential begin to actively escape the potential well, and that this can be mimicked 
using an effective potential less attractive and shorter-ranged than the bare 
potential \cite{schwarz-linek2012}. 
%However, no direct test of this assumption was made.  



To investigate these phenomena we set $\varepsilon\!=\!1.4$, which ensures a phase separated 
passive state \cite{barker_henderson1976}, and consider the evolution of the 
effective potential as a function of $Pe$. 
In Fig.\ref{figure2}a we show that as $Pe$ is increased from zero to the value $18$ 
both the depth and range 
of the effective potential reduce significantly, consistent with the expectation of 
Schwarz-Linek {\it et al.} 
\cite{schwarz-linek2012}. 
Spinodals within this range of $Pe$ values, identifying where the static structure factor diverges 
at zero wavevector, are shown in Fig.\ref{figure2}b. 
As $Pe$ is increased the critical point moves to higher values of $\varepsilon$ 
(cf. Figs.1 and 3 in Ref.\cite{schwarz-linek2012}). 
A passively phase-separated system will thus revert to a single phase upon increasing the activity. 
To examine this behaviour in more detail we show in Fig.\ref{figure2}c the trajectory of the 
critical point in the ($Pe,\varepsilon$) plane. 
Above the line there bulk densities for which phase separation occurs. 

%
%
\begin{figure}[!t]
\hspace*{0.2cm}\includegraphics[width=8.65cm]{figure_three.pdf}
\vspace*{-0.2cm}
\caption{{\bf Simulated phase separation}. {\bf (a)} Snapshot of a mixed system at 
$t/\tau_B = 40$, $Pe = 8$, $\varepsilon /(k_B T)=1.5$. {\bf (b)} Snapshot of a phase 
separating system at the same time and $Pe=8$, $\varepsilon=2.5$. {\bf (c)} The radial distribution function, 
$g(r)$, for $\varepsilon=1.5$ (black curve), 2 (blue curve), and 2.5 (green 
curve). {\bf (d)} As in (c) but focusing on larger distances. From the 
snapshots together with the long range behaviour of the $g(r)$ we can 
distinguish between a mixed and a phase separated system. The slow decay to the asymptotic value of unity, as shown in (d), indicates phase separation.   
} 
\label{sim_snapshot}
\end{figure}
%
%

In order to test the predicted trajectory of the critical point we have performed Brownian 
dynamics simulations at a bulk density $\rho_b\!=\!0.4$, which lies close to the critical 
density \cite{barker_henderson1976}, 
for various values of $\varepsilon$ and the $Pe$ values $5, 8$ and $10$. 
Visual inspection of the simulation snapshots reveals the existence of voids in the particle 
configurations corresponding to a phase separated state 
(see Figs.\ref{sim_snapshot}a and \ref{sim_snapshot}b for snapshots). 
%(the runtime must be sufficiently long for these voids to develop.) 
This visual impression can be made more quantitative by calculating the radial 
distribution function. Phase separating states generate a very characteristic 
slow decay of $g(r)$ (Figs.\ref{sim_snapshot}c and \ref{sim_snapshot}d) which 
provides a useful indicator. 
The open circles in Fig.\ref{figure2}c represent mixed states, whereas closed circles 
indicate phase separated statepoints. 
The phase boundary predicted by the theory is highly consistent with the simulation data. 
%thus providing evidence that using the effective potential \eqref{potential_pair} in an 
%equilibrium theory captures the correct physics.  



\subsection{Cluster phase}    
Returning to Fig.\ref{figure2}a, we find that for $Pe\!>\!18$ the effective 
potential develops a repulsive barrier, which grows in height (see inset) with increasing 
$Pe$, while the potential minimum becomes deeper. 
It is well-known that potentials with a short-ranged attraction and long-ranged 
repulsion (SALR potentials) exhibit unusual equilibrium phase behaviour, including clustering 
and microphase-separation \cite{sear1997,archer2007}. 
Although the attractive component of the potential may favour phase separation, the long range repulsion 
destabilises distinct liquid and gas phases and causes them to break up into droplets or clusters. 
This represents a non-spinodal type of phase transition, characterized by a divergence 
in the structure factor at finite wavevector. 
The appearance of a repulsive barrier in the effective potential suggests that a similar 
mechanism may be at work in passively attractive systems subject to high $Pe$ activity. 


In Fig.\ref{figure2}c we show the locus of points where the effective potential peak height attains 
a given value (we choose $0.05$ and $0.1$ for illustration). 
When these `isorepulsion curves' are viewed together with the critical point 
trajectory the resulting phase diagram is very similar to that obtained by Redner 
{\it et al.} in their simulation study of two-dimensional active Lennard-Jones particles 
(cf. figure 1 in Ref.\cite{redner2013}). 
However, a detailed study of the connection between the potential barrier and 
high $Pe$ clustering goes beyond the scope of the present work. 

%{\joe
%Although intuition suggests that reentrant clustering arises primarily from the interaction of 
%persistent motion with passive repulsion,  
%our theory suggests that the underlying mechanism may be more complicated.}

\section{Discussion}\label{discussion}  
In summary, we have shown that systems of active spherical Brownian particles 
can be mapped onto an equilibrium system interacting 
via an effective, activity-dependent many-body potential. 
The only required inputs are the bare potential, thermodynamic statepoint 
and the parameters specifying the state of activity. 
%Our strategy is to first reduce the defining  Langevin equations \eqref{full_langevin} 
%to an exact, coarse-grained Langevin equation for the particle positions, thus 
%integrating-out orientational degrees of freedom.  
%We then derive an approximate Fokker-Planck equation, from which an effective 
%potential can be inferred. 
Our theory captures the phenomenon of MIPS in repulsive systems and 
provides first-principles predictions for the activity 
dependence of the pair correlations, in very good agreement with 
Brownian dynamics simulation. 
As far as we are aware there exists no other approach capable of predicting 
from the microscopic interactions the pair correlation functions of an active system. 
Further insight into the steady state particle distribution could in principle be obtained 
by investigating the three-body correlations. 
These could be obtained by employing the effective 
potential in a higher-order liquid state integral equation theory 
(see e.g.~\cite{brader_triplet} and references therein).  


For passively attractive systems the theory rationalizes the experimental 
finding \cite{schwarz-linek2012} that increasing activity can suppress passive phase separation.  
We find that as $Pe$ is increased from zero to intermediate values the minimum of the effective 
potential becomes less deep, thus weakening the cohesion of the liquid phase. 
To the best of our  knowledge there is currently no alternative theoretical explanation 
of phase-transition suppression in active suspensions. 
It is an appealing aspect of our theory that the suppression of passive phase separation 
follows naturally from the same approach which yields activity-induced attraction for 
repulsive potentials. 
For high values of $Pe$ the appearance of a repulsive barrier in the effective 
potential suggests that the reentrant phase separation observed in simulations 
\cite{redner2013} may be interpreted using concepts of equilibrium 
clustering in SALR potential systems. 
This will be a subject of future detailed investigations. It is known that 
care must be exercised when analysing SALR potentials, as 
traditional liquid state theories can prove misleading \cite{archer2007}.


%
%For inhomogeneous systems we have generalized our approch to incorporate an effective external 
%potental, which breaks the translational symmetry. 
%Here we find that a passively repulsive system will exhibit both wetting and capillary 
%condensation as a function of increasing activity. Although the mechanisms identified appear 
%to be quite robust, it will be of great interest to see whether future simulation or experimental studies 
%can confirm the theoretical predictions. The liquid-vapour interface is found to be similar 
%to that in equilibrium, with an interfacial width which depends on $Pe$. 
%
%Now that we have derived an explicit mapping of a class of active systems onto an effective 
%equilibrium 

A key step in our development is the Fox approximation \cite{fox1986}, which yields an 
effective Markovian description of the coarse-grained equation \eqref{effective_langevin}.  
Making a Markovian approximation automatically imposes an effective equilibrium,  
however, we are aware that there exist certain situations for which this breaks down 
\cite{cates2013,cates_tailleur2014}. 
Establishing more clearly the range of validity of our approach, as well as its possible 
extensions, will be the subject 
of ongoing study. However, it is already clear that going beyond the Markovian approximation 
will be very challenging. 
Indeed, such a step may not even be desirable. 
Any kind of non-Markovian description 
would lead inevitably to a loss of the effective equilibrium picture and the physical 
intuition associated with it. 
It thus seems likely that 
practical improvements to the present approach will retain the Markovian 
description while seeking to optimize, or improve upon, the Fox approximation for certain 
classes of bare potential. Very recently, Maggi \textit{et al.} have employed 
an alternative approach to treating stochastic processes driven 
by Ornstein-Uhlenbeck noise \cite{maggi2015}. A comparison of their approach 
with the Fox method employed here would be very interesting.

With a view to further applications of our approach, we note that there has 
recently been considerable interest in active suspensions at very high 
densities 
\cite{faragebrader_arxiv,ni2013,berthier2013,szamel_berthier2015}. 
In particular, it has been found using computer simulations that activity has a 
strong influence on the location of the hard-sphere glass 
transition, dynamic correlation functions, such as the intermediate scattering 
function, and static pair correlations \cite{ni2013}. 
Within our effective equilibrium framework, increasing the activity of a passively 
repulsive system generates an effective attraction. 
We can therefore anticipate that for volume fractions just above the glass transition 
it will be possible to observe a reentrant glass transition, namely a melting of 
the glass followed by revitrification, as a function of increasing $Pe$. 
Moreover, the nontrivial evolution of the effective potential as a function of 
$Pe$ for attractive bare potentials (cf.~Fig.\ref{figure2}a) 
suggests these systems will present a rich variety of glassy states. 
Work along these lines is in progress.

Finally, we mention that a natural generalization of the present theory is to 
treat spatially inhomogeneous systems in external fields.  
Recent microscopic studies of active particles under confinement 
(e.g.~in a harmonic trap \cite{pototsky2012}) have provided considerable insight, 
however none of the existing approaches have considered effective {\it interparticle} 
interactions. 
Inhomogeneous generalization of the present theory enables the interaction between 
MIPS and external fields to be investigated on the microscopic level. 
Our preliminary investigations reveal, for example, activity-induced wetting at 
a planar substrate and capillary-condensation under confinement. 
This will be presented in a future publication.  



\begin{acknowledgments}
We thank Yaouen Fily, Ronald Fox and Paolo Grigolini for helpful correspondence. 
Funding provided by the Swiss National Science Foundation. 
PK thanks the Elitenetzwerk Bayern (ENB) and Matthias Schmidt for 
financial support.  
\end{acknowledgments}



\appendix
\section{Coarse-grained Langevin equation}\label{appendix1}
Equation \eqref{langevin_orientation} describes the orientational diffusion 
of an active particle. 
%\begin{equation}
% \frac{d}{dt}\pp(t)=\eeta(t)\times\pp(t)\,,
%\end{equation}
%(the second of Eqs.1 in the main text)
%where $\eeta(t)$ is gaussian distributed with $\langle\eeta(t)\rangle=0$ and 
%$\langle\eeta(t)\eeta(t')\rangle=2D_r\boldsymbol{1}\delta(t-t')$. 
%
The corresponding conditional probability distribution function $\ups(\pp,t\,|\,\pp_0, 
t_0)$, where $t>t_0$, obeys a Fokker-Planck equation which can be obtained by 
usual techniques \cite{gardiner},
%
\begin{equation}\label{FP_p}
 \frac{\partial}{\partial t}\ups(\pp,t\,|\,\pp_0, t_0) =D_r\R^2\ups(\pp,t\,|\,\pp_0, t_0)\,,
\end{equation}
%
where $\R\equiv(\pp\times\nabla_{\pp})$ is the intrinsic angular momentum 
differential operator. Eq.\eqref{FP_p} describes nothing but a diffusion 
process on the unit sphere. This problem is well-known when studying, 
e.g., dielectric relaxation in polar liquids 
\cite{debye1929,fatuzzo1967,zwanzig1970,berne1975}. In spherical coordinates, 
\eqref{FP_p} becomes
\begin{align}\label{FP_spher}
 \frac{1}{D_r}\frac{\partial}{\partial t}\ups(\Omega,t\,|\,\Omega_0,t_0) &=
 \left[\frac{1}{\sin \vt}\frac{\partial}{\partial 
\vt}\left(\sin \vt\frac{\partial}{\partial 
\vt}\right)\right.\notag\\
&\phantom{=}\left.
+\frac{1}{\sin^2\vt}\frac{\partial^2}{\partial\vp^2}\right]
\ups(\Omega,t\,|\,\Omega_0,t_0),
\end{align}
where we have defined $\Omega\equiv(\vt,\vp)$. 
%We recall that spherical 
%harmonics can be expressed by
%\begin{equation}\label{Ylm}
% Y_{lm} = 
%(-1)^m\sqrt{\frac{(2l+1)(l-m)!}{4\pi(l+m)!}} e^{im\vp}P_{lm}(\cos{\vt}),
%\end{equation}
%where $l\in\mathbb{N}$ and $m\in\mathbb{Z}$, such that $-l\leq m \leq l$, 
%and the associated Legendre polynomials are given by 
%\begin{equation}\label{legendre}
%P_{lm}(x)=\frac{(-1)^m}{2^ll!}(1-x^2)^{m/2}\frac{d^{l+m}}{dx^{l+m}}(x^2-1)^l.
%\end{equation}


%Because spherical harmonics form a complete set of orthonormal functions we 
%can expand the probability distribution function $\ups$ in terms of the 
%$Y_{lm}$, which then implicitly implies that we 
Assuming that $\ups$ and its derivatives are continuous on the sphere \cite{arfken2013}, 
we expand the probability distribution function $\ups$ in spherical harmonics
\begin{equation}\label{solexp}
\ups(\Omega,t\,|\,\Omega_{0},t_0)=\sum_{l=0}^{\infty}\sum_{m=-l}^{l}A_{lm}
(t\,|\,\Omega_0,t_0)Y_{lm}(\Omega)\,,
\end{equation}
where $Y_{lm}$ are the spherical harmonics and $A_{lm}$ are coefficients 
encoding the initial condition.
We also recall that spherical harmonics are 
eigenvectors of the operator 
$\R^2$ (in spherical coordinates), namely that
\begin{equation}\label{eigenvectors}
 \R^2 Y_{lm}=-l(l+1)Y_{lm}\,.
\end{equation}
Inserting \eqref{solexp} in \eqref{FP_spher} and using \eqref{eigenvectors} 
we obtain
\begin{align}\label{prev_orth}
 \sum_{l,m}\frac{\partial}{\partial 
t}&A_{lm}(t\,|\,\Omega_0,t_0)Y_{lm}(\Omega) =\notag\\
&-D_r\sum_{l,m}
l(l+1)A_{lm}(t\,|\,\Omega_0,t_0)Y_{lm}(\Omega).
\end{align}
Multiplying both sides of \eqref{prev_orth} by $Y_{l'm'}^{\ast}(\Omega)$, integrating 
over solid angle and using the orthogonality property, 
$\int d\Omega\, 
Y_{l'm'}^{*}(\Omega)Y_{lm}(\Omega)=\delta_{m,m'}\delta_{l,l'}$, yields
\begin{equation}
 \frac{\partial}{\partial 
t}A_{lm}(t\,|\,\Omega_0,t_0)=-D_rl(l+1)A_{lm}(t\,|\,\Omega_0,t_0)\,,
\end{equation}
which has the solution 
\begin{equation}
 A_{lm}(t\,|\,\Omega_0,t_0)=e^{-D_rl(l+1)(t-t_0)}a_{lm}(\Omega_0)\,,
\end{equation}
where the $a_{lm}$ are a new set of coefficients. 
The probability distribution is thus given by
\begin{equation} 
\ups(\Omega,t\,|\,\Omega_0,t_0)=\sum_{l,m}e^{-D_rl(l+1)(t-t_0)}a_{lm}(\Omega_0)Y_{lm
}(\Omega).
\end{equation}
The initial condition,
\begin{equation}
 \ups(\Omega,t_0\,|\,\Omega_0,t_0)=\delta(\Omega-\Omega_0)\,,
\end{equation}
together with the completeness relation of the spherical harmonics,
\begin{equation}
\delta(\Omega-\Omega_0)=\sum_{l=0}^{\infty}\sum_{m=-l}^{l}Y_{lm}(\Omega)Y_{l
m}^{\ast}(\Omega_0),
\end{equation}
allows the missing coefficients to be identified,
\begin{equation}
 a_{lm}(\Omega_0)=Y_{lm}^{\ast}(\Omega_0)\,.
\end{equation}
The conditional probability distribution is now fully determined as
\begin{equation}\label{sol_pdf}
\ups(\Omega,t\,|\,\Omega_{0},t_0)=\sum_{l=0}^{\infty}\sum_{m=-l}^{l}e^{
-D_rl(l+1)(t-t_0)}Y_{lm}^{\ast}(\Omega_0)Y_{lm}(\Omega).
\end{equation}
As $t\to\infty$ only the terms with $l=0$ survive. 
The steady-state distribution function is thus given by
%
\begin{equation}\label{eq_distr}
 \ups_{eq}(\Omega)=\lim_{t\to\infty} 
\ups(\Omega,t\,|\,\Omega_{0},t_0)
%=Y_{00}^{\ast}(\Omega_0)Y_{00}(\Omega)
=(4\pi)^{-1}.
\end{equation}
%

The conditional and equilibrium distributions, \eqref{sol_pdf} and \eqref{eq_distr}, 
respectively, can be used to coarse-grain the exact 
Langevin equations \eqref{full_langevin} and \eqref{langevin_orientation}. 
The approach taken is to consider the orientation vector $\pp_i(t)$ attached to 
particle $i$ as a stochastic variable and to provide its full statistical 
characterization. 
In spherical coordinates the orientation vector is given explicitly by 
\begin{align} 
\pp(t)&=(p_x(t),p_y(t),p_z(t))^{T}\notag\\
&=(\cos\vp(t)\sin\vt(t),
\sin\vp(t)\sin\vt(t),\cos\vt(t))^ { T },
\end{align}
where $\vp$ and $\vt$ are the azimuthal and polar angles respectively. Using 
\eqref{eq_distr} 
%and $d\Omega=d\vp d\vt\sin\vt$, 
we have that
\begin{align}
 \langle p_z(t)\rangle &=\int d\Omega\, \ups_{eq}(\Omega)\cos\vt=0, 
\end{align}
together with analagous results for the $x$ and $y$ components
\begin{equation}
 \langle p_x(t)\rangle = 0 = \langle p_y(t)\rangle\,.
\end{equation}
Defining the new stochastic variable by $\X_i(t)\equiv v_0\,\pp_i(t)$, its first 
moment is thus given by
\begin{equation}\label{correlations}
 \langle\X_i(t)\rangle=v_0\langle\pp_i(t)\rangle=\boldsymbol{0}.
\end{equation}

Calculation of the equilibrium correlation matrix requires the conditional 
probability distribution function given by \eqref{sol_pdf}. For example, for the $zz$ 
component, we obtain
\begin{align}
\langle p_z(t)p_z(t_0)\rangle &=
\!\!\!\int\!\!d\Omega\!\!\int\!\! 
d\Omega_0 \cos\vt \cos\vt_0 
\ups(\Omega,t\,|\,\Omega_0,t_0)\ups_{eq}(\Omega_0)\notag\\
&=
\frac{1}{3}\int d\Omega \int
d\Omega_0\sum_{l,m}e^{-D_rl(l+1)(t-t_0)}\times\notag\\
&\phantom{=}
\times Y_{10}^{\ast}(\Omega)Y_{lm}(\Omega)Y_{10}
(\Omega_0)Y^{\ast}_{lm}(\Omega_0)\notag\\
&=
\frac{1}{3}e^{-2D_r|t-t_0|},
\end{align}
where we have expressed the cosine functions in terms of spherical harmonics, 
$Y_{10}=\sqrt{3/(4\pi)}\cos\vt=Y_{10}^{\ast}$, and used the orthogonality 
property. 
%In order to avoid to precise which of the 
%time arguments is prior to the other, we simply replaced 
%$(t-t_0)$ by $|t-t_0|$ in the exponential. 
Calculations for the $xx$ and $yy$ 
components are performed in the same spirit. 
% namely by expressing a 
%product of trigonometric functions as a suitable combination of spherical 
%harmonics. 
We thus obtain
\begin{equation}
\langle p_x(t)p_x(t_0)\rangle=\frac{1}{3}e^{-2D_r|t-t_0|}=\langle 
p_y(t)p_y(t_0)\rangle,
\end{equation}
whereas off-diagonal components of the correlation matrix are all zero. 
We can thus conclude that 
\begin{equation}\label{correlations2}
 \langle\X_i(t)\X_j(t')\rangle=v_0^2\langle\pp_i(t)\pp_j(t')\rangle=
 \frac{v_0^2}{3}e^{-2D_r|t-t'|}\boldsymbol{1}\delta_{ij}.
\end{equation}
It has been shown \cite{ghosh2012} that the probability distribution 
function \eqref{sol_pdf} can be well approximated by an expression which 
generalizes the planar Gaussian function to the sphere. 
The new noise function $\X_i(t)$ is thus approximately Gaussian 
distributed with zero mean and exponentially decaying correlations. The 
coarse-grained Langevin equation \eqref{effective_langevin} thus describes 
a stochastic process with additive colored-noise.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%




\section{Approximate Fokker-Planck equation} \label{appendix2}
To derive from \eqref{effective_langevin} an approximate Fokker-Planck equation we 
apply the functional calculus methods of Fox \cite{fox1986}.  
% The 
% coarse-grained Langevin equation \eqref{} describes a non-Markovian process 
% because the stochastic variable $\X$(t) has a finite correlation time $\tau_p$. 
% As a consequence, one can only expect to derive an \emph{approximate} FP 
% equation corresponding to \eqref{}(Langevin colored) which is quite a 
% challenge. Numerous approaches have been proposed [cite Grigolini, vanKampen, 
% Fox, Hanggi. Actually these ref. are all nicely cited in: Venkatesh and 
% Patnaik(PRE 48, 2402, (1993))] in the past decades to handle with 
% this problem. To our knowledge, the most robust one has been proposed by Fox 
% [cite Fox] and is founded on functional calculus. This is the way we will follow 
% here.
We address the one-dimensional case before generalizing to higher dimension. 
Consider the stochastic differential equation
\begin{equation}\label{process}
 \dot{x}(t)= F(x) + g(x)\chi(t)\,,
\end{equation}
where $F(x)$ and $g(x)$ may be nonlinear functions in $x$. If $g(x)=1$, the 
process is then called {\it additive}, otherwise it is called 
{\it multiplicative}. The noise function $\chi(t)$ is by definition Gaussian distributed 
with zero mean. Its second moment determines whether it is a 
{\it white} or {\it colored} noise. As we are interested here in the case of 
additive colored noise we set $g(x)=1$.

In the framework of functional calculus, the Gaussian nature of $\chi(t)$ is expressed by the 
following probability distribution functional,
\begin{equation}\label{pdf}
 P[\chi]=N e^{-\frac{1}{2}\int ds\int ds' \chi(s)\chi(s')K(s-s')}\,,
\end{equation}
where the function $K$ is the inverse of the $\chi$ correlation function and 
the normalization constant is expressed by a path integral over $\chi$,
\begin{equation}\label{N}
 N^{-1}=\int D[\chi] e^{-\frac{1}{2}\int ds\int ds' 
\chi(s)\chi(s')K(s-s')}\,.
\end{equation}
The first and second moments of $\chi$ are given by
\begin{align}
 \langle\chi(t)\rangle &= 0\,,\label{m1}\\
 \langle\chi(t)\chi(s)\rangle&=C(t-s)\,.
\label{m2}
\end{align}
Recalling that the functional derivative may be defined according to
\begin{equation}\label{f_deriv}
 \frac{\delta 
I[\phi]}{\delta\phi(t')}= 
\left.\frac{d}{d\lambda}I[\phi(t)+\lambda\delta(t-t')]\right|_{\lambda=0}\,,
\end{equation}
we now derive two useful identities. The first concerns the functional 
derivative of the probability distribution functional,
\begin{align}
 \frac{\delta P[\chi]}{\delta \chi(t)} &= \frac{\delta N}{\delta 
\chi(t)}e^{-\frac{1}{2}\int ds\int ds'\chi(s)\chi(s')K(s-s')}\notag\\
&\phantom{=} + 
N\frac{\delta}{\delta \chi(t)}e^{-\frac{1}{2}\int ds\int 
ds'\chi(s)\chi(s')K(s-s')}\notag\\
&=
-P[\chi]\int ds K(t-s)\chi(s)\,,
\label{dp/df}
\end{align}
where, using \eqref{f_deriv} and \eqref{m1}, it can be easily shown that $\delta 
N/\delta \chi(t)=0$. The second identity demonstrates the inverse relation between the 
functions $K$ and $C$. 
The second functional derivative of $P[\chi]$ yields
\begin{align}\label{2nd_deriv}
 \frac{\delta^2P[\chi]}{\delta\chi(t')\delta\chi(t)}&=
 P[\chi]\left\{\!\int\!\!ds'\!\!\!\int\!\!ds 
K(t'-s')K(t-s)\chi(s')\chi(s)\right.\notag\\
 &\phantom{=}
 \left.\phantom{\int}\!\!- K(t-t')\right\}\,,
\end{align}
where use of \eqref{dp/df} has been made. Using \eqref{2nd_deriv} and \eqref{m2} 
together with the normalization $\int D[\chi]P[\chi]=1$, leads to
\begin{align}
0 &= \int D[\chi]\frac{\delta^2P[\chi]}{\delta\chi(t')\delta\chi(t)}\notag\\
&=\int\!\!ds'K(t'-s')\int \!\!ds K(t-s)C(s-s')-K(t-t')\,,
\end{align}
which implies that
\begin{equation}\label{delta}
 \int \!\!ds K(t-s)C(s-s')=\delta(t-s')\,. 
\end{equation}

The solution to the stochastic process described by \eqref{process}, namely the 
probability distribution functional for $x(t)$, is given by the formal expression
\begin{equation}\label{sol}
 P(y,t) = \int D[\chi]P[\chi]\delta(y-x(t))\,.
\end{equation}
Taking the time derivative of \eqref{sol} yields
\begin{align}\label{towardsFP}
 \frac{\partial}{\partial t}P(y,t) &= 
-\frac{\partial}{\partial y}[F(y)P(y,t)]\notag\\
&\phantom{=}
-\frac{\partial}{\partial y}\int D[\chi]\delta(y-x(t))P[\chi]\,\chi(t)\,.
\end{align}
The product $P[\chi]\,\chi(t)$ appearing in the second term can be rewritten in 
the following way
\begin{align}\label{PX}
 P[\chi]\,\chi(t) &= P[\chi]\int ds \delta(t-s)\chi(s)\notag\\
 &=
 -\int ds' C(t-s')\frac{\delta P[\chi]}{\delta \chi(s')}\,,
\end{align}
where we have used \eqref{delta} and \eqref{dp/df}. Inserting 
\eqref{PX} back into the second term of \eqref{towardsFP} and integrating by 
parts gives us
\begin{align}\label{45}
 &\int D[\chi]\delta(y-x(t))P[\chi]\chi(t)\notag\\
 &=
 -\int\!\! ds' C(t-s')\!\!\int\!\! D[\chi]\left(\frac{\partial}{\partial 
y}\delta(y-x(t)) \right)\frac{\delta x(t)}{\delta \chi(s')}P[\chi], 
\end{align}
which serves as the exact starting point for Fox's approximation scheme \cite{fox1986}. 

In order to progress further we need to calculate $\delta x(t)/\delta 
\chi(s')$. Applying the functional derivative with respect to $\chi(t')$ on 
\eqref{process} yields a first-order differential equation,
\begin{align}
  \frac{d}{dt}\frac{\delta x(t)}{\delta \chi(t')} &=
  \frac{\delta \dot{x}(t)}{\delta\chi(t')}
  =F'(x)\frac{\delta x(t)}{\delta \chi(t')} +\delta(t-t')\,,
\end{align}
the solution of which is 
\begin{align}\label{dx/dX}
 \frac{\delta x(t)}{\delta \chi(s')} &= \int_{0}^{t}\!\!\!ds\, 
e^{\int_{s}^{t}d\tilde{s}F'(x(\tilde{s}))}\delta(s-s')\notag\\
&=e^{\int_{s'}^{t}ds F'(x(s))}\Theta(t-s')\,,
\end{align}
where $\Theta$ is the Heaviside step function, which we define here as follows
\begin{displaymath}
\Theta(t-s') = \left\{ \begin{array}{ll}
1, & t>s'\\
\frac{1}{2}, & t=s'\\
0, & t<s'\,.
\end{array} \right.
\end{displaymath}
Using \eqref{dx/dX} in \eqref{45} we can rewrite \eqref{towardsFP} in an alternative 
form
\begin{align}\label{46}
 &\frac{\partial}{\partial t}P(y,t) = 
 -\frac{\partial}{\partial y}[F(y)P(y,t)]\notag\\
 &
 +\frac{\partial^2}{\partial 
y^2}\!\!\left\{\!\int_{0}^{t}\!\!\!\!ds'C(t-s')\!\!\!\int\!\!\! 
D[\chi]P[\chi]e^{\int_{s'}^{t}ds F'(x(s))}\delta(y-x(t))\!\!\right\},
\end{align}
which already begins to resemble a Fokker-Planck-type equation. 
However, because of the 
non-Markovian nature of $\int_{s'}^{t}ds F'(x(s))$ appearing in the exponential of 
\eqref{46}, it is apparent that a reduction of this term to an expression 
containing $P(y,t)$ is not possible. An approximation is required. 

The colored noise of interest here is characterized by an 
exponentially decaying correlation function \eqref{correlations}. 
In the literature on non-Markovian processes the time-correlation functions 
are generally notated as follows
\begin{equation}\label{correl}
 C(t-s)=\frac{D}{\tau}e^{-\frac{|t-s|}{\tau}}\,,
\end{equation}
with a diffusion coefficient $D$ and a correlation time $\tau$. 
In order to retain some coherence with the existing literature we will here 
employ the standard 
notation of \eqref{correl} and only use the relation of the parameters in \eqref{correl} 
to those of \eqref{correlations2} at the end of the calculation. 
%Rapid correlation decay implies a small value of $\tau$. 

Returning to \eqref{46}, we 
first perform a change of variable, $t'\equiv t-s'$, in the time-integral,
\begin{equation}\label{nonmarkov}
 \int_{0}^{t}ds'C(t-s')e^{\int_{s'}^{t}ds F'(x(s))}
 =\int_{0}^{t}dt'C(t')e^{\int_{t-t'}^{t}ds F'(x(s))},
\end{equation}
and then expand the time integral over $F'$ in terms of $t'$, 
\begin{equation}\label{approx}
 \int_{t-t'}^{t}ds F'(x(s))\approx F'(x(t))t'-F''(x(t))\dot{x}(t)\frac{t'^2}{2}.
\end{equation}
Neglecting the $t'^2$ term in \eqref{approx} enables the integral in \eqref{nonmarkov} 
to be evaluated
\begin{align}\label{large_t}
 &\int_{0}^{t}ds'C(t-s')e^{\int_{s'}^{t}ds F'(x(s))} \approx 
 \int_{0}^{t}dt'C(t')e^{F'(x(t))t'}\notag\\
 &=\frac{D}{\tau}\int_{0}^{t}dt'e^{-t'(-F'(x(t))+\frac{1}{\tau})}
 \approx \frac{D}{1-\tau F'(x(t))}\,,
\end{align}
where we used \eqref{correl}, and the second approximation results from 
assuming a sufficiently large $t$. We can finally put \eqref{large_t} back into 
\eqref{46} to obtain an approximate Fokker-Planck equation,
\begin{align}\label{approxFP}
 \frac{\partial}{\partial t}P(y,t) &\!=\! 
 -\frac{\partial}{\partial y}[F(y)P(y,t)]+\!D\frac{\partial^2}{\partial 
y^2}\!\!\left(\!\frac{1}{1-\tau F'(y)}P(y,t)\!\right).
\end{align}
This is Fox's result for the approximate Fokker-Planck equation corresponding to 
the non-Markovian process \eqref{process}. Eq.\eqref{approxFP} 
implicity defines a Markovian process, which approximates the non-Markovian 
process of physical interest. However, the question of whether this represents 
the best approximation remains a subject of debate. 
We note that equation \eqref{approxFP} has also been derived by Grigolini {\it et al.} 
\cite{faetti1988} using alternative methods which do not make any assumptions of a short correlation 
time. 

The one-dimensional Fokker-Planck equation \eqref{approxFP} can be generalized without 
much difficulty to describe a three-dimensional system of $N$ particles. 
The dynamics of interest is described by the stochastic 
equation \eqref{effective_langevin}. 
We now adapt the standard notation used above to that employed in the main text, 
namely $P(y,t)\rightarrow \Psi(\boldsymbol{r}^N,t)$, 
$\tau\rightarrow\tau_p=1/(2D_r)$ and $D\rightarrow v_0^2/3$, and recall that 
$D_a=v_0^2/(6D_r)$ and $\zeta^{-1}=\beta D_t$ for the friction coefficient in 
\eqref{effective_langevin}. 
Making the appropriate replacements enables us to write 
the three-dimensional generalization of \eqref{approxFP} 
\begin{align}\label{1st_smol}
 \frac{\partial}{\partial t}&\Psi(\boldsymbol{r}^N,t)=
 -\sum_{i=1}^{N}\nabla_i\cdot
 D_t\left[\beta 
\boldsymbol{F}_i(\boldsymbol{r}^N)-\nabla_{i}\right]\Psi(\boldsymbol{r}^N,
t)\notag\\
 &
-\sum_{i=1}^{N}\nabla_i\cdot \left[-D_a\nabla_{i}\left(\frac{1}{1-\frac{
D_0\nabla_i\cdot\beta\boldsymbol{F}_i(\boldsymbol{r}^N)}{2D_r}}\Psi(\boldsymbol{
r}^N,t) \right)\right].
\end{align}
A simple rearrangement of terms in \eqref{1st_smol} leads directly to  
equations \eqref{smol_eq}-\eqref{force_i}
in the main text.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%

\section{Integral equation theory}\label{integralequation} 
To calculate the steady-state radial distribution function, $g(r)$,  
from the effective pair potential \eqref{potential_pair} we employ an equilibrium 
liquid state integral equation developed by Madden and Rice \cite{madden1980}. 
This soft mean-spherical approximation (SMSA) exploits the Weeks-Chandler-Anderson 
splitting of the pair potential \cite{WCA} into attractive and repulsive contributions,  
$u(r)=u_{\rm rep}(r) + u_{\rm att}(r)$, where the repulsive part is given by
\begin{align}
u_{\rm rep}(r)=
\begin{cases}
u(r)-u(r_{\rm min}) & \;\;\;\;r < r_{\rm min}
\\
0 &\;\;\;\; r > r_{\rm min}
\end{cases}
\end{align}
and the attractive part is given by
\begin{align}
u_{\rm att}(r)=
\begin{cases}
u(r) & \;\;\;\;r > r_{\rm min}
\\
u(r_{\rm min}) &\;\;\;\; r < r_{\rm min}
\end{cases}
\end{align}
where $r_{\rm min}$ is the position of the potential minimum. 
The total correlation function, $h(r)=g(r)-1$, is related to the shorter 
range direct correlation function, 
$c(r)$, by the Ornstein-Zernike equation \cite{hansen_mcdonald1986}
%
\begin{align}
h(r)=c(r) + \rho_b\int d{\rr}' h(|\rr - \rr'|)c(r').
\end{align}
%
The SMSA approximation is given by the closure relation 
%
\begin{align}\label{SMSA}
c(r)=(1 - e^{\beta u_{\rm rep}(r)})g(r) - \beta u_{\rm att}(r). 
\end{align} 
%
For the Lennard-Jones potential the closure relation \eqref{SMSA} has been shown 
to provide results for $g(r)$ which are superior to both Percus-Yevick (PY) and 
Hypernetted Chain (HNC) theories \cite{madden1980}. 
Moreover, the SMSA theory predicts a true spinodal line 
in the parameter space, namely a locus of points for which the static structure 
factor, $S(k)=(1-\rho_b\tilde{c}(k))^{-1}$, diverges at vanishing wavevector. 
This behaviour is a consequence of the assumed asymptotic form of the direct correlation 
function, $c(r)\sim -\beta u_{\rm att}(r)$. 
Other standard integral equation theories, such as PY and HNC, do not exhibit a 
complete spinodal line, but rather a region within which the theory breaks down 
(`no solutions region')
\cite{braderIJTP}.


\section{Brownian Dynamics Simulations} \label{simulations}
In order to benchmark our theoretical predictions we perform Brownian dynamics 
simulations of $N$ particles, randomly initialized without overlap. 
The system is confined to a periodic cubic box, the size of which 
is determined by the number density according to 
$L^3=N/\rho_b$, where $L$ is the side length. 
The Langevin equations of motion, \eqref{full_langevin} and \eqref{langevin_orientation},  
are integrated via a standard Brownian 
dynamics scheme \cite{allen_tildesley} with a constant time step of $\delta t / 
\tau_B = 10^{-5}$.\!
%and a fixed magnitude of the potential of $\varepsilon =1$, where 
%$k_B T = 1$ defines the fundamental unit of energy and $k_B$ gives the 
%Boltzmann constant. 
Both the translational and rotational noise are 
Gaussian random variables with a standard deviation of $\sigma_t\!=\!(2 D_0 T)^{\frac{1}{2}}$ 
and  $\sigma_r\!=\!(2 D_r T)^{\frac{1}{2}}$, respectively. 

For the soft repulsive potential to be considered in this work, $\beta u(r)=r^{-12}$, 
we employ $N\!=\!2000$ particles. 
The potential is truncated and shifted at $r_{\rm cut} / d = 2$. 
To provide good statistics for the static quantities the 
simulations are carried out for $10^6$ time steps, sampling every 1000 
steps, which is equivalent to a total run time of $t_{\rm tot} / \tau_B = 10$ and a 
sampling rate of $\tau_B / t_{sample} = 100$.
For the second system we will consider, the Lennard-Jones system, 
$\beta u(r)\!=\!4\epsilon(r^{-12}\! -r^{-6})$,
we simulate a larger 
system of 5000 particles. 
%randomly 
%initialized without overlap, in a periodic cubic box, which size is determined 
%by the number density according to $L^3=N/\rho$. 
The integration time of the 
equations of motion is the same as in the repulsive system, as is the 
cut-off radius. In this case, the runtime is $10^7$ and the particle positions are sampled 
every $10^4$ steps. 
%The points in the phase diagram Fig. 
%2c are carried out at a fixed density of $\rho = 0.4$ and a fixed $\tau = 
%0.025$, while the interaction parameter $\varepsilon$ and the Perclet number $Pe$ 
%are varied.
 
 


\begin{thebibliography}{10}
  
\bibitem{lyderic2010}
J. Palacci, C. Cottin-Bizonne, C. Ybert, and L. Bocquet   
Phys.Rev.Lett. {\bf 105} 088304 (2010).

\bibitem{erbe2008}
A. Erbe, M. Zientara, L. Baraban, C. Kreidler and P. Leiderer, 
J.Phys.:Condens.Matt. {\bf 20} 404215 (2008).

\bibitem{howse2007}
J. Howse, R.A.L. Jones, A.J. Ryan, T. Gough, R. Vafabakhsh, and R. Golestanian, 
Phys.Rev.Lett. {\bf 99} 048102 (2007).

\bibitem{flagella2005}
R. Dreyfus, J. Baudry, M.L. Roper, M. Fermigier, H.A. Stone and J. Bibette, 
Nature {\bf 437} 862 (2005).

\bibitem{palacci2013}
J. Palacci, S. Sacanna1, A.P. Steinberg, D.J. Pine and P.M. Chaikin, 
Science {\bf 339} 936 (2013).

\bibitem{collective_motion2012} 
T. Vicsek and A. Zafieris,  %Collective motion,
Physics Reports {\bf 517} 71 (2012).

\bibitem{ramaswamy2010} 
S. Ramaswamy,  %The Mechanics and Statistics of Active Matter,
Annual Review of Condensed 
Matter Physics 
 {\bf 1} 323 (2010).

\bibitem{romanczuk2012}
P. Romanczuk, M. B\"ar, W. Ebeling, B. Lindner, and L. Schimansky-Geier, 
Eur.Phys.J. Special Topics {\bf 202} 1 (2012)

\bibitem{cates_bacteria2012}
M.~E. Cates,  %Diffusive transport without detailed balance in motile bacteria,
Rep.Prog.Phys. {\bf 75} 042601 (2012).

\bibitem{cates_tailleur2014}
M.~E. Cates and J. Tailleur, %Motility-Induced Phase Separation,
Annual Review of Condensed Matter Physics, 
{\bf 6} 219 (2015).

\bibitem{tailleur2008}
J. Tailleur and M.~E. Cates, %Statistical Mechanics of Interacting Run-and-Tumble Bacteria, 
Phys.Rev.Lett.  {\bf 100} 218103 (2008). 

\bibitem{fily2012}
Y. Fily and M.~C. Marchetti,  %Athermal Phase Separation of Self-Propelled Particles with No Alignment, 
Phys.Rev.Lett. {\bf 108} 235702 (2012).

\bibitem{stenhammar2013}
J. Stenhammar, A. Tiribocchi, R.J. Allen, D. Marenduzzo 
and M.E. Cates, 
Phys.Rev.Lett.  {\bf 111} 145702 (2013).

\bibitem{redner2013} 
G.~S. Redner, A. Baskaran and M.~F. Hagan,  %Reentrant phase behavior in active colloids with attraction, 
Phys.Rev.E {\bf 88} 012305 (2013).

\bibitem{levis_berthier2014}
D. Levis and L. Berthier,  %Clustering and heterogeneous dynamics in a kinetic Monte Carlo 
%model of self-propelled hard disks,
Phys.Rev.E  {\bf 89} 062301 (2014).

\bibitem{schwarz-linek2012}
J. Schwarz-Linek {\it et al.}, %Phase separation and rotor self-assembly in active particle suspensions,
Proceedings of the National Academy of Sciences (USA) {\bf 109} 4052 (2012).

\bibitem{bocquetPRX_2014}
F. Ginot, I. Theurkauff, D. Levis, C. Ybert, L. Bocquet, L. Berthier 
and C. Cottin-Bizonne, 
Phys.Rev.X {\bf 5}, 011004 (2014). 

\bibitem{wysocki2014}
A. Wysocki and R.G. Winkler and G. Gompper,
%Cooperative Motion of Active Brownian Spheres in Three-Dimensional Dense Suspensions, 
EPL {\bf 105} 48004 (2014). 

\bibitem{cates2013}
M.~E. Cates and J. Tailleur,  %Motility-Induced Phase Separation, 
EPL {\bf 101} 20010 (2013).

\bibitem{bialke2013}
J. Bialk\'e, H. L\"owen and T. Speck, 
EPL {\bf 103} 30008 (2013).

\bibitem{cates_kardar}
A.~P. Solon, Y. Fily, A. Baskaran, M. E. Cates, Y. Kafri, M. Kardar and 
J. Tailleur, 
arXiv:1412.3952 (2014).

\bibitem{takatori}
S.C. Takatori, W. Yan, and J.F. Brady, 
Phys.Rev.Lett. {\bf 113} 028103 (2014).

\bibitem{stenhammar2014}
J. Stenhammar, D. Marenduzzo, R.J. Allen and M.E. Cates, 
Soft Matter {\bf 10} 1489 (2014).

\bibitem{vankampen1976}
N.~G. van Kampen,  %Stochastic differential equations,
Physics Reports {\bf 24}: 171 (1976).

\bibitem{vankampen1998}
N.~G. van Kampen,  %Remarks on Non-Markov Processes,
Brazilian Journal of Physics {\bf 28} 90 (1998).

\bibitem{faetti1988}
S. Faetti, L. Fronzoni, P. Grigolini and R. Mannella, 
J.Stat.Phys. {\bf 52} 951 (1988).

\bibitem{fox1986}
R.~F. Fox,  %Functional-calculus approach to stochastic differential equations,
Phys.Rev.A {\bf 33} 467 (1986).

\bibitem{fox_weak1986}
R.~F. Fox,  %Uniform convergence to an effective Fokker-Planck equation for weakly colored noise,
Phys.Rev.A {\bf 34} 4525(R) (1986).

\bibitem{pototsky2012}
A. Pototsky and H. Stark,  %Active Brownian particles in two-dimensional traps,
EPL {\bf 98} 50004 (2012).

\bibitem{madden1980}
W. Madden and S. Rice, % The mean spherical approximation and effective pair potentials in liquids,
J.Chem.Phys. {\bf 72} 4208 (1980).

\bibitem{barker_henderson1976}
Useful tables of simulation data are to be found in: 
J. Barker and D.~H. Henderson,
Reviews of Modern Physics {\bf 48} 587 (1976).

\bibitem{sear1997}
R. Sear and W. Gelbart, %Microphase separation versus the vapour-liquid transition in 
%systems of spherical particles
J.Chem.Phys. {\bf 110}: 4582 (1997).

\bibitem{archer2007}
A.~J. Archer and N. Wilding, %Phase behaviour of a fluid with competing attractive 
%and repulsive interactions
Phys.Rev.E {\bf 76}: 031501 (2007). 

\bibitem{brader_triplet}
J.M. Brader, 
J.Chem.Phys. {\bf 128} 104503 (2008).

\bibitem{maggi2015}
C. Maggi, U.M.B. Marconi, N. Gnan and R. Di Leonardo,
arXiv:1503.03123v2

\bibitem{faragebrader_arxiv}
T.F.F. Farage and J.M. Brader, 
Arxiv: 1403.0928 

\bibitem{ni2013}
R. Ni, M.A. Cohen Stuart and M. Dijkstra,  
%Pushing the glass transition towards random
%lose packing using self-propelled hard spheres, 
Nature Communications {\bf 4} 2704 (2013).

\bibitem{berthier2013}
L. Berthier and J. Kurchan,  
%Non-equilibrium glass transitions in driven and active matter, 
Nature Physics {\bf 9} 310 (2013).

\bibitem{szamel_berthier2015}
G. Szamel, E. Flenner and L. Berthier, 
Arxiv: 1501.01333

\bibitem{gardiner}
C. Gardiner, {\it Handbook of stochastic methods} (Springer, Berlin, 1985).

\bibitem{debye1929}
P. Debye {\it Polar Molecules}, (The Chemical Catalog Com-
pany, Inc., 1929).

\bibitem{fatuzzo1967}
E. Fatuzzo and P.R. Mason,  %A theory of dielectric relaxation in polar liquids 
Proc. Phys. Soc. {\bf 90} 741 (1967).

\bibitem{zwanzig1970}
T.-W. Nee and R. Zwanzig,  
%Theory of Dielectric Relaxation in Polar Liquids  
J. Chem. Phys. {\bf 52} 6353 (1970).

\bibitem{berne1975}
B.J. Berne,  %A self-consistent theory of rotational diffusion
J. Chem. Phys. {\bf 62} 1154 (1975).

\bibitem{arfken2013}
G.B. Arfken, H.J. Weber, and F.A. Harris, 
{\it Mathematical Methods for Physicists} (Elsevier, 2013).

\bibitem{ghosh2012}
A. Ghosh, J. Samuel, and S. Sinha,  %A Gaussian for diffusion on the sphere
EPL {\bf 98} 30003 (2012).

\bibitem{WCA}
J. Weeks and D. Chandler and H. Anderson, % Role of Repulsive Forces in Determining 
%the Equilibrium Structure of Simple Liquids,
J.Chem.Phys. {\bf 54} 5237 (1971).

\bibitem{hansen_mcdonald1986}
J.-P. Hansen and I.~R. McDonald, 
{\it Theory of simple liquids} (Academic press, London, 1986).

\bibitem{braderIJTP}
J.M. Brader, %Solution of the Ornstein-Zernike equation in the critical region, 
International Journal of Thermophysics {\bf 27} 394 (2006).

\bibitem{allen_tildesley}
M.~P. Allen and D.~J. Tildesley (1991) 
{\it Computer Simulation of Liquids} 
(Oxford University Press, 1991).




%\bibitem{enculescu2011}
%M. Enculescu and H. Stark,  %Active Colloidal Suspensions Exhibit Polar Order under Gravity,
%hys.Rev.Lett. {\bf 107} 058301 (2011).

%\bibitem{angelani2011}
%L. Angelani and A. Costanzo and R. Di Leonardo,  %Active ratchets
%EPL {\bf 96} 68002 (2011).

%\bibitem{ai2013}
%B. Ai, Q. Chen, Y. He, F. Li, and W. Zhong, 
%Phys. Rev. E {\bf 88} 062129 (2013).

 

%\bibitem{archer_evans_spinodal2004}
%A.~J. Archer and R. Evans (2004) Dynamical density functional theory and its 
%application to spinodal decomposition,
%J. Chem. Phys. {\bf 121}: 4246

%\bibitem{reinhardt2012}
%J. Reinhardt and J.~M. Brader (2012) Dynamics of localized particles from density functional theory,
%Phys.Rev.E {\bf 85}: 011404

%\bibitem{speck2014}
%T. Speck {\it et al.} (2014) Effective Cahn-Hilliard Equation for the Phase Separation of 
%Active Brownian Particles, 
%Phys. Rev Lett- {\bf 112}: 218304 

%\bibitem{hansen_mcdonald1986}
%J.-P. Hansen and I.~R. McDonald, 
%{\it  Theory of simple liquids} (Academic press, London, 1986).

%\bibitem{WCA}
%J. Weeks and D. Chandler and H. Anderson (1971) Role of Repulsive Forces in Determining 
%the Equilibrium Structure of Simple Liquids,
%J.Chem.Phys. {\bf 54}: 5237


%\bibitem{rosenfeld89}
%Y. Rosenfeld (1989) Free-energy model for the inhomogeneous hard-sphere fluid mixture 
%and density-functional theory of freezing,
%Phys.Rev.Lett. {\bf 63}: 980

%\bibitem{evans_tarazona1984}
%R.~Evans and P.~Tarazona (1984) Theory of condensation in narrow capillaries, 
%Phys. Rev. Lett. {\bf 52}: 557

%\bibitem{bob92}
%R. Evans (1992) Ch.3 in {\it Fundamentals of inhomogeneous fluids} 
%(Marcel Dekker, 1992). 

%\bibitem{schmidt_capillary2004}
%D.~G.~A.~L. Aarts, M.~Schmidt, H.~N.~W. Lekkerkerker (2004) 
%Direct Visual Observation of Thermal Capillary Waves, 
%Science {\bf 304}: 847

%\bibitem{speck_interface2014}
%J.~Bialk\'e, H. L\"owen and T.~Speck (2014) 
%Negative interfacial tension in phase-separated active suspensions, 
%ArXiv :1412.4601 






\end{thebibliography}




\end{document}































1504.00192.extracted_bib
{lyderic2010 J. Palacci, C. Cottin-Bizonne, C. Ybert, and L. Bocquet Phys.Rev.Lett. {105 088304 (2010).
1504.00192.extracted_bib
{erbe2008 A. Erbe, M. Zientara, L. Baraban, C. Kreidler and P. Leiderer, J.Phys.:Condens.Matt. {20 404215 (2008).
1504.00192.extracted_bib
{howse2007 J. Howse, R.A.L. Jones, A.J. Ryan, T. Gough, R. Vafabakhsh, and R. Golestanian, Phys.Rev.Lett. {99 048102 (2007).
1504.00192.extracted_bib
{flagella2005 R. Dreyfus, J. Baudry, M.L. Roper, M. Fermigier, H.A. Stone and J. Bibette, Nature {437 862 (2005).
1504.00192.extracted_bib
{palacci2013 J. Palacci, S. Sacanna1, A.P. Steinberg, D.J. Pine and P.M. Chaikin, Science {339 936 (2013).
1504.00192.extracted_bib
{collective_motion2012 T. Vicsek and A. Zafieris, %Collective motion, Physics Reports {517 71 (2012).
1504.00192.extracted_bib
{ramaswamy2010 S. Ramaswamy, %The Mechanics and Statistics of Active Matter, Annual Review of Condensed Matter Physics {1 323 (2010).
1504.00192.extracted_bib
{romanczuk2012 P. Romanczuk, M. B\"ar, W. Ebeling, B. Lindner, and L. Schimansky-Geier, Eur.Phys.J. Special Topics {202 1 (2012)
1504.00192.extracted_bib
{cates_bacteria2012 M.~E. Cates, %Diffusive transport without detailed balance in motile bacteria, Rep.Prog.Phys. {75 042601 (2012).
1504.00192.extracted_bib
{cates_tailleur2014 M.~E. Cates and J. Tailleur, %Motility-Induced Phase Separation, Annual Review of Condensed Matter Physics, {6 219 (2015).
1504.00192.extracted_bib
{tailleur2008 J. Tailleur and M.~E. Cates, %Statistical Mechanics of Interacting Run-and-Tumble Bacteria, Phys.Rev.Lett. {100 218103 (2008).
1504.00192.extracted_bib
{fily2012 Y. Fily and M.~C. Marchetti, %Athermal Phase Separation of Self-Propelled Particles with No Alignment, Phys.Rev.Lett. {108 235702 (2012).
1504.00192.extracted_bib
{stenhammar2013 J. Stenhammar, A. Tiribocchi, R.J. Allen, D. Marenduzzo and M.E. Cates, Phys.Rev.Lett. {111 145702 (2013).
1504.00192.extracted_bib
{redner2013 G.~S. Redner, A. Baskaran and M.~F. Hagan, %Reentrant phase behavior in active colloids with attraction, Phys.Rev.E {88 012305 (2013).
1504.00192.extracted_bib
{levis_berthier2014 D. Levis and L. Berthier, %Clustering and heterogeneous dynamics in a kinetic Monte Carlo %model of self-propelled hard disks, Phys.Rev.E {89 062301 (2014).
1504.00192.extracted_bib
{schwarz-linek2012 J. Schwarz-Linek {et al., %Phase separation and rotor self-assembly in active particle suspensions, Proceedings of the National Academy of Sciences (USA) {109 4052 (2012).
1504.00192.extracted_bib
{bocquetPRX_2014 F. Ginot, I. Theurkauff, D. Levis, C. Ybert, L. Bocquet, L. Berthier and C. Cottin-Bizonne, Phys.Rev.X {5, 011004 (2014).
1504.00192.extracted_bib
{wysocki2014 A. Wysocki and R.G. Winkler and G. Gompper, %Cooperative Motion of Active Brownian Spheres in Three-Dimensional Dense Suspensions, EPL {105 48004 (2014).
1504.00192.extracted_bib
{cates2013 M.~E. Cates and J. Tailleur, %Motility-Induced Phase Separation, EPL {101 20010 (2013).
1504.00192.extracted_bib
{bialke2013 J. Bialk\'e, H. L\"owen and T. Speck, EPL {103 30008 (2013).
1504.00192.extracted_bib
{cates_kardar A.~P. Solon, Y. Fily, A. Baskaran, M. E. Cates, Y. Kafri, M. Kardar and J. Tailleur, arXiv:1412.3952 (2014).
1504.00192.extracted_bib
{takatori S.C. Takatori, W. Yan, and J.F. Brady, Phys.Rev.Lett. {113 028103 (2014).
1504.00192.extracted_bib
{stenhammar2014 J. Stenhammar, D. Marenduzzo, R.J. Allen and M.E. Cates, Soft Matter {10 1489 (2014).
1504.00192.extracted_bib
{vankampen1976 N.~G. van Kampen, %Stochastic differential equations, Physics Reports {24: 171 (1976).
1504.00192.extracted_bib
{vankampen1998 N.~G. van Kampen, %Remarks on Non-Markov Processes, Brazilian Journal of Physics {28 90 (1998).
1504.00192.extracted_bib
{faetti1988 S. Faetti, L. Fronzoni, P. Grigolini and R. Mannella, J.Stat.Phys. {52 951 (1988).
1504.00192.extracted_bib
{fox1986 R.~F. Fox, %Functional-calculus approach to stochastic differential equations, Phys.Rev.A {33 467 (1986).
1504.00192.extracted_bib
{fox_weak1986 R.~F. Fox, %Uniform convergence to an effective Fokker-Planck equation for weakly colored noise, Phys.Rev.A {34 4525(R) (1986).
1504.00192.extracted_bib
{pototsky2012 A. Pototsky and H. Stark, %Active Brownian particles in two-dimensional traps, EPL {98 50004 (2012).
1504.00192.extracted_bib
{madden1980 W. Madden and S. Rice, % The mean spherical approximation and effective pair potentials in liquids, J.Chem.Phys. {72 4208 (1980).
1504.00192.extracted_bib
{barker_henderson1976 Useful tables of simulation data are to be found in: J. Barker and D.~H. Henderson, Reviews of Modern Physics {48 587 (1976).
1504.00192.extracted_bib
{sear1997 R. Sear and W. Gelbart, %Microphase separation versus the vapour-liquid transition in %systems of spherical particles J.Chem.Phys. {110: 4582 (1997).
1504.00192.extracted_bib
{archer2007 A.~J. Archer and N. Wilding, %Phase behaviour of a fluid with competing attractive %and repulsive interactions Phys.Rev.E {76: 031501 (2007).
1504.00192.extracted_bib
{brader_triplet J.M. Brader, J.Chem.Phys. {128 104503 (2008).
1504.00192.extracted_bib
{maggi2015 C. Maggi, U.M.B. Marconi, N. Gnan and R. Di Leonardo, arXiv:1503.03123v2
1504.00192.extracted_bib
{faragebrader_arxiv T.F.F. Farage and J.M. Brader, Arxiv: 1403.0928
1504.00192.extracted_bib
{ni2013 R. Ni, M.A. Cohen Stuart and M. Dijkstra, %Pushing the glass transition towards random %lose packing using self-propelled hard spheres, Nature Communications {4 2704 (2013).
1504.00192.extracted_bib
{berthier2013 L. Berthier and J. Kurchan, %Non-equilibrium glass transitions in driven and active matter, Nature Physics {9 310 (2013).
1504.00192.extracted_bib
{szamel_berthier2015 G. Szamel, E. Flenner and L. Berthier, Arxiv: 1501.01333
1504.00192.extracted_bib
{gardiner C. Gardiner, {Handbook of stochastic methods (Springer, Berlin, 1985).
1504.00192.extracted_bib
{debye1929 P. Debye {Polar Molecules, (The Chemical Catalog Com- pany, Inc., 1929).
1504.00192.extracted_bib
{fatuzzo1967 E. Fatuzzo and P.R. Mason, %A theory of dielectric relaxation in polar liquids Proc. Phys. Soc. {90 741 (1967).
1504.00192.extracted_bib
{zwanzig1970 T.-W. Nee and R. Zwanzig, %Theory of Dielectric Relaxation in Polar Liquids J. Chem. Phys. {52 6353 (1970).
1504.00192.extracted_bib
{berne1975 B.J. Berne, %A self-consistent theory of rotational diffusion J. Chem. Phys. {62 1154 (1975).
1504.00192.extracted_bib
{arfken2013 G.B. Arfken, H.J. Weber, and F.A. Harris, {Mathematical Methods for Physicists (Elsevier, 2013).
1504.00192.extracted_bib
{ghosh2012 A. Ghosh, J. Samuel, and S. Sinha, %A Gaussian for diffusion on the sphere EPL {98 30003 (2012).
1504.00192.extracted_bib
{WCA J. Weeks and D. Chandler and H. Anderson, % Role of Repulsive Forces in Determining %the Equilibrium Structure of Simple Liquids, J.Chem.Phys. {54 5237 (1971).
1504.00192.extracted_bib
{hansen_mcdonald1986 J.-P. Hansen and I.~R. McDonald, {Theory of simple liquids (Academic press, London, 1986).
1504.00192.extracted_bib
{braderIJTP J.M. Brader, %Solution of the Ornstein-Zernike equation in the critical region, International Journal of Thermophysics {27 394 (2006).
1504.00192.extracted_bib
{allen_tildesley M.~P. Allen and D.~J. Tildesley (1991) {Computer Simulation of Liquids (Oxford University Press, 1991). %
1504.00192.extracted_bib
{enculescu2011 %M. Enculescu and H. Stark, %Active Colloidal Suspensions Exhibit Polar Order under Gravity, %hys.Rev.Lett. {107 058301 (2011). %
1504.00192.extracted_bib
{angelani2011 %L. Angelani and A. Costanzo and R. Di Leonardo, %Active ratchets %EPL {96 68002 (2011). %
1504.00192.extracted_bib
{ai2013 %B. Ai, Q. Chen, Y. He, F. Li, and W. Zhong, %Phys. Rev. E {88 062129 (2013). %
1504.00192.extracted_bib
{archer_evans_spinodal2004 %A.~J. Archer and R. Evans (2004) Dynamical density functional theory and its %application to spinodal decomposition, %J. Chem. Phys. {121: 4246 %
1504.00192.extracted_bib
{reinhardt2012 %J. Reinhardt and J.~M. Brader (2012) Dynamics of localized particles from density functional theory, %Phys.Rev.E {85: 011404 %
1504.00192.extracted_bib
{speck2014 %T. Speck {et al. (2014) Effective Cahn-Hilliard Equation for the Phase Separation of %Active Brownian Particles, %Phys. Rev Lett- {112: 218304 %
1504.00192.extracted_bib
{hansen_mcdonald1986 %J.-P. Hansen and I.~R. McDonald, %{Theory of simple liquids (Academic press, London, 1986). %
1504.00192.extracted_bib
{WCA %J. Weeks and D. Chandler and H. Anderson (1971) Role of Repulsive Forces in Determining %the Equilibrium Structure of Simple Liquids, %J.Chem.Phys. {54: 5237 %
1504.00192.extracted_bib
{rosenfeld89 %Y. Rosenfeld (1989) Free-energy model for the inhomogeneous hard-sphere fluid mixture %and density-functional theory of freezing, %Phys.Rev.Lett. {63: 980 %
1504.00192.extracted_bib
{evans_tarazona1984 %R.~Evans and P.~Tarazona (1984) Theory of condensation in narrow capillaries, %Phys. Rev. Lett. {52: 557 %
1504.00192.extracted_bib
{bob92 %R. Evans (1992) Ch.3 in {Fundamentals of inhomogeneous fluids %(Marcel Dekker, 1992). %
1504.00192.extracted_bib
{schmidt_capillary2004 %D.~G.~A.~L. Aarts, M.~Schmidt, H.~N.~W. Lekkerkerker (2004) %Direct Visual Observation of Thermal Capillary Waves, %Science {304: 847 %
1504.00192.extracted_bib
{speck_interface2014 %J.~Bialk\'e, H. L\"owen and T.~Speck (2014) %Negative interfacial tension in phase-separated active suspensions, %ArXiv :1412.4601