PHYSICAL REVIEW E 103, 012605 (2021)

Pair correlation of dilute active Brownian particles: From low-activity dipolar correction tohigh-activity algebraic depletion wings
Alexis Poncet ,1,2 Olivier Bénichou,1 Vincent Démery ,3,4,* and Daiki Nishiguchi

2

1

LPTMC, CNRS/Sorbonne Université, 4 Place Jussieu, F-75005 Paris, France
Department of Physics, The University of Tokyo, Hongo 7-3-1, Tokyo, 113-0033, Japan
3
Gulliver, CNRS, ESPCI Paris, PSL Research University, 10 rue Vauquelin, Paris, France
4
Univ Lyon, ENS de Lyon, Univ Claude Bernard Lyon 1, CNRS, Laboratoire de Physique, F-69342 Lyon, France
2(Received 17 June 2020; revised 12 October 2020; accepted 16 December 2020; published 13 January 2021)
We study the pair correlation of active Brownian particles at low density using numerical simulations andanalytical calculations. We observe a winged pair correlation: While particles accumulate in front of an active
particle as expected, the depletion wake consists of two depletion wings. In the limit of soft particles, we obtaina closed equation for the pair correlation, allowing us to characterize the depletion wings. In particular, we
unveil two regimes at high activity, where the wings adopt a self-similar profile and decay algebraically. We also
perform experiments of self-propelled Janus particles and indeed observe the depletion wings.
DOI: 10.1103/PhysRevE.103.012605

I. INTRODUCTION

The pair correlation function has played a pivotal role in
our understanding of the structure of equilibrium liquids [1];the same is to be expected for active liquids, such as bacterialcolonies [2,3], flocks of starlings [4,5] or colloidal rollers[6–8], or assemblies of Janus particles [9–12]. It has indeed
been used to quantify the order in bacterial colonies [2] and toinfer the interactions in bird flocks [4,5]. However, its broad
utilization is precluded by the lack of analytical results on
its general form, notably little is known about its regions ofpositive and negative sign, even in the homogeneous phaseand in the absence of alignment interactions [13–15].
The theoretical prediction of the pair correlation is indeed
an important challenge, which has been undertaken for minimal models without alignment, such as active Brownian or
Ornstein-Uhlenbeck particles (ABP or AOUP). So far, evenfor these minimal models, the study of the pair correlation
focused mostly on the explanation of the phase separation that
occurs at high activity [16–18]. Two main approaches havebeen followed. The first consists in computing the angularaverage of the pair correlation due to an activity represented
by a persistent translational noise. An attractive term has beenfound and interpreted as a higher tendency to phase separate[15,19,20]. However, this approach does not, by nature, retain
the angular dependence of the correlations. The second approach consists of a quantitative prediction of the polar paircorrelation using a closure of the many body Smoluchowski
equation leading to the effective velocity, whose decrease withdensity may explain the phase separation [13,14,21–23]. Thismethod requires a numerical solution of nonlinear equations
and does not provide explicit predictions. Finally, up to now
there is no analytical characterization of the pair correlation.

*vincent.demery@espci.psl.eu

2470-0045/2021/103(1)/012605(13)In this article, we address the global shape of the paircorrelation of ABP in the dilute and homogeneous regime[17]. We first observe a winged pair correlation in numericalsimulations of ABP: While particles accumulate in front of anactive particle as expected, the depletion wake consists of two
depletion wings (Fig. 1). We then resort to a linearized Dean
equation to obtain a closed equation for the pair correlation.
Solving this equation under different limits, we unveil threedifferent regimes for the correlations, which we organize on aphase diagram. In two regimes, the wings adopt a self-similar
profile and decay algebraically. Last, we measure the paircorrelation in experiments of self-propelled Janus particles
and qualitatively observe the depletion wings.
II. MODEL AND NUMERICAL SIMULATIONSFirst, we perform numerical simulations of ABP in dimension two interacting via a pairwise harmonic potential
V (r) = 2 (1 − r/a)2 for r < a, where  sets the interaction
strength and a is the particle diameter [Fig. 1(a), see Appendix
A1 for details]. The position ri (t ) and orientation θi (t ) of theparticle i followṙi = −γ −1 ∇iθ̇i =
2Dr νi ,V (ri − r j ) + U êθi +√2Dηi ,

(1)

j=i

(2)

where U is the propulsion velocity, D is the translational
diffusion coefficient, Dr is the rotational diffusion coefficient,êθ is the unit vector with polar angle θ , and ηi and νi arenormalized Gaussian white noises. The diameter a of the particles and the friction coefficient γ can be set to one through a
rescaling; the parameters , U , D, and Dr are given in arbitrary
units and only their relative values are important. We use a
low area fraction φ  0.04 to focus on two-body effects and012605-1©2021 American Physical Society

ALEXIS PONCET et al.

PHYSICAL REVIEW E 103, 012605 (2021)

FIG. 2. Phase diagram showing the different regimes as a function of Péclet number and observation length. Blue arrows indicate
the regimes explored in the simulations [Figs. 1(b)–1(e)] and in theexperiments (Fig. 4), and the observation length ranges from thediameter a of the particles to 10a. Red arrows point to the limiting
regimes explored in Fig. 3.

FIG. 1. (a) Snapshot of numerical simulations of ABP. The
red arrows indicate the orientations of the particles. The blue
axes represent the local frame used to compute the pair correlations. (b–e) Evolution of the correlations with decreasing Dr . φ =0.04 and  : D : U = 1 : 0.1 : 10, and Dr = 10, 1, 0.1, 0.01 (Pe =10, 31.6, 100, 316). Top: simulations. Bottom: numerical integrationof Eq. (7).

remain in the homogeneous phase [16–18] (see discussion in
Appendix A3).
The polar pair correlation C(r, θ , θ  ) is defined from thedensity fieldfor the positions ri (t ) and orientations θi (t ),f (r, θ , t ) = i δ(ri (t ) − r)δ(θi (t ) − θ ):C(r, θ , θ  ) = f (0, θ ) f (r, θ  ) δ(r)δ(θ − θ  )− 1,−[ρ/(2π )]2ρ/(2π )

(3)

where ρ = 4φ/(π a2 ) is the density; the correlation of a particle with itself is removed in the second term, and the r →∞ limit is removed in the third. From rotational invariance,C(r, 0, θ  ) contains all the information in C(r, θ , θ  ). We focus
on the density of particles in the reference frame of a given
particle (blue axes in Fig. 1(a) [2]), which retains the polar
character of the correlations:
 2π
1
B(r) =
C(r, 0, θ  )dθ  .
(4)
2π 0Note that if translational and rotational degrees of freedom
decouple, for instance, in absence of propulsion, the isotropic
pair correlation is recovered: C(r, θ , θ  ) = B(r) = h(r) [1].
The polar pair correlation B(r) in the simulations for : D : U = 1 : 0.1 : 10, and Dr = 10, 1, 0.1, 0.01 is shownin Figs. 1(b)–1(e). At large rotational diffusion, Dr = 10,
the correlation is nearly axisymmetric and decays quickly[Fig. 1(b)]. At smaller rotational diffusion, Dr  1, as ex-pected, the correlation is positive in front of the particle,
indicating an accumulation of other particles. However, while
a depletion wake is expected behind the active particle, asin active microrheology [24,25] or in driven binary mixtures[26], the depletion concentrates in two “wings” on the sidesof the particle. As Dr decreases, the characteristic depletion
wings appear, their length increases and their curvature decreases.
III. THEORY AND LIMITING REGIMES

To rationalize these different regimes, we consider thelengthscales and dimensionless parameters of the problem.
In the dilute limit considered here, the structure of the paircorrelation beyond the size of a particle is controlled by theparameters U , D, and Dr , which combine into three lengthscales,

D
D
U
,
,
,
(5)
r =
U =
p =
U
Dr
Dr
whose relative values are set by the Péclet number,
Pe = √

U
=
DDr

U
r

=

p

.

(6)

U

The Péclet number varies between Pe = 10 and 316 in
Figs. 1(b)–1(e). For spherical particles the coefficients D and
Dr are related through Dr ∼ D/a2 ; here we regard them as
independent parameters to disentangle the effects of translational and rotational diffusion. The form of the pair correlation
thus depends on the Péclet number and the relative valueof the observation length r with the three lengthscales, so thatthe numerical simulations can be placed on a parameter plane(blue arrows in Fig. 2).
To account for the observed correlations and characterize
the shape of the depletion wings, we resort to a linearized

012605-2

PAIR CORRELATION OF DILUTE ACTIVE BROWNIAN …

PHYSICAL REVIEW E 103, 012605 (2021)

FIG. 3. Limiting regimes. (a)–(c) Small propulsion (U → 0). (a) Numerical correlation (φ = 0.04 and  : D : Dr : U = 1 : 0.1 : 0.5 :0.05). (b) Cut of the previous correlation at radius r0 = 1.5 (dashed black line), and theoretical prediction (gray line). (c) Rescaled plotof the first Fourier coefficient B1 (r) for φ = 0.04 and  : D : U = 1 : 0.1 : 0.05 for Dr = 0.1, 0.2, 0.5, 1 (blue to red), the gray line is thetheoretical prediction at large distance. (d), (e) No rotational diffusion (Dr = 0). (d) Numerical correlation (φ = 0.04 and  : D : Dr : U = 1 :
0.1 : 0 : 10). (e) Rescaled horizontal cuts of the correlation from (d), the gray line is the theoretical prediction (8). (f), (g) No translational
diffusion (D = 0). (f) Numerical correlation (φ = 0.04 and  : D : Dr : U = 1 : 0 : 0.1 : 10). (g) Rescaled horizontal cuts of the correlationfrom (f) according to Eq. (10).

Dean equation for the density field f (r, θ , t ) [25–30], which
is valid for weak interactions. At low density, we show in
Appendix B1 that the pair correlation satisfies2
2D∇ 2 + Dr ∂θ2 + ∂θ2 + U (êθ − êθ  ) · ∇ C = − ∇ 2V.
γ(7)Solving this equation numerically (Appendix C), we obtain an excellent agreement with the numerical simulations[Figs. 1(b)–1(e)]. We now show that this equation captures
the structure of the pair correlation by examining its limiting
regimes U → 0, Dr → 0, and D → 0.
For small propulsion velocity U , which corresponds tothe left side of the phase diagram, Eq. (7) can be solved
perturbatively as shown in Appendix B2. To order one in U ,we get a dipolar correction to the equilibrium radial correlation heq (r): B(r, θ ) = heq (r) + U cos(θ )B1 (r), where the first
Fourier coefficient B1 (r) decays exponentially over a length
U . This prediction is compared to the numerical simulationsin Figs. 3(a)–3(c); a quantitative agreement with the prediction is obtained, without any adjustable parameter. Note thatat order U 2 , an attractive correction to the equilibrium paircorrelation heq (r) is found, which is compatible with the results obtained when the activity is introduced in the form of acolored noise [15,19,20].
We turn to the high-activity regime, and focus first on
the limit Dr = 0, where the wings are entirely deployed[Fig. 3(d)]; this limit corresponds to the bottom right cornerof the phase diagram (Fig. 2). Equation (7) can be solved in
Fourier space (Appendix B3), yieldingk 2Ṽ (k) 2πdθ B̃(k) = −.
(8)
2πU 0 2 r k 2 + ik · (ê0 − êθ  )The integrand has the form of the correlations in a driven
binary mixture [26], and the integration over θ  mirrors that
here there are particles moving in all directions. At distances
below r , Eq. (8) yields a dipolar correlation, as in the lowactivity regime. Beyond r , the form of the correlations canbe obtained by performing the integral and focusing on the
singularity at the origin in Fourier space:
B(x, y) rṼ (0) 1
F
D0 y2
r

x
1/3
2/3
r |y|

,

(9)

where the explicit form of the scaling function F (u) is given
in Appendix B3. The scaling Eq. (9) can be tested by plotting rescaled cuts of the correlation at different values ofy [Fig. 3(e)]; as predicted, the cuts collapse to the scaling
function F (u).
We now address the shape of the wings in presence of afinite rotational diffusion Dr . We take the limit D = 0, whichcorresponds to the top right corner of the phase diagram(Fig. 2), where the only available lengthscale is the persistence length p (Appendix B4). Beyond p , particles loose
the memory of their orientation and their correlation vanish.
Below p , the solution of Eq. (7) is dominated by angles closeto (θ , θ  )  (0, 0) and (π , π ), and we obtain the following
scaling form:
B(x, y) r

1/3Ṽ (0) 4p
p x
G
;
4
Dr y
|y|4/3p(10)

the function G can be obtained numerically. This prediction
is in very good agreement with the numerical simulations[Figs. 3(f) and 3(g)]. Last, the transition with the scaling form
obtained in the limit Dr = 0 can be obtained by matching012605-3

ALEXIS PONCET et al.

PHYSICAL REVIEW E 103, 012605 (2021)

2/3
the widths of the two profiles, x ∼ 1/3[Eq. (9)] andry
4/3x ∼ −1/3y[Eq.
(10)],
leading
to
y
∼
=
r p
U.
p
The structure of the correlations is now characterized(Fig. 2). For Pe < 1, the activity generates a dipolar correction
to the equilibrium correlations for lengths r < U . For Pe > 1,
the dipolar correction crosses-over to depletion wings at ascale r ; the shape of the wings is given by the scaling forms
Eq. (9) for r < r < U and Eq. (10) for U < r < p ; thewings decay exponentially beyond p . The behavior of thecorrelations is summarized on the phase diagram (Fig. 2).
For weak interactions, we have obtained a quantitative
agreement between our predictions derived from Eq. (7) andthe numerical simulations. Moreover, the left-hand side ofEq. (7), which controls the structure of the correlations, is alsopresent in the two-body Smoluchowski equation [22], which
is valid in the dilute limit for any interaction strength. We conclude that the structure of the correlations predicted here holds
in the dilute limit for any interaction strength, which we checkin Appendix A2 with numerical simulations of hard particles.
We finally note that, in practice, D and Dr are not taken as
independent parameters. For spherical particles, where D ∼a2 Dr , the Péclet number defined here takes the more common
forms Pe = aU/D = U/(aDr ), and U ∼ a. At small Péclet,the dipolar correction, which decays over a length U , should
thus be barely observable. At large Péclet, depletion wings
with the shape Eq. (10) should be observed for a < r < p =aPe. Some theoretical studies assume D = 0 [31–33] and define the Péclet number as Pe = U/(aDr ); here also depletion
wings with the shape Eq. (10) should be observed below p .

IV. EXPERIMENTS

To confront our findings to real systems, we perform experiments with Janus particles, whose dynamics is similar tothat of ABP [10,34,35]. Our Janus particles are propelled bya vertical AC electric field and swim in a two-dimensional
horizontal plane [36,37] (Appendix D, Fig. 4(a) and Movie 1in the Supplemental Material). We estimate the experimental
values of U , D, and Dr in Appendix D3. This allows us toplace the experiments on the phase diagram (dark blue arrow
in Fig. 2), with a Péclet number Pe  90: Depletion wings are
expected.
The measured pair correlation is shown in Fig. 4(b)(additional figures for other values of the electric field are provided in Appendix D5). Depletion wings are present, showingthat they are a robust feature of pair correlations of selfpropelled particles at high activity. At short range, qualitative
differences with the numerical simulations can be noted, suchas depletion behind the active particle (x < −a, y = 0) before
correlations turn positive. These differences can be attributed
to unavoidable hydrodynamic or electrostatic interactions[11]. Depletion wings are present at larger distances, but theirdecay cannot be characterized due to insufficient statistics.
Whether or not these wings follow the scaling laws predicted
here or exhibit different exponents is an open question.
V. CONCLUSION

We have unveiled two regimes where active polar particles without alignment interactions have a pair correlation
with a self-similar shape and an algebraic decay characterized
by anomalous exponents. In presence of velocity-orientation
coupling [12,38,39] or alignment interactions [6,11,12], we
may expect distinct scaling laws to appear in the correlations.
As we have shown, the correlations have a rich structureeven without three-body interactions, it would therefore be
instructive to measure them in a dilute configuration first, andthen to see how they evolve as density increases.
ACKNOWLEDGMENTS

The authors thank K. A. Takeuchi for his hospitality
and useful discussions, J. Iwasawa for his help with experimental and image processing protocols, and D. Bartolo, O.
Dauchot, and A. Solon for insightful discussions. D.N. was
supported by JSPS KAKENHI Grants No. JP19K23422, No.
JP19H05800, and No. JP20K14426.
APPENDIX A: NUMERICAL SIMULATIONS
1. Description

We consider N particles (N = 5000) in a square periodic
box of size L = Nρ with ρ the density (L ∼ 300 at ρ = 0.05).
Initially the positions ri and the orientations θi of the particles
are assigned uniformly at random.
We use stochastic molecular dynamics and consider thefollowing Langevin equations:√ṙi = −∇ i
V (ri − r j ) + U êθi + 2Dηi ,
(A1)
θ̇i =

FIG. 4. Experiments. (a) Experimental image, indicating the detected positions and orientations of the Janus particles (red dots andarrows). Scale bar: 10 μm. (b) Pair correlation B(r).



j=i

2Dr νi ,(A2)

where U is the velocity, D the transitional diffusion, and Dr
the rotational diffusion. The mobility γ is set to 1. ηi and νi
are Gaussian white noises with unit variance. During a timeincrement t ( t √= 0.05 for the simulations presented√ in thearticle), the term 2Dr νi generates an increment υ 2Dr t
with υ a random number √generated from a standard normal
distribution. Similarly for 2Dηi .
We use the following soft-sphere potential, where the diameter a of a particle and the potential strength  are set tounity,
1(1 − r)2 if r  1,
(A3)
V (r) = 20otherwise.

012605-4

PAIR CORRELATION OF DILUTE ACTIVE BROWNIAN …

PHYSICAL REVIEW E 103, 012605 (2021)

FIG. 5. Simulations in the limit cases Dr = 0 and D = 0 for dilute hard particles. (a) No rotational diffusion. Correlation B(r)for ρ = 0.05,  : D : Dr : U = 50 : 0.1 : 0 : 10. (b) Rescaled cuts with the exponents predicted in the article. (c) No translational
diffusion. Correlation B(r) for ρ = 0.02,  : D : Dr : U = 50 : 0 : 0.1 : 10. (d) Rescaled cuts with the exponents predicted in thearticle.

We let the system evolve for a time teq ∼ 500 before starting to record the correlations. B(x) is then recorded with aspatial resolution x = 0.1, with one measure every tc = 5,
during an overall time period T ∼ 5 × 105 . The results areaveraged over 100 realizations of the system.
2. Simulations of dilute hard particles

The theoretical approach that we use is a linearized Dean
equation, a framework which is valid for weak interactions.
Consequently, the simulations that we compare to our theoretical results have relative values of potential strength andvelocity  : U = 1 : 10. This means that the particles are able
to interpenetrate. In this weak interaction regime, we are able
to test quantitatively the scalings that we obtain.
That being said, it is interesting to test the scaling exponents that we obtain on simulations in the regime of strong
interactions. We do it in Fig. 5 and find that indeed the cutsof the correlation function obey the predicted scalings bothin the limit of low rotational diffusion and in the limit of
low translational diffusion. This leads us to state that thescaling exponents that we found are robust. Note that thelimit curves are different from our prediction in the weakinteraction regime.

particles, and expect range of density to be wider for soft
particles. There are two ways to estimate where the dilute
regime sits. The first is to look at the low density branchof the motility-induced phase separation (MIPS) line on thephase diagram for ABP [17]. As MIPS is a many-body effect, the density should be far below this line to remain in
the dilute regime. For instance, for Pe = 100, ρMIPS  0.18
(φMIPS  0.14).
The second way is to look at measured pair correlations.
By definition, in the dilute regime the pair correlation does
not depend on the density. We plot the pair correlation for
densities ρ = 0.01, 0.02 and 0.05 for D = 0 and Pe = 100 on
Fig. 6. We observe an important difference between ρ = 0.02
and ρ = 0.05, while the difference between ρ = 0.01 andρ = 0.02 is much smaller. We conclude that the limit of thedilute regime for this Péclet number sits around ρ  0.02
(φ  0.016).
APPENDIX B: THEORETICAL COMPUTATIONS
1. Equation for the correlations
a. Dean equation for active Brownian particles

The Langevin equations for active Brownian particles readṙi = −∇ i

3. Extension of the dilute regime

Here we discuss the range of density ρ where the dilute
approximation that we use should be valid. We focus on hard
012605-5

V (ri − r j ) + U êθi +

j=i

θ̇i =



2Dr νi ,√2Dηi ,(B1)(B2)

ALEXIS PONCET et al.

PHYSICAL REVIEW E 103, 012605 (2021)

FIG. 6. Simulations in the limit case D = 0 for hard particles,  : D : Dr : U = 50 : 0 : 0.1 : 10, for densities ρ = 0.01, 0.02, 0.05.

where D is the translational diffusivity, Dr is the rotational
diffusivity, and ηi and νi are Gaussian white noises of unit
variance. In this Appendix, we set the particle diameter a = 1
and the friction coefficient γ = 1.
We define the density in position-orientation space
f (r, θ , t ) asf (r, θ , t ) =

N
∞



fi (r, θ + 2mπ , t ),dϕ = ∇ϕ · dri +fi (r, θ , t ) = δ(ri (t ) − r)δ(θi (t ) − θ ).

(B4)

We consider a smooth and fast-decaying test function
ϕ(r, θ ). By definition of fi , ∞dθ fi (r, θ , t )ϕ(r, θ ).
(B5)
ϕ(ri (t ), θi (t )) = dr−∞Then, the time derivative of ϕ(ri (t ), θi (t )) can be written in
two different ways:

 ∞
d
∂ fi
ϕ(ri (t ), θi (t )) = dr(r, θ , t )ϕ(r, θ )
dθ
dt
∂t
−∞

 ∞
= dr
dθ fi (r, θ , t )(dt )−1 dϕ(r, θ ),
−∞

(B6)
 2π

∂ϕ
1
1 ∂ 2ϕ 2dθdθi + ∇ 2 ϕ(dri )2 +∂θ
22 ∂θ 2 i∂∇ϕ · dri dθi∂θV (ri − r j ) + U êθi dt= ∇ϕ · −∇ i+

(B3)

i=1 m=−∞

J(r, θ , t ) = −D∇ f (r, θ , t ) − f (r, θ , t )with the differential dϕ given by the Itô formula [40],

(B7)

j

∂ 2ϕ
.
(B8)
∂θ 2
The differentials are computed from Eqs. (B1) and (B2) (weassume ∇V (0) = 0). Performing integrations by part and recalling that ϕ is arbitrary, one obtains
 N


∂ fi
∂ 2 fi2= D∇ fi + Dr 2 + ∇ fi
∇V (r − r j (t ))
∂t
∂ θ
j=1
+ Ddt∇ 2 ϕ + Dr dt

∂
2Dr ( fi νi ). (B9)∂θUsing Eq. (B3) and rearranging the noises like Dean [27],we finally obtain the following Dean equation for f (r, θ , t ),− U eθ · ∇ fi −√2D∇ fi · ηi −∂∂f (r, θ , t ) = −∇J(r, θ , t ) −K (r, θ , t ),∂t∂θwith the currentsdθ (∇V ∗ f )(r, θ , t ) + f (r, θ , t )U êθ − f 1/2 (r, θ , t )η(r, θ , t ),

(B10)

(B11)

0

∂
f (r, θ , t ) − f 1/2 (r, θ , t )ν(r, θ , t ).
(B12)
∂θ

The spatial convolution is defined by ( f ∗ g)(r) = dr f (r )g(r − r ). η and ν are Gaussian white noises with correlations
K (r, θ , t ) = −Drηα (r, θ , t )ηβ (r, θ , t ) = 2Dδ αβ δ(r − r )δ(θ − θ  )δ(t − t  ),(B13)

ν(r, θ , t )ν(r, θ , t ) = 2Dr δ(r − r )δ(θ − θ  )δ(t − t  ).

(B14)

particles. We writeb. Linearized Dean equation

The Dean equation for active Brownian particles is nonlinear with multiplicative noise. It is thus very difficult totackle. Our approximation consists in linearizing it aroundan homogeneous density ρ which is the average density ofρf (r, θ , t ) =
+
2π



ρ
φ(r, θ , t ).
2π

The field φ is assumed to be of order 1 in ρ.

012605-6(B15)

PAIR CORRELATION OF DILUTE ACTIVE BROWNIAN …

PHYSICAL REVIEW E 103, 012605 (2021)

At the lowest order, the Dean equation [Eq. (B10)] becomes linear with additive noise.

c. Equation for the correlations

We define the following correlations(B17)
C(r1 , r2 , θ1 , θ2 ) = φ(r1 , θ1 )φ(r2 , θ2 ),1
C(r, θ , θ  ) = {C(0, r, θ , θ  ) − δ(r)δ(θ − θ  )}.
ρ
(B18)



∂φ
∂2
2
= D∇ + Dr 2 − U êθ · ∇ φ
∂t
∂θ
 2π
ρ
∂ν
. (B16)+dθ  (∇ 2V ∗ φ)(θ  ) + ∇ · η +2π 0∂θOne checks that the definition of C(r, θ , θ  ) is consistent with Eq. (1) in the main text. We use Itô calculus to compute the time
evolution of C:
C(r1 , r2 , θ1 , θ2 , t + δt ) − C(r1 , r2 , θ1 , θ2 , t )= φ(r1 , θ1 , t )δφ(r2 , θ2 , t ) + δφ(r1 , θ1 , t )φ(r2 , θ2 , t ) + δφ(r1 , θ1 , t )δφ(r2 , θ2 , t ).
Computing the terms from the linearized Dean equation [Eq. (B16)], one shows that ∂t C(r1 , r2 , θ1 , θ2 ) = D ∇12 + ∇22 + Dr ∂θ21 + ∂θ22 − U êθ1 · ∇ 1 + êθ2 · ∇ 2 C(r1 , r2 , θ1 , θ2 )
 2π

ρ
+
dθ  ∇12V ∗ C(r1 , r2 , θ  , θ2 ) + ∇22V ∗ C(r1 , r2 , θ1 , θ  )2π 0
+ 2D∇ 1 ∇ 2 + 2Dr ∂θ1 ∂θ2 δ(r1 − r2 )δ(θ1 − θ2 ).

(B19)

(B20)



We use the invariance of the system by translation and write the equation in terms of C(r, θ , θ ),∂t C(r, θ , θ  ) = 2D∇ 2 + Dr ∂θ2 + ∂θ2 + U (êθ − êθ  ) · ∇ C(r, θ , θ  ) + 2∇ 2V (r) 2πρdθ  ∇ 2V ∗ [C(r, θ , θ  ) + C(r, θ  , θ  )].
+
2π 0

(B21)

The conventions are such that the pair correlation function C(r, θ , θ  ) is normalized with respect to the density ρ. Focusing
on the low density regime (ρ → 0), we can neglect the convolution of C with the potential V . The equation that we consider is∂t C(r, θ , θ  ) = 2D∇ 2 + Dr ∂θ2 + ∂θ2 + U (êθ − êθ  ) · ∇ C(r, θ , θ  ) + 2∇ 2V (r).
(B22)
ρ→0

Note that for a passive system, the solution of this equationis the direct correlation function (the one involved in the
Ornstein-Zernike equation). In the following, we are onlyinterested in the stationary correlations which satisfy the following linear partial differential equation:
2D∇ 2 + Dr (∂θ + ∂θ  ) + U (êθ − êθ  ) · ∇ C(r, θ , θ  )= −2∇ 2V (r).

(B23)

Let us write Eq. (B23)
in Fourier space, using the convention C̃(k, θ , θ  , t ) = dre−ikrC(r, θ , θ  , t ):[−2Dk 2 + Dr (∂θ + ∂θ  )]C̃(k, θ , θ  )= 2k 2Ṽ − iU k · (êθ − êθ  )C̃(k, θ , θ  ).

The passive correlation (U = 0) does not depend on theangles, the order 0 of the equation above leads to

By rotational invariance, C(r, θ , θ  ) = C(R−θ · r, 0, θ  −θ ). Like in the article, we define the “profile seen by a particle
in its reference frame,”1
B(r) =2π

 2πC(r, 0, θ  )dθ  .

(B24)

0

C̃ (0) (k) =−Ṽ (k).
D

2. Low activity−2Dk 2 + Dr ∂θ2 + ∂θ2 C̃ (1) (k, θ , θ  )1, one can expand C in power of U := −ik · (êθ − êθ  )C̃ (0) (k),C̃ (k, θ , θ ) = ik · (êθ − êθ  )D̃ (k),(1)

C(x, θ , θ  ) = C (0) (x) + UC (1) (x, θ , θ  )+ U 2C (2) (x, θ , θ  ) + . . .

(B27)This is the usual random phase approximation solution for
direct correlations. At order 1, the equation to solve and its
solution areAt low activity, U

(B26)

(B25)
012605-7

D̃(1) (k) =

(1)

(B28)
(B29)

(0)

C̃ (k).
2Dk 2 + Dr

(B30)

ALEXIS PONCET et al.

PHYSICAL REVIEW E 103, 012605 (2021)

Next, we solve the second order,−2Dk 2 + Dr ∂θ2 + ∂θ2 C̃ (2) (k, θ , θ  )= 2k 2 (1 − êθ · êθ  )D̃(1) (k),Ṽ (k)C̃(k, θ , θ ) =DṼ (k)
B̃(k) =
D




(B31)

U2
1
−1 + 2 2D 2k + U−2
−1 +

U2
1
D2 2k 2 + U−2−D̃(1) (k)k 2 D̃(1) (k)(êθ · êθ  ).
+D
Dk 2 + DrC̃ (2) (k, θ , θ  ) =(B32)

At the end of the day, the expansions at order U 2 ofC̃(k, θ , θ  ) and B̃(k) are

U ik · (êθ − êθ  ) U 2
11−+ 2 (êθ · êθ  )− 2−22
D 2k 2 + U−2
D
2k + Uk + U−2Uikx−,
2π D 2k 2 + U−2


,

(B33)
(B34)

√
with the typical length scale U = Dr /D.
Our goal is now to look at large distances, that is to say, small wave number k. We assume that Ṽ (k) is regular at 0 (short-range
potential) and we make the substitution Ṽ (k) → Ṽ (0). We define the function G̃(k) and its inverse Fourier transform G(r) as
1

G̃(k) =

2k2 + U−2⇔

G(r) =

r
1
,
K0 √
4π
2 U

(B35)

with K0 the modified Bessel function of the second kind of order 0. In real space, the expansion of B(r) becomesṼ (0)
Ṽ (0) ∂G(r)Ṽ (0)+ U 2 3 G(r) − U,
DD
2π D2 ∂x√Ṽ (0) DrṼ (0)
Ṽ (0)
r
r
+U √
+ U2
B(r, θ ) ∼ −
K0 √
K1 √
3
3/2
r→∞
D
4π D
2 U
2π 2D
2 U
B(r) ∼ −

(B36)

r→∞

where K1 is the modified Bessel function of the second kindof order 1. We detail the meaning of the terms. The first one
in the passive correlation. The second term, scaling as U 2 is apositive contribution to the isotropic part of the correlations.
As stated in the article, this has been argued, in the literature,to account for the motility-induced phase separation. The lastterm is the dipolar correlation, at order U , on which we focused in the article. Eq. (B37) corresponds to the predictionin Figs. 3(b) and 3(c). At large distance, the Bessel functions
decay exponentially,

π −r
e .
K0 (r) ∼ K1 (r) ∼
(B38)
r→∞
r→∞
2r
Both the dipolar contribution (cos θ ) and the additional
isotropic part thus decay exponentially over the length scaleU.

C̃(k, θ , θ  ) =−2k 2Ṽ (k)2Dk 2 − iU k · (êθ − êθ  )

.

D

,

2
−2 k 2
2k 2 − i −1
k
+
x
r
r

(B41)

with the characteristic length scale r = D/U . kx = k · êx ,ky = k · êy , where êx and êy are the unit vector along thehorizontal and vertical axes of the plane.
We consider the limit of large distance, that is to say k →0. We assume that Ṽ (k) is regular so that we can replace it byṼ (0). Furthermore, an analysis of leading terms gives ky2 ∼ kx3 .
We keep up these leading terms:−2Ṽ (0)kx2

B̃(k) ∼k→0

D

−2 k 2 − 4i −1 k 3
r
y
r
x

.

(B42)

We can inverse Fourier transform with respect to ky :B(kx , y) = −2Ṽ (0) r 2 
kx K0 |y| −4i r kx3 .
πD

(B43)

We then Fourier transform with respect to kx : B(x, y) =(2π )−1 dkx eikx x B(kx , y). We perform the changes of variables q = −( r y2 )1/3 kx , w = xu1/3 ( r y2 )−1/3 and obtainṼ (0) 1
x
F 1/3
,
2/3
D y2
r |y|
 ∞

1
e−iqw q2 K0 (2 iq3 )dq.
F (w) = 2
π −∞

(B40)

This solution is very similar to the one for a binary mixtureof particles forced, respectively, by U êθ and U êθ  [26]. We
integrate over θ  to obtain B, this givesk 2Ṽ (k) 2πdθ B̃(k) = −πD
2k 2 − i −10
r k · (êx − êθ  )(B37)−2k 2Ṽ (k)=3. No rotational diffusion

We now focus on the limit of no rotational diffusion Dr =0. We easily obtain the Fourier transform C(r, θ , θ  ) from
Eq. (B23),

2D∇ 2 + U (êθ − êθ  ) · ∇ C(r, θ , θ  ) = −2∇ 2V (r), (B39)cos θ ,

B(x, y) ∼

(B44)
(B45)

This is the scaling form mentioned in the article and used inFig. 3(e).

012605-8

PAIR CORRELATION OF DILUTE ACTIVE BROWNIAN …

PHYSICAL REVIEW E 103, 012605 (2021)

FIG. 7. Numerical resolution of Eq. (B48) with Ṽ (0) = U = Dr = 1, for k = 104 . (a) Real part of C. (b) Imaginary part of C. The solution
concentrates around (γ , γ  ) = (0, 0) and (π , π ). The numerical resolution consists in discretizing γ and γ  and then solving the linear systemcorresponding to Eq. (B48).
4. No translational diffusion

We set D = 0. Looking at distances large compared to theparticle diameter, we replace Ṽ (k) by Ṽ (0). The equation we
consider is 2
2k 2Ṽ (0)∂θ + ∂θ2 + i p k · êθ − êθ  C̃(k, θ , θ  ) =
,
Dr
(B46)
with the persistence length p = U/Dr .
At small p k (large distance compared to p ), the following
development can be obtained:C̃(k, θ , θ  ) =−Ṽ (0) 
1 + i p k · (êθ − êθ  )Dr+ 2p k 2 êθ · êθ  + O(( p k)3 ) .

(B47)

This is a hint that C̃(k, θ , θ  ) is analytical around k = 0,
meaning that C(r, θ , θ  ) decays fastly (e.g., exponentially) at
distances large compared to p .
We now consider distances below p but still large compared to the particle diameter. We define the angles γ andγ  in the reference frame of k: k · êθ = k cos γ and k · êθ  =k cos γ  , where k is the norm of k. We obtain 2
2k 2Ṽ (0),∂γ + ∂γ2  + i p k(cos γ − cos γ  ) C̃(k, γ , γ  ) =Dr(B48)
and we study it in the regime p k
1. A numerical resolution
at constant k (see Fig. 7) shows that C̃(k, γ , γ  ) concentrates
around the two points (γ , γ  ) = (0, 0) and (π , π ). We focus
on (0, 0) around which the equation readsi2Ṽ (0) 2∂γ2 + ∂γ2  − p k(γ 2 − γ 2 ) C̃(k, γ , γ  ) =
k .
2
Dr
(B49)
We realize that we can inject the following scalings:
2Ṽ (0)( p k)3/2 H̃[γ ( p k)1/4 , γ  ( p k)1/4 ].
Dr
pk 1
(B50)
From Eq. (B49), the function H̃ (u, v), for u and v unbounded,is independent of k and is the solution of the linear partialC̃(k, γ , γ  ) ∼differential equation
i∂u2 + ∂v2 − (u2 − v 2 ) H̃ (u, v) = 1.
2

(B51)

Around (γ , γ  ) = (π , π ), the scalings areC̃(k, γ , γ  ) ∼pk

1

2Ṽ (0)( p k)3/2 H̃ ∗ ((γ − π )( p k)1/4 ,Dr(γ  − π )( p k)1/4 ),(B52)

with H ∗ the complex conjugate of H.
The scaling for B̃ around γ = 0 is
 2π
1B̃(k, γ ) =dγ C̃(k, γ , γ  )
2π 0
∼

pk

HB (u) =Ṽ (0)
1 π Dr( p k)5/4 H̃B (γ ( p k)1/4 ), (B53) ∞dv H (u, v).

(B54)

−∞

And around γ = π , one checks thatB̃(k, π − γ ) ∼pkṼ (0)
1 π Dr( p k)5/4 H̃B∗ ((π − γ )( p k)1/4 ),(B55)

with HB∗ the complex conjugate of HB . We now switch frompolar coordinates (k, γ ) to cartesian coordinates (kx , ky ). We
approximate
+k if γ  0
kx = k cos γ ,−k if γ  πkγ
if γ  0
ky = k sin γ .
(B56)k(π − γ ) if γ  π
As we consider small angles (|γ |1 or |π − γ |1), thekx .
values of kx and ky are such that ky
Using the two expressions Eqs. (B53) and (B55), we obtaina scaling form for the Fourier transform B̃,012605-9B̃(kx ≷ 0, ky ) =1/4Ṽ (0)p ky,( p |kx |)5/4 H̃B±
π Dr
|kx |3/4(B57)

ALEXIS PONCET et al.

PHYSICAL REVIEW E 103, 012605 (2021)

with H̃B+ = H̃B used when kx > 0 and H̃B− = H̃B∗ used whenkx < 0. Finally, we perform the Fourier inversion
1
1
B(x, y) =dkx eikx x
dky eiky y B̃(kx , ky ), (B58)2π2πwhere HB± is the inverse Fourier transform of H̃B± , one has
HB− (a) = (HB+ )∗ (−a). This form corresponds to Eq. (10) of
the article. The rescaled cuts are plotted on Fig. 3(g) withthe prediction (gray line) computed from a rescaling of thenumerical solution of Eq. (7).

first with respect to ky , then with respect to kx . Using theappropriate changes of variables, we obtain a scaling form forB(x, y),1/3Ṽ (0) 4pp x,
Gπ Dr y4|y|4/3
 ∞1
G(w) =
dzeiwz z2 HB+ (z3/4 )
2π 0
 0
1
+
dzeiwz z2 HB− (|z|3/4 ),2π −∞B(x, y) =(B59)

APPENDIX C: NUMERICAL INTEGRATION
OF THE THEORETICAL EQUATION
1. Equation in terms of three parameters(B60)

We consider the time-dependent equation [Eq. (B22)] forthe correlations C at low density,∂t C(r, θ , θ  ) = 2D∇ 2 + Dr ∂θ2 + ∂θ2 + U (êθ − êθ  ) · ∇ C(r, θ , θ  ) + 2∇ 2V (r).

(C1)

In polar coordinates, we write r = reφ . The later equation depends on four coordinates (plus time): (r, φ, θ , θ  ). By performing
a rotation of angle θ , the symmetries allow us to reduce the problem to three parameters (r, α, β ),α = φ − θ,

β = θ − θ.

(C2)

We write x = reα = (x, y), C is then a function only of x and β. Its time evolution is given by


∂
∂
2
∂t C(x, β ) = 2D∇ + Dr Langles + U (1 − cos β ) − sin βC(x, β ) + 2∇ 2V (x),
∂x
∂y
∇2 =

∂2
∂2
+
,
∂x 2
∂y2

Langles = y2

(C4)

2
∂
∂
∂2
∂2
2 ∂
−x
−y
+
x
−
2xy
2
2
∂x
∂y
∂x∂y
∂x
∂y

It is important to note that B is given by the integration
over β, 2π
dβ C(x, β ).
(C6)B(x) =
0

2. Numerical integration

We consider the domain (x, y, β ) ∈ [−xmax , xmax ] ×[0, ymax ] × [−π , π ]. We discretize it with steps x in x,
y = x in y and β in β.
We start from C(x, y, β, t = 0) = 0 and integrate Eq. (C3)in time using an explicit Euler scheme with time step t. The
differential operators are evaluated using finite differences
valid at order ( x)2 and ( β )2 ,C(x +∂C(x, y, β ) =
∂x
∂2
C(x +
C(x, y, β ) =∂x 2
and so on.
The potential is
V (x) =

120x) − C(x −2 xx) + C(x −( x)2(1 − x)2x)+ 2 −y

∂
∂
∂
∂2
+x
+2 2.
∂x
∂y ∂β
∂β

(C5)

The boundary conditions are as follows:(1) Periodic boundary conditions for β;(2) C(−xmax , y, β ) = C(+xmax , y, β ) = C(x, ymax , β ) =0;(3) We use the symmetry relation C(x, −y, −β ) =C(x, y, β ) to impose the additional points C(x, − y, β ) =C(x, y, −β ).
The numerical integration in time converges to the stationary solution at sufficiently large time. When the increment
on C over a time step t is small enough, we output thestationary solution C eq (x, β ) and its integral over β, Beq (x).
This method is used on Figs. 2(a)– 2(d).

APPENDIX D: EXPERIMENTS,x) − 2C(x)if x  1,otherwise,

and its Laplacian is evaluated on the grid.

(C3)

1. Experimental system, (C7)(C8)

The experimental system is the one used in Ref. [37]. We
used a suspension of Janus particles of diameter a = 3.17 ±0.32 μm sandwiched between two ITO electrodes separated
by a spacer of size H = 50 μm. The particles are suspended
in a sodium chloride solution of concentration 10−4 mol/L.
The use of sodium chloride diminishes the temporal variationof the system, which enables long observation required for
calculating pair correlations. Furthermore, because the addition of sodium chloride decreases the width of the screening

012605-10

PAIR CORRELATION OF DILUTE ACTIVE BROWNIAN …

PHYSICAL REVIEW E 103, 012605 (2021)

FIG. 8. Determination of the experimental parameters. 670 particles are tracked for 35 frames on average (minimum: 20 frames).
(a) Histogram of the average velocities of the particles, the vertical black line is the mean. (b) Mean-square displacement. The estimatederror bars are smaller than the symbol size. (c) Mean-square angular displacement. Error bars represent standard errors estimated from the
distribution of θ 2 for each given t, the solid line is a linear fit over 0  t  1.5 s.

double layer (Debye length) and consequently the strengthof the induced hydrodynamic flow [41], both the electrostatic and hydrodynamic interactions between the particles
are diminished. This keeps the interactions indiscernible atthe frequency range we used here [37]. For the correlations
presented in the main text, we applied an electric field of
frequency f = 5 kHz and amplitude 2 × 105 Vpp /m (voltage
10 Vpp ) in the vertical direction. In this range of frequency,
the Janus particles move in the direction of the uncoatedhemisphere due to the induced-charge electrophoresis [36]and exhibit neither of electrostatic attractive interactions nor
strong hydrodynamic interactions [37]. Moreover, no polarorder is observed even at high density [36]; these experimental
optimization justifies the theoretical modeling by active Brownian particles as a first approximation.
Videos of the system were captured by a CMOS camera(Baumer, LXG-80, 3000 × 2400 pixels, 8 bit grayscale) at theframerate 10 fps mounted on an inverted microscope (Olympus, IX70) equipped with a 40× objective lens (LUCPLFLN,
NA=0.60). A green filter is inserted between the sample andthe halogen lamp to increase the contrast of the hemispheresof the Janus particles. The acquisition length is 14 min (8400
frames).
We append an experimental movie at 10 Vpp , captured at
10 fps. The field of view of this movie is cropped to 585 px ×
409 px (70 μm × 49 μm) for visibility. The movie is played
at the real speed.

2. Image analysis

Particles were detected using the Hough circle transform
algorithm implemented in the OpenCV library [42]. The positions of the particles are the centers of the circles. For eachdetected circle, we compute the center of mass of the pixels
within it (the weights are the pixels’ values) [12]. The orientation of the particle is defined as the direction of the vector
between the center of mass and the center of the circle. We
will later give an estimate of the precision of this measure.

3. Experimental parameters

We now estimate the parameters of our experimental systems.
a. Particle diameter. We captured the images with the resolution 0.12 μm/px. In the obtained images of the particle

with the diameter a = 3.17 μm, the particle diameter appearsa  26 px.
b. Density. There are on average 487 particles in a 3000 ×2400 pixels (∼360 μm × 290 μm) image. This gives the487
 0.0047 μm−2 and thenumber density ρ  360 μm×290
μm
2
area fraction φ = ρπ (a/2)  0.04.
c. Velocity. We use Trackpy [43] to obtain the trajectories.
The instantaneous velocities can be obtained by applying a
Savitzky-Golay filter [44]. From the instantaneous velocities,
we find an average velocity U = 56 ± 7 px/s [standard deviation given for different particles, see Fig. 8 (a)]. Thus,U  6.7 μm/s, that is to say of the order of two particle
diameters per second.
Alternatively, one can compute the mean square displacement as a function of time [Fig. 8(b)]. At the lowest orderin time,  x 2  = (U t )2 . We find U  55 px/s  6.6 μm/s.
This is consistent with the previous result.
d. Rotational diffusion. To measure the rotational diffusion
coefficient Dr , we compute the mean square angular displacement (MSAD) as a function of time. At short time, we expect2
 θ 2  = 2θerr+ 2Dr t,(D1)

where θerr is the error made in the detection of the orientationof a given particle on a given frame. At times larger than themean free time τfree  ( √1ρ − a)/U  1.7 s, hydrodynamic
or electrostatic interactions may affect the orientation of theparticles during collisions, resulting in a larger angular displacement. The MSAD is plotted on Fig. 8(c); fitting the
MSAD with Eq. (D1) for 0  t  1.5 s leads to Dr =0.123 ± 0.004 s−1 and θerr = 8.7 ± 0.3◦ .
e. Translational diffusion. This is the hardest quantity toevaluate. We give an estimate based on the Stokes-Einstein
relation,
D=

kB T
,
6π η(a/2)β(D2)

where η = 1.0 × 10−3 Pas−1 is the viscosity of water, T 
300 K is the temperature, and β is a correction factor due tothe proximity of the bottom electrode. Assuming that Faxen’slaw (Eq. (7-4.28) of Ref [45]) is valid for h (distance to thewall) of the order of a, we obtain β  3 (one checks that the012605-11

ALEXIS PONCET et al.

PHYSICAL REVIEW E 103, 012605 (2021)

FIG. 9. Experimental pair correlations for different values of the voltage: 6 Vpp , 8 Vpp , and 10 Vpp (left to right). The amplitudes of theelectric field are 1.2 × 105 , 1.6 × 105 , and 2 × 105 Vpp /m.

thermal fluctuations are negligible). At the end of the day,D

5. Variation of electric field

300 × 1.38 × 10−23
6π × 1.0 × 10−3 × 1.585 × 10−6 × 3

 0.05 μm2 s−1  (0.07a)2 /s.

(D3)

4. Correlations

Once we have detected the positions and the orientations,
we consider a given frame. For each particle far from theedges of the image, we consider every other particle andcompute its position in the reference frame of the orientationof the first particle. We put the result in bins of size x =y = 0.1a. After processing all the frames and normalizing
the bins, we obtain the correlation [Fig. 4(b) of the main text].

[1] J.-P. Hansen and I. R. McDonald, Theory of Simple Liquids, 3rded. (Academic Press, London, 2006).
[2] H. P. Zhang, A. Be’er, E. L. Florin, and H. L. Swinney, Collective motion and density fluctuations in bacterial colonies, Proc.
Natl. Acad. Sci. USA 107, 13626 (2010).
[3] D. Nishiguchi, K. H. Nagai, H. Chaté, and M. Sano, Long-range
nematic order and anomalous fluctuations in suspensions ofswimming filamentous bacteria, Phys. Rev. E 95, 020601(R)(2017).
[4] A. Cavagna, A. Cimarelli, I. Giardina, G. Parisi, R. Santagati,
F. Stefanini, and M. Viale, Scale-free correlations in starling
flocks, Proc. Natl. Acad. Sci. USA 107, 11865 (2010).
[5] A. Cavagna, I. Giardina, and T. S. Grigera, The physics offlocking: Correlation as a compass from experiments to theory,
Phys. Rep. 728, 1 (2018).
[6] A. Bricard, J.-B. Caussin, N. Desreumaux, O. Dauchot, and
D. Bartolo, Emergence of macroscopic directed motion in
populations of motile colloids, Nature (London) 503, 95(2013).
[7] D. Geyer, A. Morin, and D. Bartolo, Sounds and hydrodynamics of polar active fluids, Nat. Mater. 17, 789 (2018).
[8] D. Geyer, D. Martin, J. Tailleur, and D. Bartolo, Freezing a
Flock: Motility-Induced Phase Separation in Polar Active Liquids, Phys. Rev. X 9, 031043 (2019).
[9] J. Palacci, C. Cottin-Bizonne, C. Ybert, and L. Bocquet,
Sedimentation and Effective Temperature of Active Colloidal
Suspensions, Phys. Rev. Lett. 105, 088304 (2010).
[10] F. Ginot, I. Theurkauff, D. Levis, C. Ybert, L. Bocquet, L.
Berthier, and C. Cottin-Bizonne, Nonequilibrium Equation of

We varied the strength of the electric field in the experiments setting the voltage to 6 Vpp and 8 Vpp in addition to the
10 Vpp data that we report in the article. The velocity U of
the particles is known to scale as the square of the strengthof the electric field [36], thus the Péclet number ranges from
Pe = 32 (6 Vpp ) to Pe = 90 (10 Vpp ). The experimental correlations B(r) are show on Fig. 9. In all three cases thequalitative behavior is the same: The depletion is maximal in
two wings behind the particle. Experiments with 10 Vpp havebeen performed first, and when the system ages particles start
to stick to the substrate, explaining the appearance of positive
concentric circles in the correlations for 6 and 8 Vpp .

State in Suspensions of Active Colloids, Phys. Rev. X 5, 011004(2015).
[11] J. Zhang, R. Alert, J. Yan, N. S. Wingreen, and S.
Granick, Active phase separation by turning toward the crowd,
arXiv:2011.03175 (2020).
[12] J. Iwasawa, D. Nishiguchi, and M. Sano, Algebraic correlations
and anomalous fluctuations in ordered flocks of Janus particles
fueled by an AC electric field, arXiv:2011.14548 (2020).
[13] J. Bialké, H. Löwen, and T. Speck, Microscopic theory for thephase separation of self-propelled repulsive disks, Europhys.
Lett. 103, 30008 (2013).
[14] A. P. Solon, J. Stenhammar, R. Wittkowski, M. Kardar, Y. Kafri,
M. E. Cates, and J. Tailleur, Pressure and Phase Equilibria
in Interacting Active Brownian Spheres, Phys. Rev. Lett. 114,
198301 (2015).
[15] U. Marini Bettolo Marconi and C. Maggi, Towards a statisticalmechanical theory of active fluids, Soft Matter 11, 8768 (2015).
[16] M. E. Cates and J. Tailleur, Motility-induced phase separation,
Annu. Rev. Condens. Matter Phys. 6, 219 (2015).
[17] P. Digregorio, D. Levis, A. Suma, L. F. Cugliandolo, G.
Gonnella, and I. Pagonabarraga, Full Phase Diagram of Active
Brownian Disks: From Melting to Motility-Induced Phase Separation, Phys. Rev. Lett. 121, 098003 (2018).
[18] J. U. Klamser, S. C. Kapfer, and W. Krauth, Thermodynamic
phases in two-dimensional active matter, Nat. Commun. 9, 5045
(2018).
[19] T. F. F. Farage, P. Krinninger, and J. M. Brader, Effective
interactions in active Brownian suspensions, Phys. Rev. E 91,
042310 (2015).

012605-12

PAIR CORRELATION OF DILUTE ACTIVE BROWNIAN …

PHYSICAL REVIEW E 103, 012605 (2021)[20] É. Fodor, C. Nardini, M. E. Cates, J. Tailleur, P. Visco, and F.
van Wijland, How Far from Equilibrium is Active Matter? Phys.
Rev. Lett. 117, 038103 (2016).
[21] S. C. Takatori, W. Yan, and J. F. Brady, Swim Pressure: Stress
Generation in Active Matter, Phys. Rev. Lett. 113, 028103(2014).
[22] A. Härtel, D. Richard, and T. Speck, Three-body correlations
and conditional forces in suspensions of active hard disks, Phys.
Rev. E 97, 012606 (2018).
[23] J. Bickmann and R. Wittkowski, Predictive local field theory forinteracting active Brownian spheres in two spatial dimensions,J. Phys.: Condens. Matter 32, 214001 (2020).
[24] A. Meyer, A. Marshall, B. G. Bush, and E. M. Furst, Laser
tweezer microrheology of a colloidal suspension, J. Rheol. 50,
77 (2006).
[25] V. Démery, O. Bénichou, and H. Jacquin, Generalized Langevinequations for a driven tracer in dense soft colloids: Construction
and applications, New J. Phys. 16, 053032 (2014).
[26] A. Poncet, O. Bénichou, V. Démery, and G. Oshanin, Universal
Long Ranged Correlations in Driven Binary Mixtures, Phys.
Rev. Lett. 118, 118002 (2017).
[27] D. S. Dean, Langevin equation for the density of a system
of interacting Langevin processes, J. Phys. A: Math. Gen. 29,
L613 (1996).
[28] F. D. C. Farrell, M. C. Marchetti, D. Marenduzzo, and J.
Tailleur, Pattern Formation in Self-Propelled Particles with
Density-Dependent Motility, Phys. Rev. Lett. 108, 248101(2012).
[29] D. S. Dean and R. Podgornik, Relaxation of the thermal Casimirforce between net neutral plates containing Brownian charges,
Phys. Rev. E 89, 032117 (2014).
[30] V. Démery and D. S. Dean, The conductivity of strong electrolytes from stochastic density functional theory, J. Stat.
Mech.: Theory Exp. (2016) 023106.
[31] A. P. Solon, M. E. Cates, and J. Tailleur, Active Brownian
particles and run-and-tumble particles: A comparative study,
Eur. Phys. J.: Spec. Top. 224, 1231 (2015).
[32] U. Basu, S. N. Majumdar, A. Rosso, and G. Schehr, Active
Brownian motion in two dimensions, Phys. Rev. E 98, 062121(2018).
[33] L. Caprini, U. Marini Bettolo Marconi, and A. Puglisi,
Spontaneous Velocity Alignment in Motility-Induced Phase
Separation, Phys. Rev. Lett. 124, 078001 (2020).

[34] J. R. Howse, R. A. L. Jones, A. J. Ryan, T. Gough, R.
Vafabakhsh, and R. Golestanian, Self-Motile Colloidal Particles: From Directed Propulsion to Random Walk, Phys. Rev.
Lett. 99, 048102 (2007).
[35] T. Mano, J.-B. Delfau, J. Iwasawa, and M. Sano, Optimal runand-tumble–based transportation of a Janus particle with active
steering, Proc. Natl. Acad. Sci. USA 114, 2580 (2017).
[36] D. Nishiguchi and M. Sano, Mesoscopic turbulence and local
order in Janus particles self-propelling under an ac electric field,
Phys. Rev. E 92, 052309 (2015).
[37] D. Nishiguchi, J. Iwasawa, Hong-Ren Jiang, and M. Sano, Flagellar dynamics of chains of active Janus particles fueled by an
AC electric field, New J. Phys. 20, 015002 (2018).
[38] J. Deseigne, O. Dauchot, and H. Chaté, Collective Motion ofVibrated Polar Disks, Phys. Rev. Lett. 105, 098001 (2010).
[39] O. Dauchot and V. Démery, Dynamics of a Self-Propelled Particle in a Harmonic Trap, Phys. Rev. Lett. 122, 068002 (2019).
[40] B. Oksendal, Stochastic Differential Equations: An Introduction
with Applications (Springer Science & Business Media, Berlin,2013).
[41] T. M. Squires and M. Z. Bazant, Breaking symmetries in
induced-charge electro-osmosis and electrophoresis, J. Fluid
Mech. 560, 65 (2006).
[42] The OpenCV Reference Manual, OpenCV, 2nd ed. (2019), http://www.opencv.org.cn/opencvdoc/2.3.1/opencv2refman.pdf.
[43] D. Allan, C. van der Wel, N. Keim, T. A. Caswell, D. Wieker,
R. Verweij, C. Reid, Thierry, L. Grueter, K. Ramos, apiszcz,
zoeith, R. W. Perry, F. Boulogne, P. Sinha, pfigliozzi, N.
Bruot, L. Uieda, J. Katins, H. Mary, and A. Ahmadia, softmatter/trackpy: Trackpy v0.4.2 (Zenodo, 2019), https://doi.org/
10.5281/zenodo.3492186.
[44] P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T.
Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser,
J. Bright, S. J. van der Walt, M. Brett, J. Wilson, K. Jarrod
Millman, N. Mayorov, A. R. J. Nelson, E. Jones, R. Kern,
E. Larson, C. J. Carey et al., SciPy 1.0: Fundamental algorithms for scientific computing in Python, Nat. Methods 17, 261(2020).
[45] J. Happel and H. Brenner, Wall effects on the motion of a single
particle, in Low Reynolds Number Hydrodynamics: With Special
Applications to Particulate Media, Mechanics of Fluids andTransport Processes, Vol. 1 (Springer, Dordrecht, 1983), pp.
286–357.

012605-13