A Virial expansion for multicomponent systems

Depletion effects in colloid-polymer solutions


The surface tension, the adsorption, and the depletion thickness of polymers close to a single nonadsorbing colloidal sphere are computed by means of Monte Carlo simulations. We consider polymers under good-solvent conditions and in the thermal crossover region between good-solvent and behavior. In the dilute regime we consider a wide range of values of , from (planar surface) up to -50, while in the semidilute regime, for ( is the polymer concentration and is its value at overlap), we only consider and 2. The results are compared with the available theoretical predictions, verifying the existing scaling arguments. Field-theoretical results, both in the dilute and in the semidilute regime, are in good agreement with the numerical estimates for polymers under good-solvent conditions.

I Introduction

The study of the fluid phases in mixtures of colloids and nonadsorbing neutral polymers has become increasingly important in recent years; see Refs. (1); (2); (3); (4); (5); (6) for recent reviews. These systems show a very interesting phenomenology, which only depends to a large extent on the nature of the polymer-solvent system and on the ratio , where is the zero-density radius of gyration of the polymer and is the radius of the colloid. Experiments and numerical simulations indicate that polymer-colloid mixtures have a fluid-solid coexistence line and, for , where (1) -0.4, also a fluid-fluid coexistence line between a colloid-rich, polymer-poor phase (colloid liquid) and a colloid-poor, polymer-rich phase (colloid gas). On the theoretical side, research has mostly concentrated on mixtures of neutral spherical colloids and polymers in solutions under good-solvent or conditions. In the former case, predictions for the colloid-polymer interactions have been obtained by using full-monomer representations of polymers (for instance, the self-avoiding walk model was used in Refs. (7); (8)), field-theoretical methods (9); (10); (11), or fluid integral equations (12). Moreover, some general properties have been derived by using general scaling arguments (13); (14); (15). At the point, the analysis is simpler, since polymers behave approximately as ideal chains. These theoretical results have then been used as starting points to develop a variety of coarse-grained models and approximate methods, see Refs. (16); (7); (5) and references therein, which have been employed to predict colloid-polymer phase diagrams.

In this paper, we consider polymer solutions that either show good-solvent behavior or are in the thermal crossover region between good-solvent and conditions. We study the solvation of a single colloid in the solution, assuming that the monomer-colloid potential is purely repulsive. We determine the distribution of the polymer chains around a single colloidal particle, which is the simplest property that characterizes polymer-colloid interactions. We investigate numerically, by means of Monte Carlo simulations, how it depends on the quality of the solution, which is parametrized (17); (18) in terms of the second-virial combination , where is the second virial coefficient and is the zero-density radius of gyration. This adimensional combination varies between 5.50 in the good-solvent case (19) and zero ( point). Beside good-solvent solutions, we consider two intermediate cases: solutions such that is approximately one half of the good-solvent value, which show intermediate properties betweeen good-solvent and behavior, and solutions such that is 20% of the good-solvent value, which are close to the point. In each case we compute the polymer density profile around the colloid. These results are used to determine thermodynamic properties, like the surface tension and the adsorption (20), which are then compared with the available theoretical predictions. Note that an analysis of the polymer depletion around a colloid in the thermal crossover region was already performed in Refs. (21); (22), but without a proper identification of the universal crossover limit (23). Here, we wish to perform a much more careful analysis of the crossover behavior, following Refs. (24); (18). We focus on the dilute and semidilute regimes, in which the monomer density is small and a universal behavior, i.e., independent of chemical details, is obtained in the limit of large degree of polymerization. In the dilute regime, in which polymer-polymer overlaps are rare, solvation properties are determined for a wide range of values of , from 0 up to 30-50. In the semidilute regime, simulations of systems with large require considering a large number of colloids, which makes Monte Carlo simulations very expensive. Hence, we only present results for and 2.

The paper is organized as follows. In Sec. II we define the basic quantities we wish to determine. First, in Sec. II.1 and II.2 we introduce the surface tension, the adsorption, and the depletion thickness, and discuss their relation to the density profile of the polymers around the colloid. In Sec. II.3 we discuss the low-density behavior and the relation between solvation properties and colloid-polymer virial coefficients that parametrize the (osmotic) pressure of a polymer-colloid binary solution in the low-density limit. Sec. II.4 discusses the different behavior that are expected as a function of and density and gives an overview of theoretical and numerical predictions. Sec. III summarizes our polymer model and gives a brief discussion on how one can parametrize in a universal fashion the crossover between the good-solvent and the behavior (for more details, see Ref. (18)). In Sec. IV we present our results for the dilute regime, while finite-density results are presented in Sec. V. Sec. VI discusses a simple coarse-grained model in which each polymer is represented by a monoatomic molecule, which represents a more rigorous version of the well-known Asakura-Oosawa-Vrij model. In Sec. VII we present our conclusions. Two appendices are included, one explaining how to compute the virial coefficients of a binary mixture of flexible molecules, and one discussing the small- behavior of the virial coefficients. Tables of results are reported in the supplementary material.

Ii Adsorption and depletion thickness

ii.1 Definitions

Let us consider a solution of nonadsorbing polymers in the grand canonical ensemble at fixed volume and chemical potential . Temperature is also present, but, since it does not play any role in our discussion, we will omit writing it explicitly in the following. Let us indicate with the corresponding grand potential. Let us now add a spherical colloidal particle of radius to the solution and let be the corresponding grand potential. The insertion free energy can be written as the sum of two terms, one proportional to the volume of the colloid and one proportional to its surface area (20):


where is the bulk pressure and is the surface tension. The latter quantity can be related to the adsorption defined in terms of the change in the mean number of polymers due to the presence of the colloid:


where indicates the numbers of polymers present in the solution, and are averages in the presence and in the absence of the colloidal particle, respectively. Differentiating Eq. (1) with respect to , we obtain


We also define the average bulk polymer density as


The surface tension can also be defined in the canonical ensemble, as a function of the bulk polymer density . Then, we have




for the bulk system ( as usual).

The adsorption coefficient can be easily related to the polymer density profile around the colloid. Assume that each polymer consists of monomers and define the average bulk monomer density . Then, we write


where the average is performed at chemical potential and volume , is the colloid position, and , , , are the monomer positions. If we now define the monomer-colloid pair correlation function


and the integral


we obtain


Since for , a more transparent relation is obtained by defining


in which one only integrates the density profile outside the colloid. Since , we have


In the previous discussion we have considered the monomer-colloid correlation function, but it is obvious that any other polymer-colloid distribution function could be used. In order to compare our results with those obtained in coarse-grained models (we will discuss them in Sec. VI), we will also use the pair distribution function between the colloid and the polymer centers of mass. If , , are the positions of the monomers belonging to polymer , we first define the polymer center of mass


Then, the pair distribution function between a colloid and a polymer center of mass is defined by


where the average is taken at a given value . In terms of this quantity


If we define1 , we can write a relation analogous to Eq. (12). Comparison of Eqs. (12) and (15) implies and , hence in the following we will simply refer to these quantities as and .

It is interesting to relate the pair correlation functions to the analogous correlation functions that are appropriate for a binary system consisting of polymers and colloids at polymer and colloid chemical potentials and , respectively. Indeed, one can show that, in the limit , i.e., when the colloid density goes to zero, one has


Eq. (16) allows us to relate to thermodynamic properties of the binary mixture in the limit of vanishing colloid density. For this purpose we use the Kirkwood-Buff relations between structural and thermodynamic properties of fluid mixtures (26); (25). The integral , which is relevant to determine adsorption properties, corresponds to one of the Kirkwood-Buff integrals (26); (25) defined as


where and label the different species of the mixture. The integrals can be related to derivatives of the pressure with respect to the polymer and colloid densities. For we have (26); (25)


which imply


Eqs. (12) and (5) can then be rewritten as


ii.2 Depletion thickness

Depletion effects can be equivalently parametrized by introducing the depletion thickness (29); (28); (5); (27), which is an average width of the depleted layer around the colloid. It is defined in terms of the integral as


so that


Since is only determined by , knowledge of is completely equivalent to that of the adsorption. The two quantites are related by


As we shall discuss below, for large polymer densities, hence in this limit


It is interesting to discuss the limit , in which the colloid degenerates into an impenetrable plane. Setting in Eq. (11), we obtain


For , we have , where is the pair distribution function between an impenetrable plane at and a polymer. Then, we obtain for




Taking the limit in Eqs. (24) and (25), we obtain


ii.3 Low-density expansions

For the depletion thickness and the surface quantities and can be related to the virial coefficients that parametrize the expansion of the pressure of a binary colloid-polymer system in powers of the concentrations. These relations have already been discussed in the literature (30); (31); (32); (33). They can be easily derived by using Eqs. (21) and (22). We start by expanding the pressure as


where and are the colloid and polymer concentrations and we have neglected fourth-order terms. Then, Eqs. (21) and (22) give


In the limit one should recover the results for an infinite impenetrable plane. This requires the coefficients appearing in the previous two expressions to be of order as . This is explicitly checked in App. B and allows us to write


Explicit expressions for and are reported in Appendix B.

For the depletion thickness we obtain


where we have defined the polymer volume fraction


and the adimensional combinations and , where is the zero-density polymer radius of gyration. In the limit , we should obtain the density expansion of the depletion thickness for an impenetrable plane. Using Eq. (30) we obtain


ii.4 Theoretical predictions and scaling arguments

Depletion properties have been extensively studied in the past. Here we present scaling arguments and literature results, that will be checked in the following sections by using our accurate Monte Carlo estimates.

For an ideal (noninteracting) polymer solution the insertion free energy is exactly known (9):


where is the zero-density radius of gyration. The depletion thickness follows immediately (5); (27):


For good-solvent polymers there are several predictions obtained by using the field-theoretical renormalization group. In the dilute limit , the surface tension has been determined (10) both in the colloid limit in which and in the so-called protein limit . Setting and in the results of Ref. (10), we obtain for and


Note that the dilute behavior in the colloidal regime is similar to that observed in the ideal case. The coefficients corresponding to the planar term and to the leading curvature correction are close, while the second-curvature correction is absent in the ideal case and quite small for good-solvent chains.

In the opposite limit general arguments predict (13); (10)


The constant has been estimated by Hanke et al. (10):


Eq. (5) gives then . For the depletion thickness we obtain .

Finite-density corrections have been computed by Maassen et al. (11) in the renormalized tree approximation. For they obtain


In this approximation one does not recover the correct large- behavior (42), hence we expect it to be valid only in the colloid regime. The zero-density behavior can be compared with that given in Eq. (41), which includes the leading (one-loop) correction. Differences are small, of order 5%. We expect an error of the same order for the coefficients of the density correction.

The behavior of in the semidilute regime is expected to have universal features. If the polymer volume fraction is large, we expect, on general grounds, the behavior (14)


where is an exponent to be determined and is a coefficient, which a priori can depend on . However, deep in the semidilute regime, the coil radius of gyration is no longer the relevant length scale. One should rather consider the density-dependent correlation length (34), which measures the polymer mesh size. The scaling behavior (45) should be valid for , and in this regime plays no role. Therefore, is not the relevant parameter and is independent of . To determine the exponent we can use the same argument which allows one to determine the scaling behavior of the osmotic pressure in the semidilute regime. For large we expect thermodynamic properties to depend only on the monomer density and not on the number of monomers per chain. This requirement gives (14)


where (35) . Predictions (45) and (46) can also be obtained (14) by noting that can only depend on the correlation length deep in the semidilute regime, i.e., when . Then, dimensional analysis gives


Using (34); (36); (37), we obtain again Eq. (45) with given by Eq. (46). Eq. (45) allows us to obtain the large- behavior of the adsorption and of the depletion thickness. Using Eq. (5) and the general scaling of the osmotic pressure (34); (36); (37) , we obtain


Equivalently, one could have observed that , since is the only relevant length scale. Using , we obtain . Eq. (26) implies then Eq. (48).

The large- behavior was determined in the renormalized tree-level approximation obtaining (11)


This result is fully consistent with Eq. (45), since the correction appearing in Eq. (49) vanishes for . The exponent of the -dependent correction in Eq. (49) can be easily interpreted. Consider the ratio . This quantity is adimensional, hence it is a universal function of adimensional ratios of the relevant length scales in the system. Deep in the semidilute regime the relevant polymer scale is the correlation length , hence we expect


Now we take large so that . Then, we can expand


Since , we obtain


which reproduces the behavior (49). Eq. (51) is the semidilute analogue of the Helfrich expansion in powers of that holds for . The only difference is the expansion variable: in the semidilute region, polymer size is characterized by , hence one should consider instead of .

Quantitative predictions for the large- behavior of and can be derived from Eq. (49), by using Eq. (5) and the large- behavior of . The latter can be derived from the results of Ref. (38), which give for . Thus, we obtain


In the protein limit, in which is large, beside the regime in which Eqs. (45), (51) and (53) hold, there is a second interesting regime in which one has both and . For large, these conditions are satisfied both in the dilute limit and in the semidilute region, as long as is not too large. Under these conditions, Eq. (42) holds irrespective of the polymer density. Therefore, Eq. (22) can be rewritten as


For , the pressure term can be neglected compared with the term proportional to , hence the right-hand side is density independent. This implies that the integrand that appears in left-hand side is also density independent in the density region where and is equal to . For , using the virial expansion (31) we can write


Therefore, we can identify . Moreover, vanishes for (a similar result holds for the higher-order virial coefficients). By using Eq. (5) and Eq. (42) we also predict for the adsorption


Since for large (34), this relation predicts . Note that Eq. (57) holds only for . As further increases, decreases and one finds eventually . Then, Eq. (57) no longer holds and a crossover occurs. For the asymptotic behavior sets in. Eq. (57) can be written in a more suggestive form, by noting that (34). Hence


We recover the same scaling that occurs in the dilute regime, with replacing as relevant polymer scale.

Figure 1: Different regimes for the depletion thickness in terms of and . The functions and have a regular expansion in powers of their argument.

To conclude, let us summarize the different types of behavior of the depletion thickness in the - diagram for the good-solvent case. They depend on the relative size of the three different scales that appear in the problem: the radius of gyration of the polymer, the radius of the colloid and the correlation length . In the colloid regime in which , i.e. , depletion shows two different behaviors, depending on the ratio . In the dilute regime in which the relevant scale is the radius of gyration (domain I in Fig. 1), is of order with a proportionality constant that can be expanded in powers of (Helfrich expansion). If instead (semidilute regime, domain III in Fig. 1), the relevant scale is the correlation length . The depletion thickness is proportional to with a proportionality constant that admits an expansion in powers of . Since for , the limiting behavior is independent of the colloid radius. In the protein regime in which , i.e., , depletion shows three different behaviors. In the dilute regime (domain II in Fig. 1), , i.e., is much larger than the colloid radius but much smaller than . In the semidilute regime, two different behaviors occur. If (domain IV), the role of the radius of gyration is now assumed by the correlation length and we have . Finally, as increases further, one finally finds and one observes again (domain III).

The surface tension was also computed in the PRISM approach (12), obtaining


Such an expression does not have the correct behavior for or . In the dilute regime and for small , comparison with the field-theoretical results (we shall show that they are quite accurate) shows that it only provides a very rough approximation, differences being of order 20-30%.

The adsorption was computed numerically for the planar case () in Ref. (7), obtaining


This expression allows us to compute for using the expression of the compressibility factor given in Ref. (38). In the small-density limit we obtain


while for we obtain


We can compare these expressions with the field-theory results. The leading density correction in Eq. (61) is approximately one half of that predicted by field theory, see Eq. (44), while the large- expression (62) predicts a surface tension that is 17% smaller than Eq. (49).

Finally, we mention the phenomenological expression for the depletion thickness of Fleer et al. (28); (5)


which should be only valid in an intermediate range of values of (28); (27), since it does not have the correct behavior in the limits and .

There are no predictions for polymers in the thermal crossover region. In this case, a new scale comes in, the dimension of the so-called thermal blob (34). On scales , the polymer behaves as an ideal chain, hence for the surface tension should coincide with that appropriate for an ideal chain. This implies that for any finite value of we should recover the ideal result for the surface tension, provided that is large enough. In particular, we predict


for all finite values of and . In practice, Eq. (42) holds also for finite , with the values appropriate for the ideal chain, and . If instead , we expect to observe a nontrivial crossover behavior. Its determination is one of the purposes of the present paper.

Iii Polymer model and crossover behavior

In order to determine full-monomer properties, we consider the three-dimensional lattice Domb-Joyce model (39). We consider chains of monomers each on a finite cubic lattice of linear size with periodic boundary conditions. Each polymer chain is modeled by a random walk with (we take the lattice spacing as unit of length) and . The Hamiltonian is given by


where is the Kronecker delta. Each configuration is weighted by , where is a free parameter that plays the role of inverse temperature. This model is similar to the standard lattice self-avoiding walk (SAW) model, which is obtained in the limit . For any positive , this model has the same scaling limit as the SAW model (39) and thus allows us to compute the universal scaling functions that are relevant for polymer solutions under good-solvent conditions. In the absence of colloids, there is a significant advantage in using Domb-Joyce chains instead of SAWs. For SAWs scaling corrections that decay as (, Ref. (35)) are particularly strong, hence the universal, large–degree-of-polymerization limit is only observed for quite large values of . Finite-density properties are those that are mostly affected by scaling corrections, and indeed it is very difficult to determine universal thermodynamic properties of polymer solutions for by using lattice SAWs (38). These difficulties are overcome by using the Domb-Joyce model for a particular value of (40); (24), . For this value of the repulsion parameter, the leading scaling corrections have a negligible amplitude (40); (24), so that scaling corrections decay faster, approximately as . As a consequence, scaling results are obtained by using significantly shorter chains. Unfortunately, in the presence of a repulsive surface, new renormalization-group operators arise, which are associated with the surface (41). The leading one gives rise to corrections that scale as (41), where is the Flory exponent (an explicit test of this prediction is presented in the supplementary material), hence it spoils somewhat the nice scaling behavior observed in the absence of colloids. Nonetheless, the Domb-Joyce model is still very convenient from a computational point of view. Since interactions are soft, the Monte Carlo dynamics for Domb-Joyce chains is much faster than for SAWs. We shall use the algorithm described in Ref. (38), which allows one to obtain precise results for quite long chains () deep in the semidilute regime.

The Domb-Joyce model can also be used to derive the crossover functions that parametrize the crossover between the good-solvent and -point regimes, at least not too close to the point, see Refs. (24); (18) for a discussion. Indeed, if one neglects tricritical effects, which are only relevant close to the point (42), this crossover can be parametrized by using the two-parameter model (36); (37); (23). Two-parameter-model results are obtained (43) by taking the limit , at fixed . The variable interpolates between the ideal-chain limit () and the good-solvent limit (). Indeed, for the Domb-Joyce model is simply the random-walk model, while for any and one always obtains the good-solvent scaling behavior. The variable is proportional to the variable that is used in the context of the two-parameter model. We normalize as in Refs. (24); (18), setting


Note that the crossover can be equivalently parametrized (17); (40); (44); (24); (18) by using the second-virial combination ( is the zero-density radius of gyration), which varies between the good-solvent value (19) and at the point. With normalization (66) we have for small (43); (40). The correspondence between and in the whole crossover region is given in Ref. (24).

As discussed in Ref. (24), the two-parameter-model results can be obtained from Monte Carlo simulations of the Domb-Joyce model by properly extrapolating the numerical results to . For each we consider several chain lengths . For each of them we determine the interaction parameter by using Eq. (66), that is we set . Simulations of chains of monomers are then performed setting . Simulation results are then extrapolated to , taking into account that corrections are of order (43); (40).

In this paper we have performed a detailed study of the depletion for two values of : and , which correspond to (24) and . They correspond to polymer solutions of intermediate quality. Since (19) under good-solvent conditions, we have and 0.54 for and , respectively. Hence, for we are quite close to the point, while is intermediate between the good-solvent and regimes.

In this paper we discuss depletion effects close to neutral colloids, which are modelled as hard spheres that can move everywhere in space: their centers are not constrained to belong to a lattice point. This choice is particularly convenient since it drastically reduces lattice oscillations in colloid-polymer correlation functions. Such oscillations are instead present if colloids are required to sit on lattice points, as was done in Ref. (22). Colloids and monomers interact by means of a simple exclusion potential. If and are the coordinates of a monomer and of a colloid, we take as interaction potential


Iv Dilute behavior

As we have seen in Sec. II.3, the low-density behavior of the surface tension or, equivalently, of the depletion thickness can be obtained by computing the virial coefficients and . We will thus report the computation of these two quantities and also of , which would be relevant to characterize the effective interaction between two colloids in a dilute solution of polymers. Then, we shall discuss the depletion thickness for and its first density correction.

iv.1 Virial coefficients

(FM) (FM) (SB) (SB)
1.0605(3) 1.0514(2) 4.44(1)
1.1066(1) 1.1065(3) 2.387(8)
1.1221(4) 1.1222(3) 0.771(4)
Table 1: Estimates of the universal surface combinations and . We report full-monomer (FM) and single-blob (SB) results (see Sec. VI).

To determine the virial coefficients under good-solvent conditions we have simulated the Domb-Joyce model at . We consider chains of length for , and , 24000 to derive the results corresponding to . Long chains are needed for large values of to ensure that the colloid radius is somewhat larger than the lattice spacing. Virial coefficients are determined as explained in App. A. The universal extrapolations of the finite- results for the adimensional combinations and are explicitly reported in the supplementary material. In the case of the two-parameter model, we have considered for (for both and ) and , 30000 for (only for ). The results at the same value of have then been extrapolated taking into account the scaling corrections. Results are reported in the supplementary material. We also computed the adimensional combinations and , which parametrize the depletion thickness in the presence of an impenetrable planar surface, see Table 1. The behavior of the adimensional combinations for is discussed in detail in Appendix B. We have

Figure 2: Combination vs . For , such a combination converges to 1 for all values of . We also report the approximate large- behavior.

For large values of , we have , a behavior which can be derived by means of a blob argument (22) or from the large- behavior of , as discussed in Sec. II.4. More precisely, we predict , where is the constant parametrizing the large- behavior of , defined by Eq. (42). Note that this relation holds both in the good-solvent regime with and in the crossover regime with . As for the third virial coefficient , we have already shown that vanishes as . Hence, we expect


with for . We have been unable to predict the value of . A numerical fit of the data indicates , both in the good-solvent limit and in the crossover region, see Fig. 2. As for , a blob argument implies