# Phase Behavior of Active Swimmers in Depletants: Molecular Dynamics and Integral Equation Theory

## Abstract

We study the structure and phase behavior of a binary mixture where one of the components is self-propelling in nature. The inter-particle interactions in the system were taken from the Asakura-Oosawa model, for colloid-polymer mixtures, for which the phase diagram is known. In the current model version the colloid particles were made active using the Vicsek model for self-propelling particles. The resultant active system was studied by molecular dynamics methods and integral equation theory. Both methods produce results consistent with each other and demonstrate that the Vicsek model based activity facilitates phase separation, thus broadening the coexistence region.

###### pacs:

29.25.Bx. 41.75.-i, 41.75.LxVarious phenomena involving systems containing active particles have been of significant recent research interest (1); (2); (3); (4); (5); (6); (7); (8); (9); (10); (11); (12); (13); (14); (15); (16); (17); (18); (19). Simple examples are flocking of birds (1); (2), dynamics in a bacterial colony (11), etc. Self-propelling character of active species make such systems extremely complex. While the literature in this area gained significant volume in recent time, many basic questions related to both equilibrium and nonequilibrium statistical mechanics remain open. Examples (20) are phase behavior and criticality, various fluctuation relations, kinetics of phase transitions, etc.

Considering the current status of the subject, realistic modeling involving shape, flexibility, etc., of the constituents, is very difficult. Definition of temperature in such systems is questionable and is a subject of much interest (9); (10). These systems are nonequilibrium in nature and so, to what extent methodologies of equilibrium statistical mechanics are applicable, e.g., to study phase behavior, remains to be seen. Thus, it is extremely challenging to study active matter from theoretical and computational points of view.

Interesting experiments (15); (18); (19) reported that addition of active particles in a system can dramatically alter its phase behavior. With that objective, in this letter, we present results for the phase behavior from molecular dynamics (21) (MD) and integral equation theoretical (22) (IET) study of a model where an active species swims in an environment of passive particles.

Interactions among active particles can be of various types (15); (18). In this work we are interested in constructing a model with an effective attraction among self-propelling particles, which exhibits phase separation in the equilibrium limit. Thus, our objective is to study the influence of dynamic clustering in active particles, on the phase behavior of the corresponding passive system. Our results, in conjunction with Ref. (15), demonstrate that, based on the type of activity, phase separation can either be facilitated or suppressed. In the Vicsek type model, the tendency of active particles to form clusters enhances phase separation.

In our model, the standard interactions among particles are given via a variant of the well-known Asakura-Oosawa (AO) model (23); (24) for colloid () and polymer () mixtures. In this work, we make the colloids self-propelling motivated by the fact that bacteria or other active objects in solutions are in the size range of colloidal particles in colloidal dispersion. Further, in many circumstances (e.g., living objects) there are also depletants in the solution which may create attractions among colloidal particles. It is then an intriguing question, what is the consequence of the interplay between such “ordinary” interactions and the fact that one deals with active, rather than passive particles. Moreover, colloid-polymer mixtures are popular for the study of phase transitions in condensed matter because of the easy observability of interfacial phenomena on the single particle scale. While the original AO model (23) cannot be straightforwardly generalized (it is only defined in thermal equilibrium due to the ideal gas character of the penetrable soft spheres representing the random walk-like polymers), we use a variant (24) with nontrivial potentials. Interparticle interactions in this model (24) are described as (, being the locations of th and th particles)

(1) |

where or , while

(2) |

with , , , , and . Similar size ratios are conceivable in colloidal dispersions with biopolymers as depletants (25). Equilibrium phase behavior of this entirely passive model exhibiting entropy driven phase separation was estimated with good accuracy (24). As stated above, in the present work, we introduce activity, via the well-known Vicsek model (1); (2) as described below, to observe how the phase behavior of the AO model changes.

The Vicsek model was introduced to study self-propelled dynamics in interacting biological systems where the interaction of a particle with its neighbors is decided by the average direction of motion of the neighbors. We have used as a cut-off radius to define neighborhood. At each instant of time, an active particle was supplied with an acceleration , in addition to the one due to interatomic interactions Eqs. (1,2), in the direction of the average velocity of its neighbors. We study this model by MD as well as IET calculations.

When activity is introduced, eventually the temperature assumes a value higher than the one initially assigned. This is due to additional velocity imparted according to Vicsek rule. A reason for choosing the AO model as the passive limit is its weak sensitivity to temperature. In real systems the solvent carries off all the heat produced by active particles. This effect, in our study has been avoided by an appropriate choice of the model. Furthermore, in the AO model the polymer reservoir packing fraction, which regulates the depletant density, plays the role that inverse temperature would play in a standard system where competition of energy and entropy controls the phase behavior.

In our MD simulations, we have implemented a Langevin thermostat to control the overall temperature. There we have essentially solved the equation

(3) |

particle mass being set to be same for all particles, is a damping coefficient, is the interparticle potential, is a -correlated noise and the temperature. Eq. (3) was solved via implementation of the velocity Verlet algorithm (21) with time step in units of . Henceforth all lengths are measured in units of and energy in units of . Further, for convenience, we set , , and Boltzmann’s constant to unity. To introduce activity we have used , unless mentioned otherwise. Throughout the work the starting temperature was with . For reasons mentioned earlier, this temperature does not correspond to the effective temperature in the active case. The thermostat has the effect to maintain a steady state with an enhanced temperature of the colloids after a rather short transient. Thus no artificial temperature rescaling is needed to compensate for continuous energy pumping due to the Vicsek model. Also, we have checked using other values of , viz., and , that our results on the phase diagram are independent of these choices within statistical errors.

An IET approach (22) was adopted that uses standard Ornstein-Zernike (22) (OZ) equations relating direct pair correlations with the total correlation functions. It is well-known that asymmetric models, as the one used here pose significant difficulty in the IET approach in obtaining closure relations needed to supplement the OZ equations. To overcome this difficulty, we took a pragmatic approach and empirically determined the best combination of closures by comparing structural and thermodynamic output from the IET with the MD data at a reference state. We found that the modified hypernetted chain (22) closure for AA and the Percus-Yevick (22) closure for AB as well as BB work best and so results are presented using these closures.

In order to apply IET methods to the active system, we assume that it can be mapped onto an effective passive one, interparticle interactions chosen in such a way as to reproduce the structural properties of the active system. If this assumption holds true, then both structure and phase behavior of the active system can be studied by applying the IET formalism outlined above to the passive system, onto which the active system has been mapped. In order to test our assumption, we selected a single point on the phase plane [that corresponds to Fig.1(a)] and “inverted” the simulated radial distribution function of the active system to obtain an effective interaction potential for the passive system mimicking the active one ( and are the same as in the MD simulations). For the passive case all interactions are given by Eqs.(1,2). We then performed IET calculations for the effective passive system over the entire phase plane. The validity of our assumptions was confirmed by good agreement between IET and MD results for the phase diagram of the active system.

Via MD simulations and IET, in this letter we compare the results for structure and phase behavior of the active model with that of the corresponding passive one (24). All our MD results were obtained using periodic boundary conditions starting from random initial configurations. Both MD and IET results confirm that the miscibility gap of the active model widens.

In Fig. 1 (a) we compare the structure factors of the active and passive models, calculated as

(4) |

for . In Eq. (4), , are respectively the number of particles of species , and is the total number of particles. We estimate the phase behavior by varying the composition of and particles. Previously (24), for the passive model, the phase behavior was obtained in vs plane, being the packing fraction of species . To be more specific and , where (, being the system volume) is the density of species . The state point in Fig. 1 (a) is () (). For the passive case (24), the critical point () is at (). The symbols in this figure show the MD results and the IET calculations are shown by lines. The agreement between MD and IET calculations must be appreciated. Interestingly, the sharp rise of at small for the active model indicates the presence of long wavelength fluctuations not present in the passive system.

Fig. 1 (b) presents a comparison between the interaction potential of the passive system and the effective interaction potential obtained via inversion of the simulation data for the active system. Intriguingly, an additional attraction is present in the active case, in contrast to the model of Ref. (15). This provides an explanation for the widening of the coexistence region seen below.

In Fig. 2, from MD simulations, we show two snapshots from the passive and active models at a state point different from Fig. 1(a) and closer to the passive coexistence curve (see Fig. 4). The left snapshot shows the passive system, the right one is for the active model. It can clearly be seen that, as opposed to the passive one, the active system has nicely phase separated.

Next we estimate the phase diagram. To achieve this objective via MD simulations, we have started with a homogeneous mixture of and particles and waited for the system to reach a steady state which can be identified from the potential energy and density profiles (see Fig 3). If the snapshots showed phase separation in the steady state, we have calculated profiles for and as a function of . Staying away from the interfacial region, we extract the values of and for both coexisting phases (26). On the other hand, for the IET method, a divergence criterion for the structure factor in limit was used to obtain the spinodal.

Finally, our central result, the phase diagram, is presented in Fig. 4 in the - plane. The result from the passive system (24) is also added for comparison. It is clearly seen that the introduction of a Vicsek type activity facilitates phase separation. Simulation and IET results are in qualitative agreement for both the passive and the active case. Some quantitative disagreement that appears between the theory and simulation can possibly be attributed to the mean field nature of the IET. The IET result for the binodal was obtained following the procedure outlined in Ref. (27). As expected, it lies outside the spinodal, approaching it in the critical region. Comparing the active and passive cases, it is seen that the coexistence region opens up when activity is introduced. But strong finite-size effects and critical slowing down prevent us from obtaining points very close to the critical point reliably via MD simulations with moderate system sizes.

It is worth noting that different types of activity can lead to phase behavior which is not only quantitatively but also qualitatively different. In Ref. (15) activity leads to randomly enhanced mobility resulting in an additional effective repulsion. In a Vicsek type model, originally introduced to study swarming behavior, active particles tend to cluster. This leads to an additional attractive interaction which enhances phase separation.

So far our results refer only to . In Fig. 5 we show the inverse concentration-concentration structure factor , defined as ()

(5) |

for different . For small enough the OZ behavior (22)

(6) |

with and being the susceptibility and correlation length, is nicely visible. The stronger enhancement of with the increase of is suggestive of the fact that the phase gap widens which is due to the stronger effective attraction among active particles (28). This is directly demonstrated in the inset.

In conclusion, we have presented results for the structure and phase behavior of a physically motivated model mixture containing active particles. The results are compared with the corresponding passive systems. Our molecular dynamics simulations and integral equation theory (this being the first time in literature to have been applied for the study of active matter) are in reasonable agreement and show that the tendency to form clusters in the Vicsek model leads to an effective attraction among colloids and enhances phase separation. This result is in contrast with previous work in which the model of activity differs. Basically, depending upon how motility is introduced, qualitatively different trends result. While it needs to be seen whether the conclusions hold true in more general situations, our results should be stimulating for future experiments.

It would be interesting to study critical phenomena for such phase separating active systems, and to look at interfacial properties and hydrodynamic effects. Analogous to the observation in phase separating passive fluids, our preliminary study of an active system also suggests fluctuation induced broadening of interfaces with increasing system size. Despite many interesting works involving active matter, to the best of our knowledge, understanding of such important aspects is missing. We hope our work will influence experimentalists and theorists to pursue such problems.

Acknowledgment: SKD and SE acknowledge research visits to the Johannes Gutenberg-Universität in Mainz. SKD is also grateful for financial support from DFG (SPP1296, TR6/A5) and JGU, Germany, JNCASR and DST, India, as well as ICTP, Italy. SKD and PV acknowledge discussion with F. Schmid. PV acknowledges discussions with M. Oettel. P.V. and B.T. would also like to acknowledge the MAINZ Graduate School of Excellence.

### References

- T. Vicsek, A. Cziroḱ, E. Ben-Jacob, I. Cohen and O. Shochet, Phys. Rev. Lett. 75, 1226 (1995).
- A. Czirók and T. Vicsek, Physica A 281, 17 (2000).
- R. Dreyfus et al., Nature 437, 862 (2005).
- J.R. Howse et al., Phys. Rev. Lett. 99, 048102 (2007).
- J. Palacci, C. Cottin-Bizonne, C. Ybert, L. Bocquet, Phys. Rev. Lett. 105, 088304 (2010).
- H.P. Zhang, A. Béer, E.-L. Florin and H.L. Swinney, Proc. Natl. Acad. Sci. U.S.A. 107, 13626 (2010).
- H. Ke, S. Ye, R.L. Carroll and K. Showlater, J. Phys. Chem. A 114, 5462 (2010).
- A.C. Callan-Jones and F. Jülicher, New J. Phys. 13, 093027 (2011).
- S. Wang and P.G. Wolynes, J. Chem. Phys. 135, 051101 (2011).
- D. Loi, S. Mossa and L.F. Cugliandolo, Soft Matter 7, 3726 (2011).
- C. Valeriani, M. Li, J. Novosel, J. Arlt and D. Marenduzzo, Soft Matter 7, 5228 (2011).
- I. Buttinoni, G. Volpe, F. Kümmel, G. Volpe and C. Bechinger, J. Phys.: Condens. Matter 24, 284129 (2012).
- J. Bialké, T. Speck and H. Löwen, Phys. Rev. Lett. 108, 168301 (2012).
- J.-F. Joanny and S. Ramaswamy, J. Fluid Mech. 705, 46 (2012).
- J. Schwarz-Linek et al., Proc. Natl. Acad. Sci. U.S.A. 109, 4052 (2012).
- M.E. Cates, D. Marenduzzo, I. Pagonabarraga and J. Tailleur, Proc. Natl. Acad. Sci. U.S.A. 107, 11715 (2010).
- A. Kaiser, H.H. Wensink and H. Löwen, Phys. Rev. Lett. 108, 268307 (2012).
- J. Palacci, S. Sacanna, A.P. Steinberg, D.J. Pine and P.M. Chaikin, Science 339, 936 (2013).
- X.H. Sun, C. Danumah, Y. Liu and Y. Boluk, Chem. Eng. J. 198, 476 (2012).
- A. Onuki, Phase Transition Dynamics (Cambridge University Press, Cambridge, 2002).
- D. Frenkel and B. Smit, Understanding Molecular Simulations: From Algorithm to Applications (Academic Press, San Diego, 2002).
- J.P. Hansen and I.R. McDonald, Theory of Simple Liquids (Academic, London, 1986).
- S. Asakura and F. Oosawa, J. Chem. Phys. 22, 1255 (1954); S. Asakura and F. Oosawa, J. Polym. Sci. 33, 183 (1958); A. Vrij, Pure Appl. Chem. 48, 471 (1976); R.L.C. Vink, J. Horbach and K. Binder, Phys. Rev. E 71, 011401 (2005).
- J. Zausch, P. Virnau, K. Binder, J. Horbach and R.L.C. Vink, J. Chem. Phys. 130, 064906 (2009).
- E.A.G. Jamie et al., J. Phys.: Condens. Matter 20, 404231 (2008).
- H. Watanabe, N. Ito and C.-K. Hu, J. Chem. Phys. 136, 204102 (2012).
- J. Dzubiella, C. N. Likos and H. Löwen, J. Chem. Phys. 116, 9518 (2002).
- L. Angelani et al., Phys. Rev. Lett. 107, 138302 (2011).