Effective interactions in active Brownian particles

Abstract

We report an approach to obtain effective pair potentials which describe the structure of two-dimensional systems of active Brownian particles. The pair potential is found by an inverse method, which matches the radial distribution function found from two different schemes. The inverse method, previously demonstrated via simulated equilibrium configurations of passive particles, has now been applied to a suspension of active particles. Interestingly, although active particles are inherently not in equilibrium, we still obtain effective interaction potentials which accurately describe the structure of the active system. Treating these effective potentials as if they were those of equilibrium systems, furthermore allows us to measure effective chemical potentials and pressures. Both the passive interactions and active motion of the active Brownian particles contribute to their effective interaction potentials.

---

1 Introduction

Dispersions of active particles are non-equilibrium systems: they dissipate energy in powering an active component to the particle velocities.^1,2 The collective behaviour of active particles manifests both structural and dynamic effects, as observed in, e.g., flocks of birds, schools of fish, and suspensions of bacteria.³ Inspired by natural systems, synthetic active particle systems can be designed to follow a particular behaviour, through tuning the chemistry or physics of colloidal particles and/or their external fields.⁴ Modelling active particles as active Brownian particles (ABPs) provides a simple means to better understand their behaviour; these are Brownian particles with an added self-propulsion, the direction of which diffuses over time.^5,6

Structural effects observed in suspensions of active particles include flocking,⁷ clustering,^8,9 chain and ring formation,¹⁰ and motility-induced phase separation (MIPS), which occurs in active Brownian particles due to particles slowing down in higher density regions.^6,11–13 The physics of active Brownian particles comprises passive interactions and the particles’ active motion.¹⁴ So far, it has not proved straightforward to find a thermodynamic framework for active matter. For example, it appears that mechanical pressure is only a state variable for active particles in certain systems.^13,15,16

In the framework for ABPs, the passive interactions are described by the pair potential u(r), while the active component is controlled by parameters such as the self-propulsion force and translational and rotational diffusion coefficients. A convenient dimensionless measure of activity is the Péclet number, which compares propulsion to diffusive motion.¹⁷ The pair potential u(r) specifies the direct interaction between two particles at separation r. Together with temperature and density, depending on the form of u(r)—e.g., hard-core repulsion,¹⁸ soft repulsion¹⁹ or Lennard-Jones type interaction²⁰—it sets the equilibrium structure that would be obtained in the absence of activity.

A steady radial distribution function, g(r), can be obtained from ABPs given the particle number density, pair interaction, u(r), and the strength of self-propulsion relative to diffusion. In this way, ABPs can form a steady state structure with a reproducible g(r), in a way similar to equilibrium systems.²¹ In this work, we aim to address whether an equilibrium approach could help us to understand the non-equilibrium structure of active particles. We seek to do this by finding effective pair potentials of an equivalent passive system. We begin by noting that a variety of definitions for effective potentials have been defined in previous works, alongside the warning of Louis,²² that—in systems with density-dependent interactions—there is no single effective pair potential; the correct effective potential needs to be derived according to its desired use.

Perhaps the simplest definition of an effective pair potential is the potential of mean force, w(r).²³ Physically, the potential of mean force, w(r), represents the effective interaction felt between two particles due to both their direct pair interaction and the averaged influence of all other particles in the system. It is related to the radial distribution function through g(r) = −exp[βw(r)]. However, only in the dilute limit does an equilibrium simulation using w(r) fully reproduce the observed g(r). Microscopy experiments on active particles have been conducted with very dilute dispersions, in which case the potential of mean force is an appropriate choice of effective pair potential.²⁴

Our approach seeks effective pair potentials which, when simulated as a passive system (with only these effective pairwise interactions), can regenerate the structure of the ABPs, as described by a g(r)—and not just in the dilute limit. This contrasts with, for example, the approach of Turci and Wilding,²⁵ who also aimed to recreate the structure of the active systems. They did this by separately finding 2, 3, … N-body interactions of ABPs by simulating systems containing 2, 3, … N particles, respectively. Their effective N-body interactions then do not depend on density. Our approach instead finds only an effective pair potential u_eff(r) which alone regenerates g(r). Since ABPs do have effective higher-order interactions, this u_eff(r) accounts for the effect of these. We will therefore find our u_eff(r) to be density-dependent. In this text, u_eff(r) is used to refer specifically to our definition of effective pair potential.

Due to the complexity of the relationship between g(r) and u(r) for a passive system—related by a simulation, or approximated using perturbation theory²⁶—we do not expect there to be an analytical relationship between an effective potential u_eff(r) and the parameters (density, u(r), and strength of self-propulsion relative to diffusion) input into simulations of active particles. We note that a machine learning approach already exists to find activity parameters (one-body force) and u(r) (two-body force), tested with ABPs and a variety of potentials.²⁷ The choice of whether to find u_eff(r) or the true interaction u(r) and activity parameters, depends on whether the primary interest is to understand the system's structure or dynamics.

Effective pair potentials, u_eff(r), have been used to predict phase separation of active particles; these effective potentials have been obtained both by iterative Boltzmann inversion (IBI)²⁸ and by deriving a relationship between u(r) and u_eff(r) from the low-density limit.²⁹ For weakly active systems, effective pair potentials have been obtained using an adaptation of IBI: the passive potential is used as input, and only the active contribution to u_eff(r) is found via IBI.²¹ We note that the definition of u_eff(r) used by Evans et al.²¹ is the same as ours. Iterative Ornstein–Zernike inversion, a modification of IBI, has been used to obtain effective interactions between granular particles—another inherently non-equilibrium system—again to describe phase separation and segregation, as in ref. 30.

Inverse methods have been established for equilibrium suspensions of purely passive particles, to relate their structure—as described by the radial distribution function g(r)—to the underlying pair potential u(r). Such methods, based on molecular simulation, include IBI³¹ and inverse Monte Carlo.³² They involve an initial guess for u(r) followed by an update scheme. Other methods are based on machine learning.^33,34 A particularly efficient form of inversion is based on test-particle insertion (TPI), with the method's efficiency resulting from removing the need to rerun a Monte Carlo simulation at every iteration.³³ In addition, not regenerating particle coordinates with further Monte Carlo simulation aids convergence to the correct potential.³⁵ As with IBI, TPI is a model-free inverse method, not requiring any assumptions about the form of u(r). The main limitations of TPI, also experienced with IBI, are difficulties at high density and requiring equilibrium.³³ The test-particle insertion method involves finding the potential energy change associated with inserting hypothetical test-particles. With TPI, the difficulty at high density arises due to poor sampling of ‘successful’ insertions.³⁶

This paper applies TPI to relate the structure of active particle suspensions to effective pair potentials. The inversion raises the question if there is an analogue of Henderson's (1974)³⁷ theorem: as long as the active particles have reached a steady-state structure, described by g(r), should we expect a unique corresponding u_eff(r) at a given temperature and density?

In cases where u_eff(r) can indeed be found, we subsequently measure thermodynamic properties—such as the pressure and chemical potential—of a system of passive particles interacting via u_eff(r), which we hope will shine light on questions around the construction of a statistical mechanical framework describing active matter.³⁸

In what follows, we will outline the theory of the inverse method, the simulations of ABPs, and the potentials modelled in Section 2. The results of the inverse method as applied to ABPs are presented and discussed in Section 3, before drawing conclusions.

2 Theory

2.1 Inverse method

The radial distribution function, g(r), is the ratio of the local number density about a reference particle, ρ(r), and the bulk number density, ρ:

g_DH(r) = ρ(r)/ρ, (1)

where r is the position relative to the reference particle. Typically, g(r) is found via the distance-histogram method (denoted by the subscript DH), whereby the domain is divided into concentric rings of width Δr about the reference particle, and the local number density in each ring evaluated.

An alternative method is based on test-particle insertion, where g(r) is the ratio of the local to bulk ensemble average:

Here, we use the subscript TPI to denote that g(r) is found by the test-particle insertion method. Also, β = 1/k_BT, with k_BT being the thermal energy. In the case of an equilibrium system, Ψ is the additional potential energy from the hypothetical insertion of a test particle.^35,39 Only considering pairwise interactions between particles,

where r_ti denotes the distance between the test particle t and a real particle i, and N is the number of particles in the system. In practice, we only consider particle–test-particle distances less than a cutoff r_c beyond which it is assumed that βu(r) ∼ 0. In the numerator of eqn (2), Ψ(r) is measured with respect to a specific particle–test-particle distance r; in the denominator, Ψ is no longer position dependent. Note that, in our notation, ρ and Ψ refer to the reference density and potential, respectively, whereas ρ(r) and Ψ(r) refer to distance-dependent ones.

In this work on out-of-equilibrium ABP systems, instead of using the true thermodynamic potential Ψ, we adopt an effective additional potential energy, Ψ_eff. In eqn (2), we replace Ψ with Ψ_eff. We calculate our effective only considering effective pairwise interactions between particles.

For an equilibrium system, we note from eqn (2) that g_TPI(r) is a function of βΨ(r), which is in turn a function of βu(r). This leads to the formulation of an inverse method to find the βu(r) which makes the g(r) obtained from each method match, i.e., g_TPI(r) = g_DH(r). From Henderson's uniqueness theorem, just one such u(r) up to a constant should be obtained for a given structure, at a given temperature and density.³⁷ Extending this to steady-state ABP systems, we likewise find the βu_eff(r) that makes g(r) from the two methods match.

Here, we use the following approach: we begin an iterative predictor-corrector scheme with an initial guess,

u_eff,0(r) = −k_BT ln [g_DH(r)], (4)

where the index 0 denotes the 0th iteration. Then, the corrector proposed by Schommers (1973),⁴⁰is used to update the guess of u_eff,j(r) until convergence is achieved, where g_TPI,j(r) is found by test-particle insertion and j is the iteration number.

2.1.1 Convergence. Convergence is assessed via the metric , where r has been discretised into bins with midpoints r_i. The initial guess [eqn (4)] corresponds to u₋₁(r) = 0 inputted into the right-hand side of eqn (5); this is the potential of mean force, which we denote as w(r).

Monte Carlo (MC) simulations (details in Appendix A) of passive particles with the obtained βu_eff(r) are then run, and the obtained g_DH,MC(r) is compared with the g_DH(r) from the original simulations of the active particles. Convergence of the inverse method, and agreement between g_DH(r) and g_DH,MC(r), were checked for all examples in this paper, but are not shown for brevity.

The potential of mean force is expected to be a good description of dilute passive systems. For relatively dilute passive systems, this leads to fast convergence of the inverse method. Similarity between u(r) and w(r) for an LJ passive system, at the same density, is shown in Appendix C.

For a dilute gas, the pair potential is equal to the potential of mean force, hence g(r) = −exp[βu(r)]. For denser fluids, this can be considered the first term of an expansion of g(r) in terms of number density.⁴¹ We write g(r) = −exp[βw(r)], which is g(r) = −exp[β(u(r) + u_ind(r))], where we have split w(r) into:

w(r) = u(r) + u_ind(r), (6)

with u_ind(r) accounting for indirect interaction. We interpret some of the results for the effective pair potentials in this way, i.e., w_eff(r) = u_eff(r) + u_ind,eff(r).

2.2 Simulations

2.2.1 Potentials. The motion of ABPs combines passive interactions with an active component to their velocity. For the passive component, we model three different forms of interaction potential between particles: the Lennard-Jones (LJ), the Weeks–Chandler–Anderson (WCA) and the shoulder potential.

The LJ potential combines a steep short-range repulsion with an attractive tail;²⁰ we truncate and shift the potential at r = 2.5σ:

where r is the centre-to-centre distance between particles, and ε is the potential well depth, which we set as the energy scale. The WCA potential is a purely hard-core repulsive potential formed by truncating and shifting the Lennard-Jones potential to zero after a cutoff of r = 2^1/6σ:⁴²

The shoulder potential has two length scales: a hard repulsive core (diameter σ) surrounded by a soft repulsive shell (diameter σ_s, height ε_s), truncated and shifted at r = 2.8σ:

The stiffness of the core is described by n and the steepness of the shoulder decay is described by k₀. This work takes n = 14 and k₀ = 10/σ.⁴³ All the examples in this work use σ_s = 2.5σ and ε_s = 2ε.

2.2.2 Active Brownian particles. Trajectories of the coordinates of ABPs were obtained via 2D large-scale atomic/molecular massively parallel simulator (LAMMPS) simulation.⁴⁴ The periodic box size was adjusted to contain 2500 particles. The time evolution of the position, r_i, and orientation relative to the x-axis, θ_i, of the ith active particle in an overdamped system is described by:andwhere D_t and D_r are the translational and rotational diffusion coefficients, respectively. We assume that D_t and D_r can be varied independently in an active system.¹³ Activity arises via a self-propulsion force, F_a, acting along the orientation vector n_i, at angle θ_i to the x-axis. The components of the thermal forces, ξ_i and ξ_i,θ, are white noise with zero mean and correlations 〈ξ^α_i(t)ξ^β_j(t′)〉 = δ_ijδ_αβδ(t − t′), where α and β are the x- and y-components, and 〈ξ_i,θ(t)ξ_j,θ(t′)〉 = δ_ijδ(t − t′), with t denoting time.

We vary the inputted u(r) [eqn (10)], density and Péclet number, Pe = 3F_aM/σD_r, where M = D_t/k_BT is the mobility. The Péclet number gives the ratio of advective active transport to diffusive transport.¹³ In this work, Pe is changed by varying F_a, with fixed M and D_r. We obtain snapshots of the generated particle coordinates for each combination. It is ensured that the simulation has run for sufficient time to guarantee that the system is in steady state. We average over the steady-state snapshots and obtain g_DH(r). This is inverted according to Section 2.1 to obtain βu_eff(r) for each set of input parameters.

The (reduced) number density is computed as ρ = σ²N/A, where N = 2500 is the number of particles contained within the box area A. We simulated runs of length 10⁷ time steps, with a time step of Δt/τ = 10⁻⁵, where , with m being the mass of a particle. Every 10 000th frame is saved for analysis (corresponding to 1000 saved frames). Up to the first 200 of these saved frames are discarded, due to not yet being at steady state.

2.3 Calculations

2.3.1 Effective chemical potential. For a system to be in thermodynamic equilibrium, besides thermal and mechanical equilibrium, it is also required to be in chemical equilibrium. This is obtained when the chemical potential, μ, of each species is constant throughout the system. For ABPs, it is known that there is no unique definition of state variables, such as pressure and temperature,^16,45 and similarly for chemical potential.⁴⁶ However, test-particle insertion leads to a simple method to calculate an effective chemical potential, μ_eff:³⁹

μ_eff = μ° + k_BT ln ρ + μ_ex, (12)

where we assume that there exists a reference chemical potential, μ°. The excess contribution is given by

μ_ex = −k_BT ln(〈exp[−βΨ_eff]〉). (13)

This effective chemical potential is the chemical potential of a passive system with pairwise interactions described by βu_eff(r). We note that the ensemble average 〈exp[−βΨ_eff]〉 was already calculated for the denominator of eqn (2) (replacing Ψ with Ψ_eff in the equation), and so we directly obtain βμ_ex from βu_eff(r). The ensemble average is calculated as if we had an equilibrium system, with particles interacting via u_eff(r). Since the box size was kept the same throughout a given LAMMPS simulation, we follow the approach for finding a canonical ensemble average. We note that there exist alternative approaches to defining chemical potentials for active systems, constructed such that the gradients in the chemical potential correspond to flows.^6,17,46,47 In contrast, our definition of effective chemical potential is related to the system structure (not dynamics).

2.3.2 Pressure.
2.3.2.1 Effective pressure. The pressure of a passive system described by u_eff(r) can be calculated. This approach does not give the real pressure of the system (see Section 2.3.2.2), since it does not account for the active contribution to the pressure.⁴⁸

The traditional method for calculating a pressure from a pair potential would be through the virial equation.⁴⁹ However, this would require calculating the derivative du_eff(r)/dr. Since the inversion only obtains binned values for u_eff(r), the calculated derivative will be noisy, leading to questions about the best smoothing and uncertainty in the calculated effective pressure, P_eff.

We adopt an alternative method that does not involve calculating the virial, and so avoids the disadvantages associated with evaluating the derivative du_eff(r)/dr. The test-volume method requires evaluating the energy changes with hypothetical changes in the system volume, whereby the distances between particles are proportionally expanded or contracted.^50,51 In 2D, we seek the limit as ΔA → 0 of

where ΔU is the difference in potential energy between the system at areas A + ΔA and A. The potential U is computed as

The ensemble average 〈exp(−βΔU)〉 is found by averaging over the snapshots. Note that we use the particle coordinates as generated from the ABP simulations in this calculation—as opposed to using a passive simulation with u_eff(r). This is because the effective pressure is found by probing how the effective free energy changes with small changes in volume. This probes how u_eff(r) changes with small changes in r. However, we find a discretised u_eff(r), which has two consequences: (1) making P_eff sensitive to the discretisation (so we will present error bars for P_eff based on this), and (2) making the coordinates from a passive simulation with u_eff(r) not as accurate as the original coordinates. For the latter reason, we use the original ABP coordinates, to obtain as accurate a measure of P_eff as we can.

As further justification for finding the effective chemical potential as per our definition, the relationship between P_eff and μ_eff is discussed in light of the Gibbs–Duhem equation in Appendix D. Evaluating μ_eff therefore offers another route to finding P_eff, which might be less sensitive to discretisation. Moreover, we note that we find μ_eff using the original ABP coordinates, since our inverse method naturally gives us Ψ_eff based on these coordinates. It is therefore consistent to also use these original ABP coordinates for the calculation of P_eff.

2.3.2.2 Total pressure. The pressure of an active system can be used to probe its phase behaviour.^48,52 Whilst there is some debate as to how to compute the pressure, P, for active systems,¹⁶ we adoptwhere here the force between particles i and j, F_ij, acts along r_ij, such that F_ij·r_ij = −du(r)/dr|_rijr_ij.⁵³ The first term is the ideal contribution. The second term is the contribution from activity, termed the ‘swimming pressure’.^54,55 The final term is the contribution from the interaction potential (i.e., the usual virial term). Our main aim is to acquire a measurement of pressure for the ABPs, in order to compare with our effective pressures. We adopt the definition that we have used in previous work.⁵³

3 Results and discussion

3.1 Effective pair potentials

A series of systematic studies are carried out, varying (1) the underlying passive pair potential (for the shoulder, WCA and LJ potentials), (2) Péclet number (for the WCA potential), and (3) density (for the shoulder potential). The particle number density was set to ρ = 0.3 in all cases, except for the study varying density. This was chosen to be sufficiently dense to show some structure, and not too dense as to pose difficulty for test-particle insertion. Example snapshots from both the LAMMPS (i.e., molecular dynamics, MD) simulations and the MC simulations run using the obtained βu_eff(r) are presented for interest in Fig. 1 (a subset) and in Appendix G (full set). Visually good agreement is obtained in the structure of the MD and MC snapshots, as expected, since g_DH(r) has been regenerated with low error. Differences in higher-order structural correlations are unlikely to be visible by eye.

Fig. 1 Zoomed snapshots of the MD (top row) and MC (bottom row) simulations for the studied potentials for selected densities and Péclet numbers. From left to right: From the shoulder potential study at varying ρ, we present ρ = 0.4 (at Pe = 120, see Fig. 5). Then, for the study varying potential type, we present the shoulder, WCA and LJ potentials (all at ρ = 0.3 and Pe = 300, see Fig. 2). (Snapshots for the complete set of studied parameters are shown in Appendix G.)

3.1.1 Varying potential types. We first observe that the inverse method successfully finds βu_eff(r) with different potential types: LJ, WCA and shoulder potentials. For each potential, Fig. 2 compares the passive u(r) used in the LAMMPS simulations [eqn (10)] with βu_eff(r) obtained by the inverse method.

Fig. 2 Varying potential type, with fixed Pe = 300 & ρ = 0.3, for (a) LJ, (b) WCA and (c) shoulder potentials. Each plot shows the passive potential u(r) and βu_eff(r) from the inverse method.

Qualitatively similar shapes are obtained for u_eff(r) for the WCA potential as for the LJ potential. In part, this is because at Pe = 300, the activity is much more important than the underlying passive interactions.⁵⁶ However, this similarity is interesting, since WCA u(r) (Fig. 2b) is entirely repulsive, whilst LJ u(r) (Fig. 2a), in addition to its short-range repulsion, includes an attractive tail. We conclude that the activity is responsible for the attractive well in u_eff(r) for the WCA potential, resulting in similarity with the LJ potential's u_eff(r). This attractive well in u_eff(r) was also found by ref. 21, which used a pseudo-hard sphere u(r). Evans et al.²¹ find that weak activity promotes percolation, even before MIPS takes place. For both LJ and WCA, we observe small jumps in u_eff(r) around r = 2σ. From corresponding plots of g(r) (see Appendix B), which have a secondary peak around r = 2σ, it appears that these bumps are merely from the structure.⁵⁷ Activity must be responsible for the slight long-range repulsion in u_eff(r) around r = 2.2σ for the LJ system.²¹ The presence of this jump in the shoulder potential's u_eff(r) might be hidden by the two length scales in u(r).

Interestingly, whilst the shoulder potential is entirely repulsive, u_eff(r) shows two attractive wells, around r = σ and r = σ_s. These effective attractive regions must be due to the activity—an emergent phenomenon of repulsive ABPs.¹³ For passive systems, multi-well potentials have been found to give rich assembly results in simulations, in both 2D and 3D.^58,59 We have previously mapped a phase diagram for 2D shoulder potential ABPs.⁵³ The phase diagram has MIPS and homogeneous regions. Interestingly, within each phase, a rich range of structures is observed, including monomers, dimers, trimers, chains and hexagonal structures. This richness arises from the shoulder potential having two length scales—which manifests as two attractive wells in u_eff(r).

Appendix B demonstrates the convergence of the inverse method for each potential type, additionally showing the very good agreement between g_DH(r) from the original LAMMPS simulation and g_DH,MC(r) from MC simulations of passive particles with the obtained βu_eff(r). This is a key result: it shows that βu_eff(r) is able to regenerate the structure, and so, in this respect, is truly an effective passive potential for the active Brownian particles. It is an interesting result in itself that there is a βu_eff(r) that regenerates g(r) in an equilibrium Monte Carlo simulation. Furthermore, we consider what thermodynamic properties, such as the pressure, might be encompassed in βu_eff(r).⁴⁸ This is explored in Sections 3.2 and 3.3.

3.1.2 Varying the Péclet number. The effect of varying Péclet number is also studied: this shows the influence of varying activity.

For the shoulder potential, at Péclet numbers lower than 120, particles do not have enough energy to penetrate the soft repulsive shell, and so they tend to behave as particles with larger diameter (σ_s = 2.5σ). The effective packing fraction based on this larger diameter is greater than one, leading to ‘frustrated’ structures.⁵³ This results in particles forming snake-like paths of high local density. For the LJ potential, at Péclet numbers lower than 300, there was too much particle clustering—i.e., again, areas of local high density. This was reflected in peaks in g(r) of high magnitude (e.g., >10). These areas of high local density mean that test-particle insertion does not work well. This is due to the difficulty of ‘successfully’ inserting a particle in a region that is already very dense.⁶⁰ This leads to either convergence not being achieved, or, if achieved, a noisy βu_eff(r).

We therefore choose to present the effect of the Péclet number using the WCA potential (see Fig. 3), which does not have the above problems and allows the study of a range of Pe. Note that there is no universal density limit above which test-particle insertion fails; it depends on the effective potential. In Appendix E, we discuss the second virial coefficient and what that might tell us about the phase behaviour of the ABPs.

Fig. 3 Varying Péclet number for the WCA potential, with fixed ρ = 0.3: (a) Plot of βu_eff(r) at different Pe. (b) Plot of the magnitude of βu_eff(r)'s well depth as a function of Pe.

Plots of βu_eff(r) for the WCA potential with a selection of Pe numbers (Pe = 6, 30, 60, 120 & 300) at ρ = 0.3 are shown in Fig. 3a. For the WCA potential at each Pe tested, excellent agreement between g_DH(r) and g_TPI(r) is obtained (not shown). We learn that the inverse method is able to work for different Pe. We additionally plot the magnitude of the well depths in u_eff(r) as a function of further Pe values in Fig. 3b.

The magnitude of the well depths in u_eff(r) are observed to initially increase as Pe is increased (from Pe = 6 to Pe = 22). However, at higher Pe (Pe > 22), the magnitude of the well depths in u_eff(r) decreases as Pe is increased. Whilst the apparent attraction in u_eff(r) is attributed to the activity, increasing activity appears to decrease the particle clustering. This is because ABPs tend to slow at high density for steric reasons, but increased activity reduces the effect of this slowing down.^21,61 In the limit of Pe → 0, we would recover the u(r) for a purely passive simulation of the shoulder potential. In Appendix F, we present how the mean free time between particle collisions varies with Pe, finding similar trends to those observed with u_eff(r).

3.1.3 Effective potential of mean force. Taking the case of ABPs with a LJ potential shown in Fig. 2a, the effective potential of mean force, w_eff(r) [eqn (4)], is plotted in Fig. 4a, as an example to compare with u_eff(r). The potentials w_eff(r) and u_eff(r) appear very similar; this contributes to the rapid convergence of the inverse scheme when using w_eff(r) as an initial guess.⁶² In Fig. 4b, we compare the target g_DH(r) found from the original data, with that from an MC simulation using w_eff(r). The similarity between these g(r) reflects that between u_eff(r) and w_eff(r). The first peak in the target g_DH(r) is well reproduced. The second peak is higher in g_DH,MC(r) (generated using w_eff(r)) than in g_DH(r), as a result of the second well in w_eff(r). A third peak appears in g_DH,MC(r) which is absent in g_DH(r), as higher-order structural correlations are more sensitive to small differences in the interaction potential.

Fig. 4 Demonstration of the closeness of βw_eff(r) to βu_eff(r), using the example of ABPs with a LJ potential, with Pe = 300 & ρ = 0.3: (a) Comparison of βw_eff(r) with βu_eff(r). (b) Comparison of g_DH(r) from the original data with g_DH,MC(r) from a MC simulation using βw_eff(r). These g(r) should not agree perfectly, as they correspond to different u(r), though they are similar due to the closeness of βw_eff(r) to βu_eff(r).

Therefore, it seems to be the case for relatively dilute active Brownian particles, that the effective potential of mean force explains most of the structure. The example in Fig. 4 is at high Pe (Pe = 300); i.e., activity dominates the particles’ motion. We learn that, for active Brownian particles at ρ = 0.3, effective direct interactions, which are predominantly a result of the activity, contribute most to the effective pair potential. Interestingly, activity indirectly leads to structure, but emerges as a direct contribution to the effective pair potential.

3.1.4 Varying the density. Finally, the effect of varying particle density is examined, using the shoulder potential as an example, at Pe = 120.

The shoulder potential is an interesting example because of its two length scales. The results are presented in Fig. 5. This is a particularly relevant study, as neither the u(r) or activity change in the LAMMPS simulations of the ABPs. Nevertheless, we observe differences in u_eff(r) purely due to change in density. As the density is increased, the ratio of the peak in g(r) around r = σ to that around r = σ_s increases. This suggests that, at high density, more particles are forced within the soft-shell region due to steric restrictions. The density-dependence of g(r) is unsurprising, since this is also observed in passive systems, due to structure increasing as density increases.⁶³ However, in passive systems, the u(r) obtained from the inversion method clearly does not depend on the density—since this u(r) should have converged to the original u(r) used in the Monte Carlo simulation that generated the particle coordinates. In contrast, the u_eff(r) obtained for the ABPs do depend on the density, since the emergent effects of activity depend on the density. This is via density-dependent particle velocities: whilst density-dependence of the particle velocities is not directly coded into the model, collisions are expected to slow particles down at high density.^11,64 We see that the activity leads to non-pairwise additivity.

Fig. 5 Shoulder potential at varying density, with fixed Pe = 120. (a) Plot comparing g_DH(r) and g_TPI(r), based on the converged βu_eff(r) from the inverse method. (b) Plot comparing βu_eff(r) with βw_eff(r). The colour code in the key of (a) also applies to (b).

Looking at u_eff(r), at all densities, attractive wells appear at the two characteristic distances. The depth of the well around r = σ increases slightly with increasing density. More generally, considering the region between the two wells, the apparent attractiveness of the interactions increases with increasing density. This corresponds with the increased penetration of the soft-shell observed in g(r) with increasing density. The depth of the second attractive well, around r = σ_s, is less affected by the density.

Fig. 5b additionally shows w_eff(r) at each density. This makes it clear that similarity between u_eff(r) and w_eff(r), as for passive systems, decreases with increasing density. At ρ = 0.049, there is very good agreement between u_eff(r) and w_eff(r). As the density increases, the agreement worsens because w_eff(r) changes more with increasing ρ than u_eff(r) does. This reflects our understanding that effective indirect interactions u_ind,eff(r) increase with increasing density.

3.2 Effective chemical potential

Following the calculation method given in Section 2.3.1, we plot either β(μ_eff − μ°) = ln ρ + βμ_ex (Fig. 6a, varying ρ) or βμ_ex (Fig. 6b, fixed ρ) for the ABPs.

Fig. 6 Plots of effective chemical potential as a function of (a) ρ (shoulder potential with Pe = 120) and (b) Pe (shoulder, WCA and LJ potentials, ρ = 0.3).

Considering the shoulder potential at varying density, Fig. 6a shows that the ideal chemical potential increases with density, whilst the effective excess chemical potential decreases; the overall effective chemical potential shows a maximum.

At fixed density, ln ρ is constant, so Fig. 6b plots the effective excess chemical potential only. For the WCA potential, there is a non-monotonic trend in the effective chemical potential as Pe is varied; this likely reflects the non-monotonic behaviour in the well depths in u_eff(r) (Fig. 3b). In addition to the WCA potential, Fig. 6b also plots results for the LJ and shoulder potentials. The effective chemical potentials for the LJ, WCA and shoulder potentials at Pe = 300 are close to each other. This is expected: the passive pair potentials become less important as the activity is increased, and the systems with the different potentials become increasingly similar to each other. For the shoulder potential, the effective chemical potential is higher at Pe = 300 than at Pe = 120; higher chemical potential corresponds with a less negative test-particle insertion energy [when averaged as in eqn (13)].

3.3 Pressure

3.3.1 Effective pressure. In this section, using eqn (14), we find the pressures of passive systems with the same u_eff(r) as the ABPs (using the ABP coordinates, as explained in Section 2.3.2). The results are presented in Fig. 7a and c. Strictly speaking, u_eff(r) is a function of ρ. Whilst this is ignored in calculating P_eff, the result should be reasonably accurate, since the density-dependence of u_eff(r) is not too large (see Fig. 5b).²² As with the effective chemical potential, the effective pressure appears similar between the LJ, WCA and shoulder potentials at high Pe (Pe = 300). For the shoulder potential at varying densities (Fig. 7a), the effective pressure shows a maximum, as did the effective chemical potential (Fig. 6a). Error bars are drawn in Fig. 7a and c based on the standard deviation of results from different bin widths. Due to the large error bar, the sign of P_eff at ρ = 0.4 is uncertain. We stress that there is still less variance with bin width using this method than with the virial method.

Fig. 7 Effective pressure with error bars [(a) and (c)], and total pressure, divided into its components of ‘ideal + interaction’ and ‘swimming’ pressures [(b) and (d)]; both as a function of ρ [shoulder potential with Pe = 120, (a) and (b)] and Pe [shoulder, WCA and LJ potentials, ρ = 0.3, (c) and (d)]. The key in (b) also applies to subfigure (a), and, likewise, the key in (d) also applies to subfigure (c).

For a passive system, the Gibbs–Duhem equation at constant temperature can be written as dμ/dρ = (1/ρ)dP/dρ.^47,65 By comparing the data shown in Fig. 6a and b with those in Fig. 7a and c, Appendix D attempts to test the Gibbs–Duhem equation, both for varying ρ (shoulder potential) and Pe (WCA potential). However, we observe in Fig. 3 and 5 that u_eff(r) varies with ρ and Pe. Louis²² explains the challenges of deriving different thermodynamic properties, such as μ_eff and P_eff, from density-dependent effective pair potentials. Whilst we expect eqn (13) for μ_eff to still be exact, due to its similarity with the method used to find u_eff(r) in the first place, we noted above that this is not the case with our calculation of P_eff. In future work, the role of many-body effects should be considered, in determining whether the Gibbs–Duhem relationship holds for the effective thermodynamics of ABPs. Nevertheless, u_eff(r) does not vary too much with ρ or Pe (see Fig. 3 and 5), which likely explains the resemblance between Fig. 6a and 7a, and between Fig. 6b and 7c. We note that alternative flow-based approaches for defining chemical potentials for active systems obey the Gibbs–Duhem equation under local density approximations.^46,47

3.3.2 Total pressure. We compare these effective pressures with the total pressure, P, of the active systems. Note that we substitute k_BT/ε = 0.01 into eqn (16), as this was used in these LAMMPS simulations—this only affects the ideal gas pressure term, which we find to be small compared to the other terms, in the presence of significant activity.

In Fig. 7b and d, we plot the total pressure and its constituents: the swimming pressure and the passive pressure (the sum of the ideal and interaction terms). The general trend is for the swimming and (to a lesser extent) passive pressures (and hence the total pressure) to increase with increasing ρ and Pe. Consequently, at higher ρ and Pe values, the total pressure is order of magnitudes larger than the effective pressure. It is expected from its coefficient that the swimming pressure should increase with Pe (Fig. 7d), in addition to ṙ_i also depending on F_a [eqn (10)]. The summation over the N particles explains why the swimming pressure increases with ρ (Fig. 7b). The passive pressure likely increases with density due to particles being forced closer together (Fig. 5a), such that more particles lie at distances with large −du(r)/dr values. For the WCA potential, the passive pressure increases with Pe (Fig. 7d). Correspondingly, it can be seen in Fig. 3a that the position of the minimum in u_eff(r) shifts to smaller r values as Pe increases. The passive pressure for the shoulder potential simulations was similar at both Pe = 120 and Pe = 300.

Conclusions

We have newly demonstrated that, for ABPs, a steady state g(r) can be inverted to obtain βu_eff(r). This works for different potential types, Péclet numbers and densities tested. The u_eff(r) combines the underlying u(r) and activity. The influence of these factors on u_eff(r) was demonstrated for a variety of potentials, with systematic changes in potential types, Pe and density. The inverse method converges quickly because, at low density, u_eff(r) resembles w_eff(r), which is used as the initial guess. Through u_eff(r), we see the non-additive effects of activity. The purely repulsive shoulder and WCA potentials give particularly insightful u_eff(r), since their u_eff(r) contain attractive wells: this is a clear demonstration of the emergent attractive behaviour of ABPs. Furthermore, the shape, such as the depth, of these attractive wells offers a means to quantify this emergent behaviour. Varying density, whilst keeping the u(r) and activity parameters used in the LAMMPS simulation constant, produces density-dependent u_eff(r) that offer a way to quantify density-dependent emergent activity effects.

We therefore believe that this work presents another way to characterise and understand the behaviour of active particle suspensions—in particular, their structure. Whilst u_eff(r) could be viewed as a proxy for g(r), this work raises the question, is there any more meaning to u_eff(r)? Here, we have used u_eff(r) to calculate μ_eff and P_eff, which may help to address thermodynamic questions about active systems. Further work could explore if u_eff(r) could give any meaningful insight into the particle dynamics, and any quantitative information on the activity parameters, given the u(r) used in the numerical simulation. It has been demonstrated with passive systems that the test-particle insertion inverse method can be extended to many-body interactions, i.e., finding uⁿ(r), where n denotes the n-body interaction.⁶⁶ This could also be explored with ABPs, finding .²¹ We expect these higher-order interaction terms to be non-zero, since we have found u_eff(r) to be density-dependent.

Another route for future exploration, currently underway, is using our inverse method with mixtures of active and passive particles. Effective interactions between active particles in active–passive mixtures have been previously found in the dilute limit,²⁴ whilst our method offers an approach applicable also beyond this limit, and to find effective interactions between all combinations of particle types. In addition, our method could be applied to other classes of active matter besides ABPs, such as run-and-tumble particles.⁶⁷ We expect the method to work for any active system that obtains steady-state structures and does not have areas of high local density. This includes systems with non-reciprocal interactions. It is worth mentioning that, in real systems, active particles may experience non-reciprocal interactions, either originating from hydrodynamic interactions⁶⁸ or differences in surface mobility towards chemical gradients.⁶⁹ Our test-particle insertion method ignores any potential non-reciprocity in the interactions; it just finds an effective (reciprocal) pair potential from the structure. We still expect our method to be able to find a u_eff(r), and thus also a μ_eff and P_eff. These might show how asymmetries in effective interactions shift collective steady states.

Author contributions

Clare R. Rees-Zimmerman: conceptualization, data curation, formal analysis, investigation, methodology, project administration, software, validation, visualization, writing – original draft, writing – review & editing. C. Miguel Barriuso Gutierrez: data curation, formal analysis, investigation, methodology, project administration, software, validation, writing – review & editing. Chantal Valeriani: conceptualization, funding acquisition, investigation, methodology, project administration, resources, supervision, writing – review & editing. Dirk G. A. L. Aarts: conceptualization, formal analysis, funding acquisition, investigation, methodology, project administration, resources, supervision, validation, writing – review & editing.

Conflicts of interest

There are no conflicts to declare.

Data availability

The code for the test-particle insertion method can be found at https://github.com/creeszimmerman/TPI.

Supplementary information (SI) is available. See DOI: https://doi.org/10.1039/d5sm00706b.

Acknowledgements

C. R. R. Z.'s work is funded by a Junior Research Fellowship from Christ Church, University of Oxford. The authors would like to thank José Martín-Roca (Universidad Complutense de Madrid) for checking the active pressure calculations. C. V. acknowledges the funding awards IHRC22/00002 and PID2022-140407NB-C21 from MCIN/AEI/10.13039/501100011033 and FEDER, EU. The authors thank Francisco Alarcón (University of Guanajuato) for useful discussions.

References

1. U. Khadka, V. Holubec, H. Yang and F. Cichos, Nat. Commun., 2018, 9, 3864.
2. N. R. Smith, Phys. Rev. E, 2023, 108, L022602.
3. R. S. Negi, R. G. Winkler and G. Gompper, Soft Matter, 2022, 18, 6167–6178.
4. C. W. Shields and O. D. Velev, Chem, 2017, 3, 539–559.
5. A. Zöttl and H. Stark, Annu. Rev. Condens. Matter Phys., 2023, 14, 109–127.
6. J. Stenhammar, D. Marenduzzo, R. J. Allen and M. E. Cates, Soft Matter, 2014, 10, 1489–1499.
7. L. Caprini and H. Löwen, Phys. Rev. Lett., 2023, 130, 148202.
8. R. M. Navarro and S. M. Fielding, Soft Matter, 2015, 11, 7525–7546.
9. B. M. Mognetti, A. S. Šarić, S. Angioletti-Uberti, A. Cacciuto, C. Valeriani and D. Frenkel, Phys. Rev. Lett., 2013, 111, 245702.
10. V. Telezki and S. Klumpp, Soft Matter, 2020, 16, 10537–10547.
11. M. E. Cates and J. Tailleur, Annu. Rev. Condens. Matter Phys., 2015, 6, 219–244.
12. S. A. Mallory, C. Valeriani and A. Cacciuto, Annu. Rev. Phys. Chem., 2018, 69, 59–79.
13. J. Martin-Roca, R. Martinez, L. C. Alexander, A. L. Diez, D. G. A. L. Aarts, F. Alarcon, J. Ramírez and C. Valeriani, J. Chem. Phys., 2021, 154, 164901.
14. A. Mauleon-Amieva, M. Mosayebi, J. E. Hallett, F. Turci, T. B. Liverpool, J. S. van Duijneveldt and C. P. Royall, Phys. Rev. E, 2020, 102, 032609.
15. G. Junot, G. Briand, R. Ledesma-Alonso and O. Dauchot, Phys. Rev. Lett., 2017, 119, 028002.
16. A. P. Solon, J. Stenhammar, R. Wittkowski, M. Kardar, Y. Kafri, M. E. Cates and J. Tailleur, Phys. Rev. Lett., 2015, 114, 198301.
17. J. Stenhammar, A. Tiribocchi, R. J. Allen, D. Marenduzzo and M. E. Cates, Phys. Rev. Lett., 2013, 111, 145702.
18. J. K. Percus and G. J. Yevick, Phys. Rev., 1958, 110, 1–13.
19. F. H. Stillinger and T. A. Weber, J. Chem. Phys., 1978, 68, 3837–3844.
20. J. E. Lennard-Jones, Proc. Phys. Soc., 1931, 43, 461.
21. D. Evans, J. Martín-Roca, N. J. Harmer, C. Valeriani and M. A. Miller, Soft Matter, 2024, 20, 7484–7492.
22. A. A. Louis, J. Phys.: Condens. Matter, 2002, 14, 9187.
23. J. G. Kirkwood, J. Chem. Phys., 1935, 3, 300–313.
24. Y. Mu, L. Lei, J. Zheng, W. Duan, Z. Wang, J. Tang, Y. Gao and Y. Wang, ACS Nano, 2022, 16, 6801–6812.
25. F. Turci and N. B. Wilding, Phys. Rev. Lett., 2021, 126, 038002.
26. M. Gottschalk, AIP Adv., 2021, 11, 045026.
27. M. Ruiz-Garcia, C. M. Barriuso G., L. C. Alexander, D. G. A. L. Aarts, L. M. Ghiringhelli and C. Valeriani, Phys. Rev. E, 2024, 109, 064611.
28. B. Trefz, S. K. Das, S. A. Egorov, P. Virnau and K. Binder, J. Chem. Phys., 2016, 144, 144902.
29. T. F. F. Farage, P. Krinninger and J. M. Brader, Phys. Rev. E:Stat., Nonlinear, Soft Matter Phys., 2015, 91, 042310.
30. G. M. Rodríguez-Liñán and M. Heinen, J. Comput. Phys., 2019, 394, 232–242.
31. D. Reith, M. Pütz and F. Müller-Plathe, J. Comput. Chem., 2003, 24, 1624–1636.
32. A. Lyubartsev, Coarse-grained modeling of biomolecules, CRC Press, Boca Raton, FL, 2017, ch. 1.
33. C. R. Rees-Zimmerman, J. Martín-Roca, D. Evans, M. A. Miller, D. G. A. L. Aarts and C. Valeriani, J. Chem. Phys., 2025, 162, 074103.
34. S. M. Kampa, F. Sammüller, M. Schmidt and R. Evans, Phys. Rev. Lett., 2025, 134, 107301.
35. A. E. Stones, R. P. A. Dullens and D. G. A. L. Aarts, Phys. Rev. Lett., 2019, 123, 098002.
36. C. R. Rees-Zimmerman, A. Heafield, D. Ellerbeck, A. E. Stones, R. P. A. Dullens and D. G. A. L. Aarts, JCIS Open, 2025, 20, 100156.
37. R. Henderson, Phys. Lett. A, 1974, 49, 197–198.
38. E. Fodor and M. C. Marchetti, Phys. A, 2018, 504, 106–120.
39. B. Widom, J. Chem. Phys., 1963, 39, 2808–2812.
40. W. Schommers, Phys. Lett. A, 1973, 43, 157–158.
41. R. Hunter and L. White, Foundations of Colloid Science, Oxford University Press, Oxford, 2nd edn, 2001.
42. J. D. Weeks, D. Chandler and H. C. Andersen, J. Chem. Phys., 1971, 54, 5237–5247.
43. N. V. Gribova, Y. D. Fomin, D. Frenkel and V. N. Ryzhov, Phys. Rev. E:Stat., Nonlinear, Soft Matter Phys., 2009, 79, 051202.
44. S. Plimpton, J. Comput. Phys., 1995, 117, 1–19.
45. A. P. Solon, Y. Fily, A. Baskaran, M. E. Cates, Y. Kafri, M. Kardar and J. Tailleur, Nat. Phys., 2015, 11, 673–678.
46. S. Paliwal, J. Rodenburg, R. van Roij and M. Dijkstra, New J. Phys., 2018, 20, 015003.
47. S. Hermann, D. de las Heras and M. Schmidt, Phys. Rev. Lett., 2019, 123, 268002.
48. M. Sanoria, R. Chelakkot and A. Nandi, Phys. Rev. E, 2021, 103, 052605.
49. M. Allen and D. Tildesley, Computer simulation of liquids, Clarendon, Oxford, 1983.
50. R. Eppenga and D. Frenkel, Mol. Phys., 1984, 52, 1303–1334.
51. V. I. Harismiadis, J. Vorholz and A. Z. Panagiotopoulos, J. Chem. Phys., 1996, 105, 8469–8470.
52. R. G. Winkler, A. Wysocki and G. Gompper, Soft Matter, 2015, 11, 6680–6691.
53. J. Martín-Roca, R. Martinez, F. Martínez-Pedrero, J. Ramírez and C. Valeriani, J. Chem. Phys., 2022, 156, 164502.
54. S. C. Takatori, W. Yan and J. F. Brady, Phys. Rev. Lett., 2014, 113, 028103.
55. A. K. Omar, Z.-G. Wang and J. F. Brady, Phys. Rev. E, 2020, 101, 012604.
56. E. Chacón, F. Alarcón, J. Ramírez, P. Tarazona and C. Valeriani, Soft Matter, 2022, 18, 2646–2653.
57. R. Jeanneret and D. Bartolo, Nat. Commun., 2014, 5, 3474.
58. M. Engel and H.-R. Trebin, Phys. Rev. Lett., 2007, 98, 225505.
59. J. Dshemuchadse, P. F. Damasceno, C. L. Phillips, M. Engel and S. C. Glotzer, Proc. Natl. Acad. Sci. U. S. A., 2021, 118, e2024034118.
60. G. C. Boulougouris, I. G. Economou and D. N. Theodorou, Mol. Phys., 1999, 96, 905–913.
61. J. Su, M. Feng, Y. Du, H. Jiang and Z. Hou, Commun. Phys., 2023, 6, 58.
62. T. C. Moore, C. R. Iacovella and C. McCabe, J. Chem. Phys., 2014, 140, 224104.
63. A. L. Thorneywork, R. Roth, D. G. A. L. Aarts and R. P. A. Dullens, J. Chem. Phys., 2014, 140, 161106.
64. S. E. Moran, P. W. A. Schönhöfer and S. C. Glotzer, New J. Phys., 2022, 24, 063007.
65. J.-L. Barrat and J.-P. Hansen, Basic concepts for simple and complex liquids, Cambridge University Press, Cambridge, 2003.
66. A. E. Stones and D. G. A. L. Aarts, J. Chem. Phys., 2023, 159, 194502.
67. I. Santra, U. Basu and S. Sabhapandit, J. Stat. Mech.: Theory Exp., 2020, 2020, 113206.
68. E. S. Bililign, F. B. Usabiaga, Y. A. Ganan, A. Poncet, V. Soni, S. Magkiriadou, M. J. Shelley, D. Bartolo and W. T. M. Irvine, Nat. Phys., 2022, 18, 212–218.
69. C. H. Meredith, P. G. Moerman, J. Groenewold, Y.-J. Chiu, W. K. Kegel, A. van Blaaderen and L. D. Zarzar, Nat. Chem., 2020, 12, 1136–1142.

---