Effective potentials induced by mixtures of patchy and hard co-solutes

The archetypal example of depletion interactions emerges when hard sphere colloids are dispersed in a solvent with non-adsorbing polymers (the so-called co-solute). In such systems, the presence of the polymers gives rise to an entropy-driven effective attraction between the colloids. More than 60 years ago, Asakura and Oosawa studied depletion interactions when an idealized polymer is considered as a co-solute.¹ In this case, the depletion interaction can be derived analytically^1,2 and the Asakura–Oosawa model has since become a widely used reference system.³ Analytic treatment of colloidal suspensions is, however, limited to a few specific cases, and the entropic depletion interaction is part of a larger family of effective interactions mediated by different kinds of co-solutes. Such effective interactions are ubiquitous in nature^4–6 and can induce phase separation, enhance crystallization and gelation, or give rise to different types of arrested states. The significant impact on the behavior of colloidal systems makes the deliberate control of effective interactions highly desirable as this would allow for the design of specific materials, which can be exploited in fields spanning catalysis, photonics, sensor and filter technology, energy storage, medicine, and food science.^7–13 

Colloidal effective interactions can be controlled not only by changing the relative size and density of the co-solute but also by choosing the proper co-solute–colloid and co-solute–co-solute interactions. For instance, the model studied by Asakura and Oosawa can be modified by using different types of colloids (e.g., hard spherocylinders^14,15 or soft spheres^16,17), different types of co-solutes (e.g., hard spherocylinders,¹⁸ star-polymers,^19,20 or self-propelling particles²¹), or different types of both (e.g., hard spherocylinders with a hard sphere co-solute,²² hard rods with a binary Lennard-Jones co-solute,²³ or microgel mixtures²⁴). Polydisperse systems have also been investigated²⁵ as well as systems in which the co-solute particles are close to a second-order critical point, thus inducing tunable long-range interactions among the colloids.^26–28 The case of specific interest for the present work is the one in which the co-solute self-assembles into a specific structure.

In investigating the use of patchy particles^29–31 as co-solutes, some of us have recently shown that if the co-solute reversibly self-assembles into non-adsorbing polymer chains, attractive oscillating two-body interactions between colloids arise.³² This is an especially intriguing result in light of several other studies, which have shown that the self-assembly of one-component systems can be controlled through designed oscillating interactions,^12,33–35 leading, for example, to the formation of specific crystals^36–40 and quasicrystals.^40–43 Yet, it is unknown how to design experimental systems to give rise to such potentials. Encouraged by the results emerging in colloidal dispersions with a self-assembling co-solute,³² in this work, we explore how the effective interaction between the colloids can be tuned by using binary co-solute mixtures, where one of the two components is capable of self-assembly. Specifically, we investigate mixtures of hard spheres (HSs) and patchy particles of equal sizes for different concentrations of the two species. In using these systems, we study whether the resulting effective potentials can be predicted from those obtained for the respective single-component co-solutes.

This paper is organized as follows: In Sec. II, we introduce the model system and the applied simulation techniques. In Sec. III A, we discuss the effective potentials mediated by the single-component co-solutes. Section III B provides an analysis of the structural properties of the co-solute mixtures in the absence of the two colloids, and in Sec. III C, we discuss the resulting effective potentials mediated by the co-solute mixtures. In Secs. III A–III C, the patchy co-solute particles considered have two attractive sites (2P). In Sec. III D, we present the same analysis, but instead of the 2P patchy co-solute, we use patchy co-solute particles with three attractive sites (3P). Finally, in Sec. IV, we summarize all results and present our conclusions.

We perform NVT Monte Carlo (MC) simulations of two colloids immersed in a binary mixture composed of hard spheres (HSs) and patchy particles of valence two (2P) or three (3P) at a packing fraction ϕ = 0.262 to evaluate the effective interactions at different species concentrations. A representative snapshot of the HS–2P co-solute mixture is shown in Fig. 1. In all cases, colloid–co-solute interaction is provided only by excluded volume and is of hard-sphere-type. The size of all different co-solute particles (HS, 2P, and 3P) is the same, and it is set to σ = 1, while the size of the two colloids is σ_c = 10σ. Consequently, σ is used as the unit of length throughout this article.

The patchy co-solute particles (2P or 3P) are modeled as hard spheres decorated with two or three attractive sites (patches), respectively. Site–site interaction is treated using the Kern–Frenkel potential,^44–46 which was originally introduced by Bol,⁴⁷ 

This potential consists of a square-well attraction,

of depth ɛ (the unit of energy) and width δ modulated by an angular function, which depends on the mutual orientation of the patchy particles. Specifically, two patches α and β on particles i and j attract each other only when the center-to-center vector $r⃗ij$ lies inside the cones of the two patches, i.e.,

where the unit vector identifying, for instance, the orientation of patch α on particle i is $n̂iα$⁠. The volume available for bonding is controlled by cos(θ_max). For all simulations, we have fixed δ = 0.119σ and cos(θ_max) = 0.894 717, ensuring the single-bond per patch condition.

After fixing the temperature T and the packing fraction ϕ, we evaluate the effective potentials V_eff for different concentrations of the two co-solutes. For most simulations, we employ a parallelepipedal simulation box with lengths L_x = 50σ and L_y = L_z = 21σ. In some cases, we have simulated a smaller number of particles using a simulation box of size 35σ × 20σ × 20σ. The box geometry always guarantees that the maximal surface-to-surface distance r of the colloids is larger than the distance at which V_eff goes to zero. Periodic boundary conditions are applied in all directions. We employ a standard umbrella sampling technique^26,27,32 to evaluate the probability P(r) of finding the two colloids at distance r, and we extract the effective interaction from the relation

where β = 1/k_BT, k_B is Boltzmann’s constant, and g(r) is the two-body radial distribution function of the colloids.

Despite the efficiency of the method, which allows us to uniformly sample all distances between the colloids, we cannot directly separate the contributions of the two co-solute species from the total effective interaction. In order to obtain this information, we calculate the effective force F_eff acting on the colloids in the following way:^48,49 at each MC simulation step, virtual displacements of size Δr are performed by the two colloids. These displacements can be positive or negative. The probability that such a displacement leads to at least one collision is

Here, Δϕ_HS is the change in the colloid–co-solute excluded volume interaction due to the virtual displacement Δr. Δϕ_HS = 0 when the virtual move results in no collision [the contribution to P_overlap(Δr) is consequently 0] and Δϕ_HS becomes infinite when the virtual move leads to at least one overlap [the contribution to P_overlap(Δr) = 1 in this case]. This enables us to write F_eff as

For a given colloid separation r, an imbalance of the overlap probability at positive Δr (⁠$limΔr→0+$⁠) and the overlap probability at negative Δr (⁠$limΔr→0−$⁠) gives rise to an effective force. If the two probabilities are equal, F_eff is zero. Notably, this way of calculating F_eff allows us to calculate the total F_eff and the individual contributions stemming from the HS and the patchy particles, respectively, by considering only overlaps of the colloids with one or the other co-solute species. In Eq. (6), this is indicated by the index i, which can be tot, HS, or 2P (or 3P, respectively).

We note that the ratio $Poverlapi(Δr)/Δr$ is insensitive to the magnitude of the displacement as long as Δr is small enough. Previous works have shown that an optimal choice of magnitude provides an average collision probability of about 5%–15%.^48,49 In this work, we aimed for a collision probability of the whole mixture of ≈10%, which implies that the collision probability for each species is lower than this value and depends on the relative amount of the species in the mixture.

Finally, we also perform MC simulations of the binary mixture of co-solutes, i.e., in the absence of the two colloids, to investigate its structural properties. To this end, we have focused on a smaller system (N = 1315) in a cubic simulation box at the same T and ϕ as used for the calculation of V_eff. From these simulations, we extract the chain length distribution P(l), i.e., the probability P of finding a chain of length l. In addition, we calculate the partial radial distribution function g_ij(r) with respect to species i and j and the concentration–concentration structure factor, which is an indicator of the demixing transition. The latter is defined as^50–52 

where x_i is the concentration of the species i and

is the partial structure factor between species i and j.

We start our discussion by reporting the effective interactions arising between two large HS colloids immersed in a single-component co-solute (either HS or 2P particles). This study is fundamental to understand and possibly predict the effective interactions mediated by mixtures of these two components. Figure 2(a) shows V_eff mediated by a co-solute of HS-type as a function of ϕ. All potentials are characterized by an attractive well for r/σ < 1, which results from the depletion of the co-solute between the colloids. In addition, on increasing ϕ, the oscillatory character of V_eff increases, signaling the emergence of correlations between the colloids due to the underlying structure of the co-solute.^53–56 

A different scenario occurs for the co-solute of 2P particles, which is shown in Fig. 2(b). Here, we focus on a low temperature value, i.e., T = 0.1, for which 2P particles are known to reversibly self-assemble into chains, whose lengths depend on the control parameters T and ϕ. In fact, the high-T region is less interesting in terms of depletion interactions since on increasing T, the bonding energy of the 2P particles becomes much smaller than the thermal energy, meaning that the particles behave essentially as HS.³² 

At low ϕ, the effective potential between the large colloids is completely attractive and, different from the standard depletion potentials, goes to zero at r ≈ 4σ. The long-range nature of the resulting potential can be understood in terms of the chains formed by the co-solute because the chains now act as depletants with an effective diameter mediating the colloid–colloid interaction. On increasing ϕ, similar oscillations as found in the effective potentials mediated by the HS co-solute appear, but the effective potentials for the 2P co-solute remain attractive at all distances. This has been explained in terms of ordering of 2P chains in the confined space between the colloids.³² In particular, it has been shown that in between the two colloids, 2P particles are correlated in orientation, reminiscent of nematization, in order to form bonds with other particles, thus satisfying both energy and entropy constraints. The resulting minima and maxima of the oscillations in V_eff are roughly found at distances of the order of the size of the depletant as in the case of the effective potential generated by the HS co-solute. We stress that the attractive oscillations shown in Fig. 2(b) are a peculiarity of the reversibly self-assembling co-solute and have not been observed in the case of effective interactions mediated by irreversibly bound polymer chains⁵⁷ since no competition between bonding and excluded volume takes place in that case.

Before analyzing the effective interactions mediated by the binary mixture of co-solutes, it is worth to investigate the co-solute behavior in the absence of large colloids. By fixing ϕ and T and mixing 2P and HS together, we obtain an extra control parameter, the 2P concentration x_2P, which determines the 2P chain length distribution P(l) and its average value ⟨l⟩, as shown in Fig. 3(a). We also note that, in the case of binary mixtures, the chain length distribution is characterized by an exponential behavior as predicted by the Flory–Stockmayer theory,⁵⁸ with an average length of the chain that increases on increasing x_2P. To gain insights into the structure of the mixture, we calculate the partial radial distribution functions with respect to all so-solute species, i.e., g_2P–2P(r) and g_HS–HS(r), and that of the mixture, g_HS–2P(r), for different concentrations. The corresponding results are reported in Fig. 3(b), revealing a striking difference between the different partial radial distribution functions. In particular, we find that the peak of g_HS–HS(r) is almost independent of concentration. The same is observed for g_HS–2P(r). On the contrary, we find that the first peak of g_2P–2P(r) decreases on increasing x_2P. To understand how the composition of the mixture affects the structure of the system, we extract the average number of neighbors using the relation

where the subscripts i and j indicate the species involved (2P–2P, HS–HS, and HS–2P) and r_min represents the distance at which the first minimum after the first peak in g_ij(r) is found. We find that the number of HS neighbors for one HS particle decreases with x_2P = 1 − x_HS since it ranges from ≈9 at x_HS = 1 (only HS particles) to ≈1 for x_HS = 0.125 at the studied ϕ. Analogously, the number of 2P neighbors for a 2P particle changes from ≈6 when x_2P = 1 to ≈2 at x_2P = 0.250. Interestingly, if we calculate the number of 2P neighbors of a 2P particle by choosing r_min = σ + δ, we find that the number of neighbors does not change much as a function of x_2P. The value is always close to 2, which is a further indication of the fact that 2P particles are found mostly bonded with other 2P particles, thus forming chains. Finally, also from g_HS–2P(r), we find that the number of 2P neighbors of one HS particle decreases with x_2P as expected. Due to the strong attraction between 2P particles, we are then interested to understand whether the two species could undergo demixing for some concentrations of 2P and HS particles.

An indication of demixing is provided by the concentration–concentration structure factor S_cc(q). In fact, close to the spinodal line, there is usually a strong increase in concentration fluctuations within the mixture, which is accompanied by a divergence of lim_q→0 S_cc(q).⁵⁰ The concentration–concentration structure factor for the 2P–HS mixture at different concentrations is reported in Fig. 3(b) (lower-right panel). We observe that on decreasing x_2P and, hence, on increasing x_HS, weak oscillations appear in S_cc(q), which progressively increase. At the same time, the low q part of the curve increases above 1 for x_2P = 0.25, signaling the appearance of inhomogeneities. A closer look at the snapshots in Fig. 3(c) helps to clarify what is happening in the binary mixtures. As discussed above, the radial distribution functions indicate that the 2P particles are always surrounded by at least two other 2P particles, thus satisfying their bonds. At high x_2P, this can be achieved while still allowing particles to be homogeneously distributed throughout the simulation box due to the high ϕ at which the system is simulated. On decreasing x_2P, this is not possible anymore, and to maximize the bonds, 2P particles form few chains embedded in the sea of HS particles, as shown in the snapshots of Fig. 3(c). The inhomogeneous distribution in the simulation box of 2P particles is indicated by the increase at small q values of S_cc(q). Although this could be interpreted as demixing, we instead relate this result to the intrinsic behavior that 2P particles have at low T and low ϕ. Indeed, it has been shown that, under these conditions, the system of 2P particles can be treated as an ideal gas of chains, rather than particles, and an approximate expression of S(q) at low q values has been formulated starting from the S(q) of a single chain.⁵⁹ Since the 2P chains are uniformly distributed throughout the box, we consider x_2P = 0.25, T = 0.1, and ϕ = 0.262 as a stable state point of the mixture.

We focus our investigation on the state point with T = 0.1 and ϕ = 0.262 since, as discussed above, for low enough T and high enough ϕ, the 2P particles have a tendency to form aligned chains. Under these conditions, the competition between the bonding and excluded volume entropy leads to non-trivial effects in terms of effective interactions. Our aim is to understand whether or not the contribution to the effective potential stemming from the two species can be simply obtained from the linear combination of the effective potentials of the single-component co-solutes. In other words, given the effective potential between two large HS colloids immersed in a mixture of the 2P and HS co-solute with concentration x_2P and x_HS = 1 − x_2P, i.e., $Vx2P,xHSeff(r)$⁠, we want to verify whether

holds. We therefore compare the effective potential of the two colloids immersed in the mixture, for different concentrations of the two species, with the corresponding ones calculated from Eq. (10). The results are shown in Fig. 4. The gray curves refer to V_eff calculated for the single-component HS and 2P co-solutes and are used as reference potentials. Note that on increasing the concentration of 2P particles, the effective potential becomes totally attractive as expected. In addition, for all the studied concentrations, almost perfect agreement is found with the potential obtained from Eq. (10), suggesting the additivity of the total effective interactions of the mixture.

On the basis of this result for the total interactions, we would also expect that each contribution to the total potential, arising from the two species in the mixture separately, coincides with the corresponding single-component one. To shed light on this point, we cannot rely on the umbrella sampling method, which does not provide the two terms separately, but we resort to the virtual displacement method explained in Sec. II. This method allows us to have an independent measurement of the effective force between the two colloids, which we compare with that obtained from the umbrella sampling method via the relation

In addition, and most importantly, the virtual displacement method provides the effective forces stemming from the two components of the mixture separately. Figure 5 (left panels) shows the total effective force mediated by the mixtures compared with the force calculated from the umbrella sampling method and that obtained from the linear combination of Eq. (10) (where, instead of V_eff, we use F_eff). Again, we find that the agreement is very good, independently of the composition of the mixture.

In addition, the central and right panels of Fig. 5 report the components of the effective force due to the HS species (⁠$βFxHSeff$⁠) and the 2P species (⁠$βFx2Peff$⁠), respectively. Both forces are compared with the effective forces stemming from the single-component co-solutes, multiplied by the concentration of the species, i.e., $xHSβFHSeff$ and $x2PβF2Peff$⁠, respectively, which are indeed the two terms in Eq. (10). Interestingly, we find that the separate contributions of the two species do not match their corresponding term in Eq. (10), even if their sum satisfies the same equation. In particular, while the HS force $βFxHSeff$ is larger than $[xHSβFHSeff]$⁠, the 2P force $βFx2Peff$ is smaller than $[x2PβF2Peff]$⁠. This result is found for all mixture compositions, but the deviations are largest at x_2P = 0.5 and significantly decrease with increasing 2P concentration. For x_2P = 0.875, they become negligible.

Since the virtual displacement method is based on counting the number of overlaps of the two colloids with the co-solute in between and outside them, the behavior of the two components of the effective force indicates that there are more HS co-solute particles between the two colloids than expected from $xHSβFHSeff$⁠. Consequently, there are less 2P co-solute particles in this region than expected from $x2PβF2Peff$⁠. This is true for almost all concentrations when the colloids are at a distance r/σ ≲ 1. For x_2P ≥ 0.5, this also occurs at larger distances. While the HS co-solute particles satisfy only entropic constraints, the 2P co-solute behavior is always determined by the competition of bonding energy and excluded volume. In particular, we expect that the lower number of 2P particles is related to a reduced ability to form bonds in between the two colloids (at least when r/σ ≲ 1). Thus, an increasing concentration of HS particles in this region makes the 2P particles stay outside the colloids where they still are able to form bonds. To test this hypothesis, we calculate the average number of bonds of 2P particles that are in contact with the colloids, i.e., considering only the 2P particles exerting a force on the HS colloids. The results for this observable are shown in Fig. 6(a) as a function of the colloids’ surface-to-surface distance r, distinguishing the 2P particles in two groups: those acting on the colloids from the outside and those exerting a force from the in between region, as determined during the virtual displacement method (cf. Sec. II). The calculation has been performed for the co-solute mixture with 50%HS–50%2P [upper panel of Fig. 6(a)], where this effect should be strongest. We compare the results to those obtained for the single-component 2P co-solute at the same packing fraction and temperature [lower panel of Fig. 6(a)]. The gray area drawn between the two curves in the upper panel of Fig. 6(a) highlights that there is a gap between the number of bonds per particle in the two regions. Although the difference is small, the gap is clearly absent in the case of the single-component 2P co-solute, indicating that HS particles influence the ability of 2P particles to form bonds in the region between the two colloids. Finally, for both cases, we also highlight (purple area) a small region for r/σ > 1 where the number of bonds formed between the two colloids is larger than that formed outside them. This is related to the nematization effect discussed in Sec. III A: when the depletion zone disappears and 2P particles can access the region between the colloids, they can satisfy both entropic and energetic constraints by reorienting their bonds in order to form parallel chains.

We repeat the analysis discussed above for the case of two HS colloids immersed in a mixture of co-solute particles made of HS and three-patch (3P) particles. Previous work³² has shown that the effective interactions generated by 3P co-solute particles at high ϕ and low T exhibit the same features as the 2P ones, i.e., 3P particles give rise to an oscillatory effective attraction between HS colloids. Again, the attractive oscillations are related to the orientation of particle bonds, which, in this case, tend to form planar sheets of 3P particles under confinement, which are orthogonal to the vector that joins the centers of the two colloids. Differently from 2P particles, the phase diagram of 3P particles contains a gas–liquid coexistence region with a critical point that, for the model parameters employed in this work, has been located at ϕ_c = 0.072 and T_c = 0.125.⁶⁰ Here, we investigate the mixture of HS and 3P particles at T = 0.125 and ϕ = 0.262 for which the effective interaction mediated by the two species separately is already known.³² Figure 6(b) shows the total effective force and the contribution to the force from 3P and HS co-solutes at two different concentrations of the two species (x_3P = 0.875 and x_3P = 0.750). Interestingly, in this case, upon increasing the concentration of HS, we find a progressive deviation from the additive behavior of the total force. In addition, the two components of the effective force display a long-range tail, which clearly signals the presence of some kind of correlation between the co-solutes. It is worth to stress that, at the two studied compositions, the mixtures are stable with respect to demixing in the absence of the two colloids but are not too far from the onset of this transition. This is shown in Fig. 6(e), where demixing clearly takes place for x_3P = 0.500, as signaled by the low q divergence of the concentration–concentration structure factor. Thus, the results obtained for the effective forces can be interpreted in terms of a local phase separation driven by the presence of the two colloids in the mixture. This hypothesis, also explaining the long-range components of the effective force, is supported by the snapshot of the mixture in the presence of the HS colloids shown in Fig. 6(c), where the region between the two colloids is devoid of 3P particles even for r/σ > 1. Consequently, only the HS co-solute particles are found in this region, indicating local demixing.

Finally, to show that such a non-additive behavior disappears far from demixing, we have studied a 3P–HS co-solute mixture at a higher temperature (T = 0.2) at the same packing fraction. The results are shown in Fig. 6(d) for two different mixture compositions. As for the 2P and HS mixtures, in this case, the total effective force coincides with that arising from the linear combination of the two single-component contributions. In addition, the effective force generated by HS co-solute particles is always higher than the expected one and it is balanced by that of the 3P particles, which is smaller than expected, in agreement with the findings for the 2P and HS mixture results, shown in Fig. 5. Also in the present case the two effects cancel out such that the total force is found to be additive.

In this work, we have investigated the effective interactions of two HS colloids immersed in a binary mixture of smaller co-solute particles. The binary mixtures considered contain HS and patchy particles (2P or 3P). 2P particles are known to reversibly assemble into chains, and a previous study on the effective interactions between colloids mediated by those particles reported that the competition between bonding and entropic constraints modulates the effective potential, giving rise to attractive oscillations between the colloids.³² This occurs mainly at low T and high ϕ conditions, where the present analysis is mostly focused. When a mixture of 2P and HS co-solutes is considered, the effective potential mediated by the co-solute shows a continuous transition from a standard effective potential found in highly asymmetric HS binary mixtures to a completely attractive potential as observed for HS colloids in the pure 2P co-solute. Despite the non-trivial effects driven by the assembly of 2P particles and their confinement close to the colloids, we find that, when two HS colloids in a mixture of HS and 2P co-solutes are considered, the effective interactions of the mixture are always additive and, thus, they can simply be estimated from the linear combination of the effective interactions mediated by the pure HS and the pure 2P co-solute. Interestingly, we find that the separate contributions to the effective force from both HS and 2P particles deviate from the expected ones, but these deviations cancel out, providing a total effective force that satisfies the linear combination. This is attributed to the fact that a larger number of 2P particles are found outside the two colloids than in the confined region between them and vice versa for the HS particles. We have shown that this is related to the fact that 2P particles are able to form more bonds outside the colloids, while HS particles are better accommodated in the confined region between the two colloids when their surface-to-surface distance is r < σ.

We have repeated the same analysis for HS–3P co-solute mixtures, finding that, for some concentrations, the linear combination does not work to a satisfactory degree. This behavior can be attributed to local demixing induced by the presence of the colloids, despite the fact that no signs of demixing have been found in the mixture without the HS colloids, even if the considered state point is close to the demixing line. When the two colloids are placed in the mixture, we find that 3P particles are fully depleted from the confined region between the colloids even for r > σ. However, if the temperature is increased, the linear combination again yields the correct result also for the HS–3P mixtures.

In conclusion, we have presented a first survey of effective interactions mediated by binary co-solute mixtures in which one component is capable of self-assembly. The present findings open the way to use additive effective interactions in the presence of binary co-solutes. Indeed, if the linear combination is valid, knowledge of the potential in the single-component co-solutes will be sufficient to predict the total effective potentials in the mixtures. Hence, this approach provides a promising route for the inverse design of new kinds of colloid–colloid interactions in silico, which is able to induce peculiar collective behavior of the colloids, using a method that is, in principle, also readily realizable in experiments. In this context, it will also be fundamental to characterize many-body effects that certainly play a role in the formation of colloidal superstructures as discussed in some recent studies.^61,62

Additionally, it has to be understood whether additivity would still be valid for binary co-solutes of different sizes or different types (e.g., inverse patchy colloids⁶³) and how different co-solute–colloid interactions would affect it. Simulations, as the ones presented here, as well as state-of-the-art DFT^56,64–67 and recent integral equation theory calculations^68,69 will surely provide insights in this direction.

P.H.H. acknowledges the Austrian Science Fund FWF for support (Erwin Schrödinger Fellowship, Grant No. J3811 N34).

The data that support the findings of this study are available from the corresponding author upon reasonable request.