Solvent-induced depletion interactions in multiparticle collision dynamic simulations

Molecular-dynamics-coupled multiparticle collision dynamic (MPC-MD) simulations have emerged to be an e±cient and versatile tool in the description of mesoscale colloidal dynamics. However, the compressibility of the coarse-grained °uid leads to this method being prone to spurious depletion interactions that may dominate the colloidal dynamics. In this paper, we review the existing methodology to deal with these interactions, establish and report depletion measurements, and present a method to avoid arti¯cial depletion in mesoscale simulation methods.


Introduction
Computer simulations have largely shown to be an essential approach to unravel the intricate behavior of many soft matter and biophysical systems, both in equilibrium and nonequilibrium conditions. In many of these systems, the behavior of the main system components, such as colloids,¯laments, or membranes, is crucially a®ected by the role of the solvent in which they are suspended. Most relevant solvent-induced e®ects are hydrodynamic interactions, depletion, screening of electrostatic interactions, or phoretic e®ects, whose strengths and dependencies can be very diverse. In order to provide an e®ective description of these solvent e®ects, and given the large separation in length and time scales, various mesoscopic simulation methods have been profusely developed in the last decades, process which is still ongoing. 1 Some of the most extensively used methods are lattice Boltzmann (LB), 2-4 dissipative particle dynamics (DPD), [5][6][7] or multiparticle collision (MPC) dynamics. [8][9][10] The establishment of these methods runs logically in parallel with their applications to di®erent types of systems, with which the methods are constantly tested, such that both advantages and pitfalls are being detected and new improved versions or application details are being developed. Of particular interest in this work is the presence of solvent-induced depletion interactions within the MPC method.
Depletion interactions are well known in soft matter physics, as they are present in various types of colloid mixtures. Depletion refers to an e®ective force between colloids, most often attractive, induced by steric interactions with other particles, typically smaller colloidal particles or polymers. 11,12 Within this context, the smallest particles are commonly referred to as the depletants. The steric interactions of the smaller with the larger colloidal particles de¯ne a depletion layer around the colloids with thickness related to the depletant size. This layer is not accessible to depletants in case that two colloids get close enough for their depletion layers to overlap. Con¯gurations with close colloids are then energetically more favorable for the depletants as the volume freely accessible to them increases. Simultaneously, the osmotic pressure on the colloids imbalances when their depletion layers overlap, since they are no longer surrounded by a symmetric concentration of depletants, resulting this in an e®ective attractive interactions between colloids. The seminal model to describe depletion is the Asakura-Oosawa (AO) model, 13 which relates the attraction depletion potential U depl to the overlap volume of the two colloids depletion layers V ov , Here, refers to the number density of depletants, k B is the Boltzmann constant, T the average temperature and d the intercolloid separation. The overlap volume V ov ðdÞ can be calculated by the di®erence in volume available to the depletants for the two large spheres at in¯nite separation and at separation distance d at which the depletion layers overlap, Therefore, the size of the small particles crucially determines the range and the strength of depletion interactions. In this way, for standard dispersions, where the solvent particles are orders of magnitude smaller than the suspended colloids, the attraction should be completely negligible. Nevertheless, in the context of mesoscopic simulations, signi¯cant solvent-induced depletion interactions might appear due to the use of coarse-grained solvents, which bring additional length scales, not so negligible in comparison with the colloid size, and that appear in the form of lattice grids, or interaction ranges. These induced attractions can then be considered as an artifact of the chosen algorithms. This problem has been identi¯ed in various previous occasions, 14,15 and it is typically circumvented by including an additional separation and/or repulsion between colloids. Being correct, this solution has potential disadvantages, since these additional interactions are not implemented on the basis of real physical grounds; and other e®ects, such as short range hydrodynamic interactions, can be also nonphysically screened away.
In this work, we focus on the MPC method with suspended colloids modeled with the molecular dynamics coupling. We quantify the strength, and the range of the depletion interactions with various choices of potentials, and summarize the existing approaches to avoid such simulation artifacts. We then present a strategy in which solvent-induced depletion interactions can be safely neglected, without incurring in any additional unphysical consequence. Finally, as an application example, the pair dynamics of thermophoretic colloids is discussed.

Method
The MPC method describes the°uid as a collection of point particles of mass m which perform alternating streaming and collision steps. During streaming, particles propel ballistically, this is without any type of collisions. In regular time intervals h, the collision step is performed. Therein, particles are binned into cubic collision cells of side length a, and interchange momentum with the particles in the same collision box, such that the total amount is conserved before and after the collision. Here, we focus on the stochastic rotation collision rule, which also conserves kinetic energy, providing a description of the solvent within the microcanonical ensemble. With this collision rule, each particle rotates its velocity relative to the center-of-mass velocity of the collision cell by an angle a around a random axis. 16 Other collision rules 17,18 would not change any of the conclusions here presented.
MPC o®ers a particle-based, coarse-grained, mesoscale description of a°uid, resolving both hydrodynamics and heat transport correctly. This°uid has, in the described formulation of the method, the equation of state of an ideal gas, due to the lack of potential interactions between solvent particles. This also means that MPC particles are point like, this is they do not have any speci¯c size. 10 The implementation of the complex structures such as colloids of polymers, varies depending on the problem. 19,20 In a large number of applications in polymer physics, 21-23 the interaction between polymer monomers and MPC solvent particles is performed in the MPC collision step, where the monomers are just considered as MPC particles with di®erent mass. This approach has repeatedly proved that hydrodynamic interactions are properly considered, 24,25 and the lack of excluded volume interactions between monomers and solvent particles makes that solvent-induced depletion interactions are absolutely not present. In some other applications, the excluded volume interactions between monomers and solvent particles is though desirable. One possible solution is to include colloid-°uid interactions by hard-sphere collision rules (bounceback or specular re°exions) for MPC particles reaching the surface of a colloid. 16,26 This coupling does practically not su®er from arti¯cial depletion forces, but it also considers practically no speci¯city of the surface collision what is in some cases a disadvantage. The most standard strategy considered in MPC has so far been to include solute-solvent potential interactions which are then solved by molecular dynamics, 27-29 such that we refer to this approach as MPC-MD. The use of explicit pair potentials leads to the in principle point-like particles of the MPC°uid being assigned a certain e®ective diameter comparable to that of the colloids. This means that the MPC particles act e®ectively as depletants. Since MPC particles do only interact with each other in a coarse-grained way and can overlap to an arbitrary extent, their character as a depletant is similar to the AO model, but with a soft interaction UðrÞ in between colloid and depletant. The depletion potential can then be calculated from Eq. (1) and 14 with ¼ 1=k B T and jr 1 À r 2 j ¼ d. This expression can be integrated numerically for most potentials as will be shown later.
The most commonly employed potentials in MPC simulations are of the Mie type, 30 UðrÞ ¼ 4"½ðs=rÞ 2n À ðs=rÞ n þ C, for r < r c and zero otherwise. Here, r is the pairwise distance, s the typical colloid radius and " and n describe the strength and the steepness of the potential. Purely repulsive interactions are obtained with C ¼ and r c ¼ 2 1=n ; and attractive interactions with C ¼ 0 and an adequate cut-o® radius which for n ¼ 6 is r c ¼ 2:5, this is the Lennard-Jones potential. 31 The parameters chosen in the simulations here presented are, unless otherwise speci¯ed, s ¼ 4, ¼ 1, n ¼ 24 for attractive potentials, and n ¼ 3 for repulsive. The MPC solvent units are determined by a ¼ m ¼ k B T ¼ 1, and the chosen parameters for these simulations are very standard in MPC applications, namely collision time h ¼ 0:1, a rotation angle ¼ 120 , and an average number of particles per collision cell ¼ 10, what corresponds to a°uid with kinematic viscosity ¼ 0:79, Schmidt number Sc ¼ 13, and Prandtl number Pr ¼ 5:3.

Solvent-Induced Depletion
In order to characterize the existing depletion forces, two colloids are¯xed at various separation distances d. The bath of MPC particles act as depletants, which allows us to directly quantify the resulting colloid-colloid depletion force as a function of distance between the colloids. Integration then yields the depletion interaction potential. This is the standard method employed for MD simulations 12 which can be also applied for MPC-MD, just using the MPC particles as depletants.
Typical depletion potentials and forces are shown in Fig. 1 for various interactions. In Figs. 1(a) and 1(b), forces and potentials are displayed for the two main types of relevant interactions, namely a soft repulsive, and a steep attractive potential, which are calculated by simulations and by the analytical approach in Eq. (3), which shows to completely capture the simulation results. Since the MPC°u id is known to behave with the ideal gas equation of state, Eq. (3) should be valid for any arbitrary choice of the MPC solvent. We con¯rm this statement by comparing simulation results at di®erent rotation angles as shown in Fig. 1(c). Interestingly, we¯nd that no di®erence for the depletion forces is found even in the case of a pure ideal gas, which is simulated by switching o® all solvent-solvent interactions setting to zero. This proves that depletion interactions are solely determined by the equation of state of the mesoscopic solvent together with its interactions with the colloids, or any other complex structure treated by MD. Notably, in an ideal-gas-like solvent such as the MPC°uid, the depletion layer is exactly as long-ranged as the colloid-°uid interaction potential.
Though the dependencies of the depletion forces on the simulation parameters can be extracted from Eq. (3), it is of interest to study them on representative examples. We show the dependence on the colloid-°uid interaction strength " in Fig. 1(d) and on the°uid density in Fig. 1(e). In both cases, the point where depletion forces vanish is identical, as determined by the range of the interaction potential. However, the functional shape of the force naturally varies with the employed parameters,  which is interesting to consider at the moment of performing an educated parameter choice. It can be seen for example, that depletion forces are stronger and already signi¯cant at smaller overlaps for higher values of " and .

Existing approaches to avoid MPC-induced depletion
The depletion forces due to an MPC-MD°uid were¯rst estimated by Padding and Louis 14 for colloids interacting with repulsive interaction with the solvent particles. In order to avoid depletion, two strategies were proposed. One of these strategies was the use of an explicit counter-potential based on AO theory, which would precisely counteract the depletion forces. Although this provided a satisfactory agreement to Browninan dynamics (BD) simulations, problems were already addressed since this approach was clearly not accounting for higher-order interactions, anisotropies in the solvent density induced by external¯elds, or the time dependence in the build-up of depletion forces. 14 The most extended strategy to avoid depletion has been the use of a separation between colloids, which we de¯ne as , additional to the sum of both colloidal radii. A large enough value of will not only avoid depletion forces, but also ensure the presence of a layer of°uid in between colloids, which can also resolve lubrication forces. 14,32 In principle, in the case of a solvent with an ideal-gas-like equation of state, like the discussed version of MPC, the depletion layer extends as much as the potential interaction layer, such that an additional separation of the size of the cut-o® radius should be su±cient, ' 2r c . Simulations of colloids described by MPC-MD with this choice of 14 compared the radial distribution functions, to simulation results obtained by BD, and found signs of a remaining depletion interaction. This could be most likely attributed the soft nature of the colloid-solvent interactions and to the role of other parameters as displayed in Figs. 1(d)-1(e). In a type of combination of the two previously discussed strategies, Thakur and Kapral 33,34 choose a relatively small additional separation distance, together with an inter-colloidal repulsion. A full phase diagram of depleted bounded states was obtained, showing that, for a high enough choice of inter-colloidal repulsion, no depleted states were observed, even when the depletion layers overlap. This occurs because the attractive depletion potential is combined in this region with the repulsive potential, which results in a compensation of both e®ects.

Proposed strategy to avoid solvent-induced depletion
The previous approaches were applied to colloids with repulsive colloid-solvent interactions, and for most applications they can be considered satisfactory, although it should be considered that the resulting interaction at short distances is di±cult to evaluate. More problematic is the approach to avoid solvent-induced attraction in the case of colloids with attractive colloid-solvent interactions where the depletion layers can be considerably larger as well as the depletion interactions. 33,35 A possible solution would be the use of a large inter-colloidal separation, which is not always desirable, since it can signi¯cantly reduce all°uid-mediated colloid-colloid interactions at short ranges, such as the frequently relevant hydrodynamic interactions, or those due to solvent intrinsic gradients such as thermophoresis, or di®usiophoresis.
The main idea we use here is looking for colloid-solvent potentials with relatively small depletion layers, such that the additional inter-colloidal separation does not interfere with the physically relevant°uid-mediated short range colloid-colloid interactions. With this spirit, we explore the e®ect of an additional hard core interaction, or displacement Á combined with the main potential as already done by Kihara 36 UðrÞ ¼ 1; r Á; which reduces to the previously employed Mie potential when Á ¼ 0. Depletion is only dependent on the colloid-°uid potential width, which in this potential is r c , while other solvent-induced e®ects such as hydrodynamic interactions show to be proportional to the colloid radius, which here is then s ¼ þ Á. Therefore, one can try to introduce the shift Á while keeping the colloid radius s constant. This will minimize the changes in the functional form of the interaction, and will certainly reduce any e®ects related to the range of the colloid-°uid interaction. In Fig. 2, the depletion forces are displayed for colloids¯xed at di®erent distances, which are explicitly measured for colloids with a¯xed radius s, but with various displacements Á. For the example of a steep attractive colloid-°uid interaction in Fig. 2(a), the introduction of displacements reduces signi¯cantly the intensity of the depletion attraction, as well as its range. For the example of the repulsive colloid-°uid    Fig. 2(b), the intensity of the depletion attraction is almost not modi¯ed, opposite to its range. This means that the required colloid-colloid separation to avoid solvent-induced depletion can be quite strongly reduced in both cases. For example, with these potentials and Á ¼ 0:5s colloids with the attractive potential could use an extra separation ¼ 0:15s, or the repulsive potential with ¼ 0:10s.

Application to Thermophoretic Self-Propelled Dimers
Simulations of thermophoretic self-propelled colloids constitute a very clear example, where the correct choice of the interaction parameters is critical in order to avoid solvent-induced depletion interactions. Thermophoresis refers to the directed motion that colloidal particles experience in the presence of a steady local temperature gradient in the surrounding°uid, which can occur towards cold (thermophobic colloids) or warm areas (thermophilic colloids). 37,38 Thermophoretic propulsion of colloidal swimmers emerges naturally in MPC-MD simulations as the result of a temperature gradient and a suitable choice of colloid-°uid interaction potentials. 29,39 Typically, thermophobic behavior is obtained with steep attractive colloid-°uid interactions and thermophilic behavior with repulsive interactions. 29,40 Thermophoretic self-propelled motion can be expected in the cases of Janus or dimeric colloidal particles, where the colloid surface has an asymmetric heat capacity, what has been studied experimentally. 41,42 and by means of simulations. 22,40 The thermophoretic properties of the nonheated part produce then a propulsion against or towards the heated part. Equally sized thermophobic microdimers have shown to hydrodynamically behave like pullers, thermophobic microdimers like pushers, and half-coated Janus particles like neutral swimmers. 40 Mesoscopic simulations of microdimers clearly indicate that these hydrodynamic behaviors can be modi¯ed and even reversed for microdimers built with beads of di®erent sizes. Thermophobic microdimers with small heated beads show for example an important lateral attraction as expected for pusher type of swimmers. 35 The collective behavior of these microdimers exhibits very interesting properties. The interplay of the hydrodynamic attraction in the direction perpendicular to the propulsion direction and the thermophoretic repulsion opposite to the propulsion give rise to the self-assembly of the microdimers in planar moving structures with a well-de¯ned orientation and hexagonal order. 43 A paradigmatic type of swimmer geometry is the thermophobic dimer made out of two same-sized colloids. In this case, the intercolloidal phoretic repulsion is combined with a hydrodynamic short-ranged lateral attraction and long-ranged repulsion. This case specially highlights the necessity of choosing appropriate colloid-°uid interactions with suitable displacements using the°ow¯eld of this swimmer, which is shown in Fig. 3(a), with the corresponding quanti¯cation of the lateral°ow velocity is shown in Fig. 3(b). These results do not really depend on the particular choice of the potentials, and in particular on the choice of Á in Eq. (4). To avoid depletion in the standard nondisplaced case, Á ¼ 0, the required colloid-colloid extra separation would be ¼ 0:5s, as illustrated in Fig. 2(a). The corresponding cut on short-range hydrodynamics is illustrated by the dashed black line in Fig. 3(b), which will strongly diminish the e®ect of short-ranged hydrodynamic interactions. If the potential is then modi¯ed to that in Eq. (4) with Á ¼ 0:5s, the intercolloidal distance to avoid depletion reduces to ¼ 0:2s, which is shown by the solid black line in Fig. 3(b), which will then not a®ect the short-ranged hydrodynamic interactions. It is therefore clear that in this example, the choice of the potential can eventually screen away the short-ranged hydrodynamics by the need to avoid solvent-induced depletion, which would be physically incorrect. The collective behavior of swimmers with both choices of Á is expected then to be signi¯cantly di®erent, 35 but only simulations with a large value of Á will provide the physically correct behavior.
Another relevant aspect is how the change of the colloid-°uid interaction range, introduced when modifying the parameter Á, a®ects the results. As pointed out, the variation is quantitative, and in Fig. 3 it appears in the absolute values of the propulsion velocity v p and°ow velocities v l . This is a consequence of the thermophoretic e®ect, which depends on the colloid-°uid interaction range. 35 The e®ect that the potential displacements have on these properties is shown in Fig. 4 for an asymmetric dimer. The decrease of the propulsion velocity with increasing Á can be clearly shown in Fig. 4(a). Meanwhile, the°ow velocities in Fig. 4(b) show only to be very slightly a®ected. Overall, an intermediate choice such as Á ¼ 0:5s proves to be a very suitable choice for the potential displacement, since the corresponding ¼ 0:2s still allows to keep short range hydrodynamics, and the decrease of the related phoretic velocities is just 10% with respect to the unperturbed ones.
In order to test the previous approach, we have performed test simulations with two dimeric swimmers. In this case, where the depletion is important, bounded states had been already reported 33,35 leading to Brownian pairs or rotating pairs, with in some cases even intertwined structures. This can be quanti¯ed in systems as one shown in Fig. 5, where the distance of the phoretic colloids of the two dimers is computed for certain time. The pair displaying a constant distance corresponds to a system with an undisplaced potential where depletion leads to a bound state, which in this case even shows an spontaneous modi¯cation of the distance. On the other hand, when parameters are chosen to avoid depletion, we do not observe the formation of any bound states for any of the investigated cases. Here, we show a pair of swimmers simulated with displaced potentials which shows a strongly oscillating separation, corresponding to a free di®usive motion, which will eventually become larger. Furthermore, simulations with a large number of colliding swimmers, as those shown in Refs. 35   Besides the already discussed additional separation between colloids , the possible choices of colloid-colloid interaction potentials are limited by the requirement that they should be strongly repulsive in the range where depletion is signi¯cant. Though this requirement is often consistent with the desired physical model, a possible improvement of the method is to combine the short-ranged repulsive with another potential, such as hard-sphere bounce-back rules, in the spirit of Ref. 40, where such a combination has been employed for the colloid-°uid interaction. Colloids could then interact with arbitrary potentials, but, as soon as their depletion layers start overlapping, the bounce-back rule will ensure that they will be separated again. Furthermore, this method would also diminish the in°uence of thermal°uctuations, which with softer potentials would allow for an eventual overlap of the depletion layer. These overlaps would be impossible when the bounce-back rule is used. It should be considered that these bounce-back rule changes the position of the colloids and can thereby induce a discontinuity into the forces. For faster moving colloids, the bounce-back can be combined with steeper repulsive interaction at short ranges; and for slowly moving colloids, the disturbance will be likely not signi¯cant in any case. It is also interesting to discuss the simulation of the so-called squirmers, 44,45 where the interaction of the active colloids and the MPC solvent particles is modeled by prescribing a surface velocity at the colloids boundary. Given the absence of a well-de¯ned interaction layer, and the frequent use of virtual particles inside the colloids, depletion e®ects as described in this work, are in principle not expected to be be present in these simulations, although other artifacts such as huge density°uctuations should be carefully accounted for. 46

Summary and Conclusions
The combination of potential interactions of colloidal particles with a compressible mesoscopic solvent may induce arti¯cial depletion interactions between colloids. We have here summarized the existing approaches to avoid such simulation artifacts with the MPC-MD method and present simulation results that con¯rm the feasibility of the existing analytic approach. Furthermore, we extend the existing methodology by introducing displaced potentials to model the colloid-°uid interactions. This provides an additional way of tuning the width of the depletion layer, which is of importance when short-ranged,°uid-mediated interactions are crucial to model the physical e®ects of interest. Regarding the optimal parameters to be employed, it is important to note that minimizing the depletion layer, surface-induced e®ects, like phoresis, are also minimized, such that an educated compromise should be taken in each case. For which it is useful to remind that although various parameters have certain in°uence in the depletion interaction, the potential displacement is the one with strongest e®ect. The present results are of importance when simulating concentrated systems by the