Re-entrant percolation in active Brownian hard disks

Non-equilibrium clustering and percolation are investigated in an archetypal model of two-dimensional active matter using dynamic simulations of self-propelled Brownian repulsive particles. We concentrate on the single-phase region up to moderate levels of activity, before motility-induced phase separation (MIPS) sets in. Weak activity promotes cluster formation and lowers the percolation threshold. However, driving the system further out of equilibrium partly reverses this effect, resulting in a minimum in the critical density for the formation of system-spanning clusters and introducing re-entrant percolation as a function of activity in the pre-MIPS regime. This non-monotonic behaviour arises from competition between activity-induced effective attraction (which eventually leads to MIPS) and activity-driven cluster breakup. Using an adapted iterative Boltzmann inversion method, we derive effective potentials to map weakly active cases onto a passive (equilibrium) model with conservative attraction, which can be characterised by Monte Carlo simulations. While the active and passive systems have practically identical radial distribution functions, we find decisive differences in higher-order structural correlations, to which the percolation threshold is highly sensitive. For sufficiently strong activity, no passive pairwise potential can reproduce the radial distribution function of the active system.


Introduction
Active Brownian particles (ABPs) are one of the canonical models for describing and studying the properties of active matter.In their basic form, ABPs are spheres or disks that are self-propelled by a force of constant magnitude.Energy is dissipated by friction, and the motion (both translational and rotational) is subject to stochastic kicks.Like all active matter, systems of ABPs are intrinsically out of equilibrium.They are typically studied under steady-state conditions defined by the level of particle motility.
Despite the simplicity of the ABP model, it exhibits rich phase behaviour.2][3] This phenomenon, known as motility-induced phase separation (MIPS), occurs at sufficiently high mean density and activity in two 1,[4][5][6] and three dimensions. 1,7It is a consequence of the symmetry-breaking induced by the activity and the collisions between particles.However, even before MIPS takes place, a significant degree of clustering occurs, 8 to the extent that macroscopic networks of particles can arise without phase separation, as already observed in equilibrium systems with conservative attractive forces. 9,10As far as we are aware, the formation of these pre-MIPS structures in purely repulsive ABPs has not been studied in detail yet.
The mean cluster size in a given system depends on conditions such as temperature and density.At the transition from a distribution of discrete clusters to a system-spanning network, the mean cluster size diverges, defining the percolation threshold.5][16] In the latter case, the addition of conducting "filler" nanoparticles to an insulating matrix causes a large and sudden increase in the conductivity of the composite material when the fillers form a network that connects opposite boundaries of a macroscopic sample.In soft matter, percolation is one of the signatures of gelation, which can influence rheological properties and the process of phase separation. 9,17,18The ubiquity of colloidal gels in the food and pharmaceutical industries therefore makes percolation a phenomenon of wide-ranging importance in soft matter. 195][26] For example, increasing the aspect ratio of hard nanoparticles from spheres to long rods (such as carbon nanotubes) causes the packing fraction at the percolation threshold to decrease dramatically, but a simple inverse relationship between aspect ratio and the critical density is only reached at very high aspect ratios. 27rcolation theory has mostly been applied to explain the physical behaviour of equilibrium systems, but there are cases where the formation of system-spanning networks is important out of equilibrium.For example, the manufacture of composite materials often involves stirring for dispersal of the filler, 28 and flow or spin processes 29 for moulding.Any non-equilibrium structure produced by these processes is "frozen" into the network if the material is then solidified.Shearing a system containing carbon nanotubes causes them to align with the direction of the flow, impairing connectivity and raising the percolation threshold. 30,31owever, any local aggregation caused by the shear can have the opposite effect. 32When considering spherical particles, where alignment is not a consideration, the percolation threshold can either increase or decrease depending on the ratio of hard-core to connectivity distance and the strength of the shear flow. 33n the case of active matter, Levis and Berthier predicted a steady enhancement of percolation (i.e., a lowering of the critical density) with increasing activity in a kinetic Monte Carlo model of ABPs, albeit with a definition of pairwise connectivity that increased in range along with the level of motility. 34More recently, Sanoria and coworkers examined percolation in a suspension of soft ABPs, demonstrating that the softness of the interactions between particles has a major influence on structure and phase behaviour. 35For disks where the repulsive force is only linear in the pairwise separation, these authors concluded that increasing activity is equivalent to increasing softness.At high activity, the dense MIPS state gives way to a percolating porous network, which is then destroyed by a further increase in activity.Most recently of all, Caporusso and coworkers have studied the phase diagram of active dumbbells in three dimensions with strong, short-ranged attraction. 36Rich phase behaviour arises in this system due to the additional ingredient of particle shape on top of activity and the dimensionality of the space.Percolation occurs both in an arrested gel phase in the passive limit (zero activity) due to the sticky conservative attraction, and at sufficiently high density and activity, where the motility of the particles breaks up large clusters.
Considering pusher-type microswimmers-a different class of active particles-a percolation transition has been detected with increasing density at a specific level of activity, not only in suspensions of squirmers but also of asymmetric dumbbells. 37n the case of run-and-tumble particles, simulations on a twodimensional lattice have shown re-entrant percolation as a function of translational diffusion rate at low reorientation rates. 38 system of passive colloidal particles in a liquid of active chiral E. coli exhibits the characteristics of a percolation phase transition, albeit with critical exponents differing from those expected of the equilibrium case. 39Colonies of gliding M. xanthus bacteria exhibit collective motion once they have self-assembled into large clusters-effectively a percolation transition arising from the bacteria's motility. 40Percolation is also important in the context of electrical signalling networks in B. subtilis, where the percolation threshold for signal transmission across a community coincides with the optimal balance between the cost of firing to individual cells and the benefit to the bacterial community. 41he aims of the present article are two-fold.Firstly, we study the mechanism of percolation before MIPS appears in a two dimensional suspension of repulsive ABPs, unravelling how the percolation threshold changes as activity is increased.Secondly, we explore whether the effect of activity on percolation can be modelled as a small perturbation of a system of passive particles.Our route to an effective passive model is an adaptation of the iterative Boltzmann inversion technique, 42,43 which is usually used to devise potentials for coarse-grained models in equilibrium, starting from more detailed atomistic simulations.

Brownian dynamics of active particles
We simulate a system of N disks of diameter σ in a twodimensional square box with edge L and periodic boundary conditions in both directions.The packing fraction cannot be precisely defined because the effective particle diameter depends on activity, while the conventional Barker-Henderson expression 44 for effective diameter relies on equilibrium conditions.We therefore use the number density ρ = N/L 2 as a parameter to change the state of the system. 3,45ctive particles move via Brownian dynamics with a selfpropelling force on each particle that has its own rotational dynamics.The translational and rotational equations of motion for this model are where is the conservative force acting on particle i and V (r i j ) is the inter-particle pair potential.k B is Boltzmann's constant and T the absolute temperature.F a is the (constant) modulus of the self-propulsion force acting along the orientation vector ⃗ n i , which forms an angle θ i with the positive xaxis, and D t and D r are the translational and rotational diffusion coefficients, respectively.The components of the random forces ⃗ ξ i and ξ i,θ are white noise with zero mean and delta-correlation: , where α, β each represent x or y components, and ⟨ξ θ (t)ξ θ (t ′ )⟩ = δ (t − t ′ ).For equilibrium systems, the two diffusion coefficients follow a Stokes-Einstein relation for spherical particles (with diameter σ ) 46 : D r = 3D t /σ 2 .8][49][50][51] In our work we will treat D t and D r as independent control variables in line with previous numerical studies. 3,6,52or the repulsive core we use the pseudo-hard sphere 53,54 (PHS) model, which has been shown to reproduce equilibrium 55 and out-of-equilibrium 56,57 characteristics of hard-sphere suspensions, while being continuously differentiable.The PHS pairwise potential has the form (3) where r is the centre-to-centre distance, σ is the particle diameter and ε sets the energy scale.We will express all lengths and energies in terms of σ and ε, respectively.Introducing m as the mass of a disk, we may also define a natural unit of time τ PHS = mσ 2 /ε.
Our goal is to study clustering and percolation before MIPS appears.To generate configurations for the percolation analysis, we simulate the system with LAMMPS 58 (Large-scale Atomic/Molecular Massively Parallel Simulator).We initiate the system with particles located on a square lattice and run for 10 6 steps with a time step of 5×10 −6 τ PHS to ensure the system reaches steady state.The simulation is run for a further 2 × 10 8 steps, saving the particles' coordinates frequently, but with sufficient time between frames to ensure that the sampled configurations are independent.As a measure of the degree of activity, we use the Péclet number Pe, i.e., the dimensionless ratio between advective and diffusive transport, defined as: where v = F a D t /k B T is the self-propulsion velocity and τ r = 1/D r the reorientation time.For the results presented in the main article, the activity is varied by changing D r while fixing The choice of these values is based on the work of Stenhammar 1 and other studies that deploy the potential in Eq. ( 3). 3 In ESI Section 1, we test the effect of controlling Pe by varying F a rather than D r .

Monte Carlo of passive particles
We will compare the structure of the active system to that of a passive (equilibrium) model with an effective interaction potential that is to be deduced in Section 3.2.The passive system is most easily simulated using standard canonical Monte Carlo methods.Starting from an initial configuration of particles located on a square lattice, we perform stochastic trial displacements and test for acceptance of the move using the Metropolis algorithm, 59 with a target of 50% of the moves being accepted.We simulate for 5 × 10 4 MC sweeps to equilibrate and then gather statistics over a further 5 × 10 5 MC sweeps, saving the configuration every 1000 sweeps to produce 500 independent configurations.

Detecting percolation
The percolation threshold is defined as the critical set of conditions at which the mean cluster size diverges.Beyond this point in the macroscopic limit, an "infinite", system-spanning cluster is always present.In our simulations, a cluster is defined as a group of connected disks.In keeping with general studies of continuum percolation, two disks are taken to be connected when their centres are separated by less than a threshold, 24 which we take to be λ σ , setting the range parameter λ = 1.3.This value of λ allows a large range of densities and Péclet number to be studied without entering the MIPS region 3 .The sensitivity of the results to λ is tested in ESI Section 2.
To detect percolating clusters, we use a strict wrapping criterion, 60,61 where a cluster is only counted as percolating if it connects to its own images through the periodic boundary conditions of the simulation cell, thereby generating an infinite structure.This is a more stringent requirement than the cluster simply having a dimension large enough to connect opposite sides of the cell.Wrapping clusters are detected by building a contiguous image of the cluster that contains exactly one image of each disk in the cluster, and then testing for pairs of particles that are not directly in contact in the isolated cluster, but become connected when the periodic boundary conditions are considered.If the cluster contains one or more such pairs, then it wraps through the periodic boundaries, otherwise it is isolated and finite. 10he percolation probability, R(L, ρ) is the probability that a randomly chosen configuration in a trajectory contains a percolating cluster in any direction in a simulation cell of linear dimension L. This is calculated by counting the number of configurations that contain a wrapping cluster and dividing by the total number of configurations in the trajectory.
In the macroscopic limit, the percolation transition is characterised by a step change in R(L, ρ) from 0 to 1 at a critical density ρ = ρ c .However finite-size effects broaden this transition in a well-defined way for the equilibrium case.One way to establish the value of ρ c is to find the density at which the curves of R(L, ρ) as a function of ρ intersect for different values of L. Alternatively, lines of constant R can be extrapolated to infinite L, where they should converge at ρ c .How this is done in practice is detailed in the next section.

Percolation of active disks
For a given value of the Péclet number and linear box dimension L, we locate the range of density over which the percolation probability R(L, ρ) rises from 0 to 1.The resulting curves for 0 ≤ Pe ≤ 40 at L = 200σ are shown in Fig. 1a.We immediately see that the introduction of a small amount of activity causes a shift in the percolation response to lower density, meaning that activity is promoting percolation.Increasing the activity further, the trend continues up to approximately Pe = 7, where the shift stops, and increasing the activity further causes the transition to widen.
For each Péclet number, we then perform simulations at four box sizes, L/σ = 100, 141, 200, 282.point at which all these curves intersect.In passive systems, 62 the width of the response R(L, ρ) is proportional to L −1/ν where the value of the critical exponent ν in two dimensions is 4/3, and curves for different L should collapse onto each other if plotted as a function of the scaled density x = (ρ − ρ c )L 1/ν .The inset of Fig. 1b shows that these equilibrium scaling and collapse properties still apply in our weakly active, non-equilibrium system of ABPs.Preservation of the equilibrium values of critical exponents has also been observed in very soft active particles in a different regime of activity in recent work by Sanoria et al. 35 To calculate the percolation threshold, we identify the density at which each percolation probability curve equals a series of probabilities, R(L, ρ) = p for 0.2 ≤ p ≤ 0.8 in steps of ∆p = 0.1, by linear interpolation between measured points on the curves.All these densities should converge on the percolation threshold ρ c as L increases, since R tends to a step function, lim L→∞ R(L, ρ) = Θ(ρ − ρ c ). Fig. 2a shows lines of density at constant R(L, ρ) = p as a function of L −1/ν with ν = 4/3.Extrapolation of the lines for different p to L −1/ν = 0 gives the percolation threshold ρ c in the macroscopic limit.We take the spread of the intercept values as a measure of the statistical uncertainty in ρ c .
We use this method to calculate the percolation threshold as a function of Péclet number, as reported in Figure 2b.As shown in the figure, a small amount of activity rapidly reduces the percolation threshold until a minimum is reached around Pe = 10, after which increasing the activity further results in an increase in the percolation threshold.This nonmonotonic behaviour creates re-entrant percolation as a function of increasing activity for 0.5365 < ρσ 2 < 0.5455.Analogous results for ABPs where activity is controlled by varying the propulsion force F a rather than the rotational diffusion constant D r are presented in the Electronic Supplementary Information (ESI).
We stop the percolation analysis at Pe = 40, which is just below the MIPS boundary. 3After phase separation, droplets of the dense phase are always fully connected, so the nature of percolation changes completely to a question of whether the droplet itself is large enough to span the simulation cell.
The mean coordination number ⟨n⟩ over all disks is far from being an invariant at the percolation transition, as shown in the inset of Figure 2b, rising by almost 30% over the range of activities studied.The rise is steepest when weak activity is introduced to the reference equilibrium model.Tracking ⟨n⟩ along a line of constant density ρ = 0.54σ −2 rather than along the percolation threshold, we still see a monotonic rise, even as the system leaves the re-entrant region and stops percolating at higher Pe (Fig. S6 of the ESI).Hence, the increasing mean coordination of individual particles reflects greater internal networking within clusters rather than connectivity of the system on a larger scale.
The evolution towards intra-cluster connections implies a growing local inhomogeneity of the fluid structure.This effect is visible as an increasingly patchy texture with activity (but without phase separation) in the snapshots in Fig. 3.The local density at a given particle can be quantified by the inverse area of its cell in a Voronoi tesselation of the system.The increasing patchiness registers as a broadening of the distribution of this local density with Péclet number (see ESI Fig. S9).
As in any case of continuum percolation, decreasing the connectivity range λ shifts the percolation threshold to higher densities, since fewer pairs qualify as neighbours.As shown in Fig. S2 of the ESI, we observe pre-MIPS re-entrant percolation down to λ = 1.2.Below this range, the minimum in the percolation threshold is driven towards the density of the MIPS transition, which then intervenes.Conversely, re-entrance persists on increasing the connectivity range at least as far as λ = 1.4,with the minimum occurring at decreasing Pe as the range lengthens.

Effective interaction potentials
To understand the re-entrant behaviour of the percolation threshold in Fig. 2b, we test whether these weakly active systems can be modelled as a perturbation of a passive, equilibrium system with an effective potential to account for the structural effects of activity.Schemes for extracting effective potentials specifically for active matter have been devised by other authors, including approximate first-principles approaches by Farage et al. 63 (further investigated by Rein and Speck 64 ) and by Cameron et al. 65  Garcia. 67 the context of ABP percolation, we want to capture the structural correlations of the pre-MIPS steady state in as much detail as possible rather than deriving an equation of state or predicting the MIPS transition.With these priorities in mind, we use a modified iterative Boltzmann inversion (IBI) method 42,43 to find effective pairwise potentials from the radial distribution functions (RDFs) measured in the active simulations.IBI is usually used to deduce coarse-grained potentials for equilibrium simulations, starting from more detailed atomistic simulations.It relies on Henderson's theorem, 68 which shows that if a RDF is the result of pairwise interactions between particles then the underlying pair potential V (r) is unique.IBI therefore compares the RDF of the true system being modelled with the RDF returned from a simulation using a trial potential, expressed in piecewise linear form.The difference between the two RDFs is used to modify the trial potential.Further iterations of comparison and modification should lead to convergence of the two RDFs and (according to Henderson) the correct pair potential.To the best of our knowledge, IBI has not previously been used to extract effective pairwise potentials for active matter.
We already know the conservative PHS contribution of Eq. ( 3) to the disk interactions, so we split our putative effective potential V eff (r) into a sum of the PHS potential and the contribution V active (r) from the activity: The unknown V active is the only part of the potential that is changed via the IBI procedure as follows: where g(r) and g target (r) are the radial distribution functions obtained from a Monte Carlo simulation using the current trial potential V eff (r) and the dynamic simulations of the active system respectively.V active (r) is tabulated on a grid with spacing δ r = 0.01σ , matching the bin width of g(r).Linear interpolation is used between the tabulated points, while V PHS (r) is evaluated directly from Eq. ( 3) to obtain V eff (r) for the Monte Carlo simulations.We truncate the potentials at r = 3σ , which is large enough to capture the important structural information and to avoid a prominent discontinuity at the cutoff.It has been shown that increasing the cutoff in IBI does not necessarily capture any additional information and can make convergence worse. 43 assist convergence, we initialise V active (r) to 0 for the lowest activity level (Pe = 1) and then seed each successively higher level of activity with the converged result of the previous case.This approach provides the best possible starting point for each optimisation and leads to better convergence than an initial approximate inversion of g(r), as usually used in IBI in the context of coarse-grained potentials.We terminate the iterative process when the total squared deviation between the target active and effective RDFs stops decreasing and reaches a decisive plateau (see ESI Fig. S3).
Figure 4a shows a comparison of g(r) for Pe = 2 at ρ = 0.54σ −2 , R(L,ρ) and the result of the modified IBI procedure.After 35 iterations, the two RDFs are indistinguishable within the residual statistical noise.Hence, IBI is indeed able to return a pair potential that accurately reproduces the RDF of the active fluid.Figure 4b shows the resulting potentials V eff (r) from the modified IBI procedure for four different Péclet numbers at ρ = 0.54σ −2 .The effective potentials share common features, with a sharp attractive well at contact, which deepens with increasing activity as expected.
There is also a shallower secondary feature at r = 2σ before a tail to 0, which is related to the second peak in g(r) at that radius.This feature already distinguishes the ABP effective potential from the potential of a simple passive fluid such as Lennard-Jones (LJ), where g(r) may also show several peaks but the potential only has a single minimum.The second and later peaks in the LJ (and similar) radial distribution functions are the result of correlations between successive shells of neighbours, and no additional features are required in the potential energy function at the positions of the peaks in g(r).IBI correctly returns a smoothly decaying LJ potential with increasing r when presented with a multi-peak g(r) to invert under passive conditions. 43In contrast, the fact that the ABP radial distribution function generates a second feature in the effective potential implies that indirect correlations of nextnearest neighbours interacting through the primary minimum of V eff (r) would not be sufficient to reproduce the detailed structure of the active fluid, even at the weakest levels of activity studied here.
Beyond Pe = 10, IBI fails to return a V eff (r) that accurately reproduces g target (r).Henderson's theorem guarantees that a given RDF originates from a unique pair potential if only pairwise interactions are at play. 68However, it does not guarantee that any RDF can be produced by a pairwise potential.Our simulations of ABPs at moderate levels of activity provide examples of RDFs for which a generating pair potential does not seem to exist.We may infer that at Pe = 10, our ABP system is already showing the effects of many-body and/or directional interactions in the RDF (which is still a pairwise correlation function).These effects cannot be fully imitated in a system with conservative, isotropic, pairwise interactions.Unfortunately, the apparent nonexistence of a faithful V eff (r) as MIPS is approached means that we cannot gain insight into that transition through the evolution of the effective potential.
Using the effective potentials for weak activity, where the RDF is accurately reproduced, we perform larger Monte Carlo simulations and calculate the percolation probabilities, which are shown in Fig. 4c (open symbols).The simulations using the effective potentials produce a decisively different percolation response from their active counterparts, even for the weakest activity studied of Pe = 1.For Pe = 1 and 2, percolation in the passive simulations is shifted to higher density, whereas for Pe = 10 it is shifted to lower density.We have also calculated other structural quantities, such as the average coordination number, radius of gyration and cluster size distribution, which show small deviations when comparing the passive and active systems.The details of these calculations and the results are shown in the ESI.
The comparison between the active and passive cases immediately shows that percolation is extremely sensitive to details of a fluid's structure.Each pair of passive and active simulations in Fig. 4c has practically identical RDFs, but distinct clustering.It is already known that RDFs can be insensitive to the precise form of the potential, 69 despite the formal uniqueness theorem. 68Here, however, we are comparing other structural aspects of fluids where the RDFs are already faithfully reproduced.The differences in percolation indicate that there must be structural differences in the two systems beyond the pairwise correlation function g(r).
Following the work of Torquato et al. and Wang et al., 69,70 we calculate the conditional nearest neighbour distribution, G V (r), which provides a signature of any many body effects in the system.2πrρG V (r) dr is the probability that, given a circular region containing no particle centres, a particle centre lies in a shell of thickness dr around this empty region.By construction, G V (r) up to the radius of a hard particle (approximately r = σ /2 in our model) is determined solely by the particle density, since the density sets the probability that a point inside a disk (which is inaccessible to the centres of other particles) has been selected.Beyond that point, G V (r) can in principle be expressed in terms of many-body correlation functions. 70ESI Fig. S10 shows that there is a distinct difference between G V (r) for the active and passive systems, despite the "structural degeneracy" of the systems at the pairwise level. 71The deviation of the G V (r) curves implies that many-body correlations exist in the structure of ABPs from very low Péclet number and therefore that the structural effects of even weak activity cannot be captured in full by an effective pairwise potential.
Setting aside the detailed mapping of active disks onto an effective passive system, we may still ask whether re-entrant percolation can arise from simple pairwise attraction at equilibrium as a function of the strength of the attraction with respect to the thermal energy.We have performed equilibrium MC simulations on square-well disks with range λ SW σ and depth ε SW , to test how clustering varies with temperature.Setting both the range of the attraction λ SW and the connectivity criterion λ equal to 1.3 (matching our ABP simulations), we find a monotonically decreasing percolation threshold with decreasing temperature until gas-liquid phase separation intervenes (see ESI Fig. S4).This qualitative behaviour persists for different values of λ = λ SW .Nevertheless, for an estimate of the square-well ranges corresponding to the activity-induced potentials in Fig. 4b, we have evaluated the second virial coefficient B 2 and determined the value of λ SW that would return the same value if the well depths are also matched.Beyond about Pe = 2, the effective square-well range settles down at λ SW = 1.23 (ESI Section S5).
Data in an early study of percolation of square-well disks 72 imply that re-entrance might arise if the range of attraction is decoupled from the connectivity distance and set to a significantly shorter value, λ SW < λ , although the authors did not draw attention to this point.We have investigated this possibility with the benefit of more powerful computers and the full treatment of finite-size scaling as deployed earlier in the present work, but have not found any combination of λ and λ SW that gives rise to re-entrance.Figure S4 of the ESI shows some example results for λ SW < λ at λ = 1.3.Hence, the re-entrance seen in ABPs does seem to be an inherent and non-trivial result of the activity.

Conclusions
Our simulations show that weak activity promotes percolation in systems of disks where the conservative interactions are purely repulsive.This initial response is readily understood in terms of the effective attraction induced by the activity, which increases the connectivity of the particles, leading to divergence of the linear cluster size at lower number density than in the underlying passive system.With our modified iterative Boltzmann inversion scheme, we have extracted the detailed form of an effective pair potential that reproduces the structure of the active fluid as far as two-body correlations, while showing that higher-order correlations are also present.
Less expected is the nonmonotonic effect of activity, which is seen as a minimum in the percolation threshold around Pe = 10 and suppression of percolation relative to that point upon further increasing the activity, despite the continually increasing mean coordination number.This effect can be related to a study of the phase behaviour of ABPs with LJ attraction by Redner, Baskaran and Hagan (RBH) 73 .In that work, strong clustering and even phase separation occur in the absence of activity at sufficiently low temperature due to the conservative LJ forces.Introducing moderate activity breaks up the clusters by driving the particles out of their potential energy minima, reversing the effects of the LJ attraction and yielding a single active phase.Stronger activity then leads to phase separation by a different mechanism, MIPS, just as it does even for hard ABPs.RBH's paper concentrates on this re-entrance of the phase behaviour.In our work, which is confined to the pre-MIPS regime, the ABPs lack conservative attraction.Hence, unlike the system studied by RBH, there is no phase separation or attraction-enhanced clustering at zero activity.Without conservative attraction to obscure it, we therefore initially see the effects of activity-induced attraction by increased clustering and percolation at small Pe.However, analogous to the observations on phase separation by RBH, 73 increased activity reverses the effects of attraction, breaking up the clusters and suppressing percolation.The key difference for the hard ABPs in our work is that the initial attraction is itself induced by weak activity before the clustering effects that it induces are partially destroyed by stronger activity.Like the LJ system studied by RBH, a MIPS transition arises for ABPs at even higher activity 3 but we have deliberately confined our study of percolation to the one-phase, pre-MIPS regime.
Effective pairwise potentials can be deduced for ABPs up to Pe = 10, enabling the radial distribution functions of the active system to be reproduced in passive MC simulations.However, comparison of the active and mapped-passive simulations shows that higher-order correlations are present even for very weak activity.These many-body effects are effectively invisible at the twobody level of description in the sense that the radial distribution functions are identical to within the statistical noise, which itself is minimal.Our use of IBI to extract effective pairwise potentials from active simulations is just one example of the broader challenge of deducing interactions from trajectories, 74,75 which may originate from simulation or experiment, in or out of equilibrium, and may involve directional or many-body contributions in addition to isotropic pairwise terms.Work towards a more general and sophisticated treatment of these interactions is currently underway.

Fig. 1 (
Fig. 1 (a) Percolation probability R(L, ρ) versus density ρ for Péclet numbers in the range 0 ≤ Pe ≤ 40 for simulations of ABP in a box of linear size L = 200σ .(b) Percolation probability as a function of density for weak activity, Pe = 1 at four different linear box sizes.The four curves intersect at the same density, corresponding to the percolation threshold.The inset shows the results of scaling the curves onto the L = 200 curve using the known scaling exponent for passive systems in two dimensions ν = 4/3.

Fig. 2 (
Fig. 2 (a) The value of density ρ for a given percolation probability p = 0.2, 0.3 . . .0.8 versus L −1/ν for activity Pe = 1 and ν = 4/3.The lines are best fits through points with the same value of p.The intercept of the lines with the vertical axis gives ρ c for infinite box size and the spread of the intercepts gives the uncertainty in the threshold.(b) The percolation threshold as a function of Péclet number.Most error bars are smaller than the symbols.The inset of shows the monotonically increasing average coordination number with Pe at the percolation threshold.

Fig. 3
Fig. 3 Snapshots of the L = 200σ system at a density of ρ = 0.540σ −2 and (a) Pe = 1, (b) Pe = 10, (c) Pe = 40.The five largest clusters are coloured red (largest of all), magenta, green, cyan and orange.Black particles belong to clusters that are not in the five largest.Only panel (b) contains a percolating cluster.
, and machine-learning approaches by Bag and Mandal 66 and by Ruiz-1-8 J o u r n a l N a me , [ y e a r ] , [ v o l .] ,

Fig. 4 (
Fig. 4 (a) The radial distribution function for active dynamics simulations at Pe = 2, ρ = 0.54σ −2 and the result of passive Monte Carlo simulations with an effective potential after 35 iterations of IBI.(b) potentials produced by IBI at different values of the Péclet number as labelled.(c) Percolation probability curves for active systems (filled symbols) and passive systems with effective potentials from IBI (open symbols), system size L = 200σ .Each colour represents a different activity level.The Monte Carlo result for passive PHS (i.e., Pe = 0) is included for reference.