As first explained by the classic Asakura–Oosawa (AO) model, effective attractive forces between colloidal particles induced by depletion of nonadsorbing polymers can drive demixing of colloid–polymer mixtures into colloid-rich and colloid-poor phases, with practical relevance for purification of water, stability of foods and pharmaceuticals, and macromolecular crowding in biological cells. By idealizing polymer coils as effective penetrable spheres, the AO model qualitatively captures the influence of polymer depletion on thermodynamic phase behavior of colloidal suspensions. In previous work, we extended the AO model to incorporate aspherical polymer conformations and showed that fluctuating shapes of random-walk coils can significantly modify depletion potentials [W. K. Lim and A. R. Denton, Soft Matter **12**, 2247 (2016); J. Chem. Phys. **144**, 024904 (2016)]. We further demonstrated that the shapes of polymers in crowded environments sensitively depend on solvent quality [W. J. Davis and A. R. Denton, J. Chem. Phys. **149**, 124901 (2018)]. Here, we apply Monte Carlo simulation to analyze the influence of solvent quality on depletion potentials in mixtures of hard-sphere colloids and nonadsorbing polymer coils, modeled as ellipsoids whose principal radii fluctuate according to random-walk statistics. We consider both self-avoiding and non-self-avoiding random walks, corresponding to polymers in good and theta solvents, respectively. Our simulation results demonstrate that depletion of polymers of equal molecular weight induces much stronger attraction between colloids in good solvents than in theta solvents and confirm that depletion interactions are significantly influenced by aspherical polymer conformations.

## I. INTRODUCTION

In the field of condensed matter physics, the Asakura–Oosawa (AO) model of colloid–polymer mixtures^{1,2} has a status akin perhaps to the van der Waals model of fluids, the Einstein model of solids, the Ising model of magnets, and the primitive model of electrolytes. The AO model, inspired by deep physical insight, first identified polymer depletion as the basic mechanism underlying effective attraction between colloidal particles induced by nonadsorbing polymers. In soft, fragile materials, depletion-induced interactions^{3–7} between mesoscopic particles typically compare in magnitude to thermal energies and thus can strongly influence self-assembly and thermodynamic phase stability. In this way, the AO model qualitatively explains the observed phase behavior of colloid–polymer mixtures, in particular, demixing into colloid-rich and colloid-poor bulk phases.

In many practical applications, such as in stabilizing foods^{8,9} and pharmaceuticals against coagulation or preventing the aggregation of proteins,^{10,11} it is important to minimize depletion-induced attraction. In other applications, such as in purifying water by promoting flocculation of colloidal impurities^{12} or in guiding the self-assembly of virus particles,^{13,14} amplifying the effects of polymer depletion is instead desirable. Depletion also contributes to macromolecular crowding and segregation of biopolymers within biological cells.^{15–22}

In its original form, the AO model depicts polymer coils as effective spheres, of fixed size defined by the radius of gyration, that are mutually penetrable, but impenetrable to colloidal particles due to excluded-volume interactions. The model reveals that depletion of polymers from the space between hard colloidal surfaces creates an imbalance in polymer concentration, and thus in osmotic pressure, that drives effective attraction between colloids. Equivalently, configurations in which excluded-volume shells of neighboring colloids overlap maximize the free volume available to polymer coils and thus are entropically favored.

Although it captures the essence of polymer depletion, the AO model omits certain important aspects of real physical systems. Most obviously, by idealizing polymer coils as effective spheres of unvarying size, the model neglects the internal degrees of freedom—structure and flexibility—of polymers in solution. In biological systems, the structure associated with folding (or misfolding) of proteins determines the function of such biopolymers with relevance for many diseases.

The realization that polymers are flexible, aspherical objects predates the AO model by at least two decades. Kuhn argued^{23} that linear polymer coils in solution can be modeled as random walks (RWs) with fluctuating shapes that approximate those of elongated, flattened ellipsoids (in their principal-axis frame). The insight that the end-to-end path of a polymer is a physical manifestation of a random walk has spurred many studies of shapes of random walks.^{24–38} As a vital example, the shapes of RNA, DNA, and proteins are important for cellular processes in the crowded environment of biological cells,^{39–44} translocation of polymers through narrow pores,^{45,46} and packaging of DNA in viral capsids.^{47}

Depletion forces and their impact on polymer crowding and phase behavior in colloid–polymer mixtures have been probed by neutron scattering,^{48–54} atomic force microscopy,^{55} total internal reflection microscopy,^{56} optical trapping,^{57–59} and turbidity measurements,^{60–62} to name but a few experimental methods. Modeling studies of colloid–polymer mixtures have used scaling and mean-field free-volume theories,^{63–71} force-balance theory,^{72} perturbation theory,^{73,74} polymer field (renormalization group) theories,^{75–82} integral-equation theories,^{83–87} density-functional theories,^{88–92} adsorption theories,^{93–95} and computer simulation of molecular^{96–111} and coarse-grained^{70,112,113} polymer models.

Previous studies have investigated depletion forces induced by aspherical depletants (e.g., rods, ellipsoids) that are fixed in size and shape.^{58,74,114,115} Recently, we explored polymer crowding and depletion forces in models of colloid–polymer mixtures, with polymers modeled as fluctuating, penetrable ellipsoids, in both *θ* solvents^{71,116,117} and good solvents,^{118} distinguished by whether polymer segments are effectively ideal (noninteracting) or nonideal, excluding volume to one another. The purpose of the present paper is to assess the interconnected influences of polymer shape and solvent quality on depletion interactions in colloid–polymer mixtures.

The remainder of this paper is organized as follows. In Sec. II, we describe modeling of linear polymer coils both as random walks and as equivalent ellipsoids that fluctuate in size and shape. The statistics governing conformational fluctuations depend on whether a coil is modeled as a self-avoiding walk, appropriate for a polymer in a good solvent, or as a non-self-avoiding random walk, corresponding to a polymer in a *θ* solvent. In Sec. III, we outline our numerical methods, based on Monte Carlo (MC) estimation of average polymer depletion volume, for computing the potential of mean force (PMF) between hard-sphere colloids induced by depletion of nonadsorbing polymers in solvents of differing qualities. Section IV presents results of our calculations of PMFs for fluctuating ellipsoidal and fixed spherical polymers in both good and *θ* solvents. Section V summarizes our study and concludes with an outlook for possible future work.

## II. MODELS

As noted above, the classic Asakura–Oosawa coarse-grained model of colloid–polymer mixtures^{1,2} idealizes nonadsorbing polymer coils as effective spheres of fixed size. The spherical polymer approximation, while incorporating an important length scale, ignores aspherical conformations and fluctuations in conformation, both of which can significantly affect depletion-induced forces. As in our earlier work on polymer crowding,^{71,116–118} we extend the AO model by representing the polymers as effective ellipsoids that fluctuate in size and shape according to random-walk statistics. In the current study, we consider only the colloid limit, in which the colloids are larger than the polymer coils, such that penetration of polymers by colloids is negligible. Although we focus here on linear homopolymers, the analysis below is easily generalized to other macromolecular architectures.^{119}

### A. Polymer coils as random walks

The size and shape of a polymer coil composed of *N* identical segments linked to form a connected chain can be characterized by a gyration tensor **T**, expressed as a matrix with elements

where *r*_{ki} denotes the *i*th component of the position vector **r**_{k} of the *k*th segment relative to the coil’s center of mass. The gyration tensor relates to the moment of the inertia tensor **I** of a rigid body via $T=Rp21\u2212I$, where **1** is the unit tensor and

is the radius of gyration of a particular coil conformation expressed in terms of the eigenvalues of **T**, Λ_{i} (*i* = 1, 2, 3). Note that *R*_{p} is invariant with respect to a change in the reference frame. The root-mean-square (rms) radius of gyration—a property of a polymer coil that is experimentally measurable (e.g., via neutron or light ray scattering)—relates to the eigenvalues via

where the angular brackets represent an ensemble average over polymer conformations.

If the ensemble average in Eq. (3) is evaluated relative to a fixed frame of reference in which the polymer coil rotates, then the average tensor describes a sphere, represented by a symmetric matrix with equal eigenvalues. On the other hand, if the average is evaluated relative to a reference frame that rotates with the principal axes of the coil and if coordinate axis labeling preserves the order of the eigenvalue magnitudes (Λ_{1} > Λ_{2} > Λ_{3}), then the average tensor describes an anisotropic object, represented by an asymmetric matrix.^{32,33} In other words, the average shape of a fluctuating random walk is spherical when viewed from a fixed frame of reference, but distinctly aspherical —shaped like an elongated, flattened bean—when viewed from the principal-axis frame.^{23,29,30}

### B. Polymer coils as fluctuating ellipsoids

The shape of a polymer coil in the principal-axis frame of reference can be fit by a general ellipsoid with principal radii proportional to the square-roots of the respective eigenvalues of the gyration tensor. In Cartesian (*x*, *y*, *z*) coordinates, the ellipsoid surface is described by

Note that the actual shape of the coil is not necessarily ellipsoidal but is approximated by an ellipsoid whose principal radii have the same proportions as those derived from the gyration tensor. For a freely jointed polymer coil of *N* segments of Kuhn length *l*, corresponding to an ideal non-self-avoiding random walk (RW),^{37} modeling conformations of a linear polymer dispersed in a *θ* solvent, the shape probability distribution is accurately fit by the analytical form^{38}

where *λ*_{i} ≡ Λ_{i}/(*Nl*^{2}) are scaled (dimensionless) eigenvalues and the fit functions are given by

with fit parameters *K*_{i}, *α*_{i}, *d*_{i}, and *n*_{i} listed in Table I. The rms radius of gyration of an ideal polymer coil is $Rg=N/6l$, while the principal radii of the ellipsoid representing a particular conformation are given by

It should be noted that the factorization ansatz of Eq. (5) is not exact since extensions of a random walk in orthogonal directions are not strictly independent. Nevertheless, conformations that significantly violate the ansatz occur only rarely for sufficiently long polymer chains.

Eigenvalue i
. | K_{i}
. | α_{i}
. | d_{i}
. | n_{i}
. |
---|---|---|---|---|

1 | 0.094 551 | 0.080 65 | 1.096 | 1/2 |

2 | 0.014 414 6 | 0.018 13 | 1.998 | 5/2 |

3 | 0.005 276 7 | 0.006 031 | 2.684 | 4 |

Eigenvalue i
. | K_{i}
. | α_{i}
. | d_{i}
. | n_{i}
. |
---|---|---|---|---|

1 | 0.094 551 | 0.080 65 | 1.096 | 1/2 |

2 | 0.014 414 6 | 0.018 13 | 1.998 | 5/2 |

3 | 0.005 276 7 | 0.006 031 | 2.684 | 4 |

For self-avoiding walks (SAWs), modeling conformations of linear polymers dispersed in a good solvent, whose segments exclude volume to one another,^{26,63,64} the rms radius of gyration is related to the segment number via *R*_{g} = *CN*^{ν}*l*, with Flory exponent *ν* = 0.588 and amplitude *C* = 0.441 08.^{35} (For an ideal polymer in a *θ* solvent, *ν* = 1/2 and $C=1/6=0.40825$.) As a common example, for polystyrene, good solvents are benzene, toluene, and chloroform, while typical *θ* solvents are cyclohexane and decalin depending on temperature.

Since the gyration tensor eigenvalues vary as *N*^{2ν}, the scaled eigenvalues are now defined as *λ*_{i} ≡ Λ_{i}/(*N*^{2ν}*l*^{2}) and are related to the principal radii via

The shape probability distribution is accurately fit by the same factorized function as for RW polymers [Eq. (5)], but with fit functions of the form

The fit parameters *a*_{i}, *b*_{i}, and *c*_{i} are tabulated in Table II. It should be emphasized that the eigenvalue distributions [Eqs. (6) and (9)] are fits to statistics from molecular simulations of linear polymer chains^{35,37,38} and reflect considerable fluctuations in polymer size and shape.

Eigenvalue i
. | a_{i}
. | b_{i}
. | c_{i}
. |
---|---|---|---|

1 | 11 847.9 | 2.355 05 | 22.3563 |

2 | 1.116 69 × 10^{9} | 3.716 98 | 148.715 |

3 | 1.068 99 × 10^{14} | 4.848 22 | 543.619 |

Eigenvalue i
. | a_{i}
. | b_{i}
. | c_{i}
. |
---|---|---|---|

1 | 11 847.9 | 2.355 05 | 22.3563 |

2 | 1.116 69 × 10^{9} | 3.716 98 | 148.715 |

3 | 1.068 99 × 10^{14} | 4.848 22 | 543.619 |

The average shape of a polymer coil can be quantified by an asphericity parameter,^{32,33}

defined such that a perfect sphere, with all eigenvalues equal, has $A=0$, while a needle-like object has $A\u22431$. Interestingly, for both RW and SAW coils, $A\u22430.54$ in uncrowded environments.^{118}

The probability distributions *P*_{i}(*λ*_{i}) for the individual eigenvalues differ somewhat from the fit functions *f*_{i}(*λ*_{i}) in Eqs. (6) and (9). Each is obtained from the parent distribution [Eq. (5)] by integrating over the other two eigenvalues, with limits imposed by the requirement of eigenvalue ordering (*λ*_{1} > *λ*_{2} > *λ*_{3}),

For comparison, Fig. 1 shows the scaled eigenvalue distributions of polymers in *θ* and good solvents. Note that, accounting for the different scaling factors—*N* for RW polymers, but *N*^{1.176} for SAW polymers—the unscaled eigenvalues are significantly larger for a SAW polymer in a good solvent than for a RW polymer in a *θ* solvent, reflecting the more extended conformations of polymers with excluded-volume interactions.

### C. Colloid–polymer mixtures

To explore the influence of aspherical polymer conformations and solvent quality on the effective interactions induced between colloidal particles by depletion of nonadsorbing polymers due to colloid excluded volume, we consider a monodisperse suspension of colloidal particles, modeled as hard spheres of radius *R*_{c}, mixed with free polymer coils, modeled as ellipsoids whose shapes fluctuate according to the statistics of random walks (Fig. 2).

The strength and range of depletion-induced interactions depend on concentration and size of the polymer coils relative to the colloids. Calibrating theoretical models to experimental systems requires an appropriate measure for the effective size of a polymer coil. Implementations of the AO model often take the effective radius of a polymer coil simply as the rms radius of gyration, defining the polymer-to-colloid size ratio as *q* ≡ *R*_{g}/*R*_{c}. For later reference, we note that, for a given size ratio *q*_{RW} of a RW polymer in a *θ* solvent, the scaling relations (Sec. II B) dictate the size ratio *q*_{SAW} of a SAW polymer of equal molecular weight (same *N*) in a good solvent,

More accurate measures for the effective size of a polymer account for the effect on the depletion layer thickness of deformation of a coil near a hard surface. We first review the simplest case of a polymer near a hard, flat wall, as addressed by Asakura and Oosawa,^{1} and then consider polymers near hard, spherical colloids.

Consider a solution of *N*_{p} polymers in a volume V containing two hard, flat, parallel plates of area *A* separated by a distance *D* much shorter than the lateral extent of the plates ($D\u226aA$). The polymer coils are free to diffuse, except for the constraint imposed by the plates. In a closed system, the potential of mean force between the plates induced by depletion of polymers from the intervening space is defined as the difference in the Helmholtz free energy *F*(*D*) at plate separation *D* and at infinite separation,

In the dilute limit of noninteracting polymers amid plates separated by distance *D*, the free energy is related to the single-polymer partition function $Z1(D)$ via

The potential of mean force between plates induced by depletion of *N*_{p} polymer coils can then be expressed as

Considering that the polymers and plates interact only via excluded-volume interactions, the single-polymer partition function is simply proportional to the fraction of volume available to a polymer,

where $x0$ and $x(D)$ represent average thicknesses of depletion layers adjacent to a surface outside of and between the plates, respectively, and angular brackets denote ensemble averages over depletant conformations (see Fig. 3). For depletants with simple geometrical shapes (e.g., spheres or ellipsoids), it can be shown that the depletion layer thickness equals the integrated mean curvature *c*_{d} of the depletant.^{120} For a smooth, convex body, *c*_{d} is defined as an average over the body’s closed surface *S* of the mean curvature,

where *R*_{1}(**r**) and *R*_{2}(**r**) are the local radii of curvature at a point **r** on the surface. To model depletants that fluctuate in size and shape, we augment this definition with an average over conformations,

Note that *c*_{d} has physical dimensions of length.

From Eq. (18) and the limiting relation

the partition function for a polymer in a system with infinitely separated plates is

Substituting Eqs. (18) and (22) into Eq. (17) and assuming a polymer solution so dilute that the excluded volume is only a small fraction of the total volume, the PMF between the plates per unit plate area is

where Π_{p} = *N*_{p}*k*_{B}*T*/*V* is the osmotic pressure of an ideal gas of polymer coils. In the limit as the plates come together (*D* → 0), where $x(D)\u21920$, the PMF induced by real polymer chains with the radius of gyration *R*_{g} approaches the exact contact value^{1,5,64}

Now, identifying $x0$ with the integrated mean curvature *c*_{d} (averaged over conformations) of an uncrowded polymer coil, modeled as a fluctuating ellipsoid (see Fig. 3), and defining *γ* ≡ *c*_{d}/*R*_{p,eff} as the coefficient of proportionality between *c*_{d} and the *effective* polymer radius *R*_{p,eff}, we have

and thus, finally,

We emphasize that *R*_{p,eff} represents the effective radius of a nonadsorbing polymer coil in the presence of hard colloidal particles, taking into account both the diffuse periphery and the aspherical shape of the coil. We conclude that an experimental system with the nominal polymer-to-colloid size ratio *q* = *R*_{g}/*R*_{c}, defined by reference to the radius of gyration of a real polymer chain, should be modeled using an *effective* size ratio

For a sphere of fixed radius (AO model), the integrated mean curvature simply equals the radius (*γ* = 1). For a fluctuating ellipsoid, on the other hand, determining *γ* is nontrivial. We computed *γ* using two independent, but equivalent, numerical methods. In the first method, based on Eq. (25), we numerically integrated the mean curvature over the surface of the ellipsoid and averaged over the polymer shape probability distribution [Eq. (5) combined with Eqs. (6) or (9) for RW or SAW polymers, respectively]. In the second method, we computed the half-width of an ellipsoid in a fixed direction and numerically averaged over orientations and the shape probability distribution. Both methods involve numerically evaluating a five-dimensional integral, and both give the same result to within numerical precision. From the shape distributions corresponding to a RW polymer (*θ* solvent) and a SAW polymer (good solvent), we find *γ*_{RW} = 0.932 54 and *γ*_{SAW} = 0.924 31, respectively.

For polymers dispersed in a suspension of hard-sphere colloids, the effective size ratio *q*_{eff} can be objectively defined by equating the free energy cost of inserting a hard sphere into a solution of polymers, as predicted by polymer field theory, with the work done to inflate a sphere in the model polymer solution. When applied to polymers obeying RW statistics and presumed spherical in shape and fixed in size (AO model), this definition yields^{75,100,121,122}

assuming that *q*_{eff} < 1, such that penetration of a polymer by a colloid can be neglected. For polymers with aspherical (ellipsoidal), fluctuating shapes dispersed with hard-sphere colloids, we modify this definition in the same manner as for polymers near a hard, flat wall by incorporating the integrated mean curvature,

This definition ensures that in the limit *q* → 0, the model recovers the exact contact value of the PMF induced by RW polymers between hard, flat plates. We emphasize that the 1/*γ* adjustment proved essential in our earlier study^{117} for achieving quantitative agreement with PMF data from molecular simulations^{96,97} and from experiments on DNA-induced depletion forces.^{57}

In contrast, for polymers that obey SAW statistics, also presumed spherical and of fixed size (AO model), field theory yields an effective size ratio^{36,122}

where *C*_{1} ≃ 3.2130, *C*_{2} ≃ 2.6073, and *C*_{3} ≃ 0.1197. For aspherical (ellipsoidal) SAW polymers, we similarly incorporate the integrated mean curvature and define

### D. Potential of mean force

Adapting Eq. (15) from plates to spheres, the potential of mean force between two colloids in thermal equilibrium at absolute temperature *T* with a solution of nonadsorbing polymers in a closed volume is defined as the change in Helmholtz free energy *F*(*r*) of the system upon bringing the particles from infinite separation to center-to-center separation *r*,

since in an isotropic fluid, the pair potential depends on only the radial coordinate. For a system in chemical equilibrium with a polymer reservoir, the PMF is defined as the change in grand potential Ω(*r*). In earlier work,^{116} we applied an alternative (but equivalent) definition, which is more appropriate in the nanoparticle limit, in which polymer coils are significantly larger than and penetrable by the colloids (nanoparticles).

The free energy varies with colloidal separation due to mechanical work performed by the colloids in changing the excluded volume of the polymer with osmotic pressure Π_{p} = *n*_{p}*k*_{B}*T*, assuming a dilute (ideal gas) polymer solution of mean density *n*_{p}. In the AO model, this work is easily evaluated,

where *A*_{ov}(*r*) and *V*_{ov}(*r*) are the cross-sectional area and volume, respectively, of the overlap region of the two excluded-volume shells, and we choose *F*(*∞*) = 0. It should be noted that, when used to model thermodynamic phase behavior, Eq. (33) must be corrected for triplet overlaps at size ratios above *q* ≃ 0.1547, especially away from the dilute colloid concentration limit. For spherical colloids and spherical polymers of fixed radius *R*_{p}, the convex-lens-shaped pair overlap region, defined by the intersection of two spherical excluded-volume shells, has volume

for 2*R*_{c} < *r* < 2(*R*_{c} + *R*_{p}) (otherwise zero). Equations (33)–(34) express the conventional AO potential.

In the case of aspherical depletants, this simple geometric approach can be adapted by calculating an *average* of the overlap volume $Vov(r)$ over an ensemble of polymer conformations (orientations and shapes). From a large sample of randomly generated conformations (microstates), $Vov(r)$ equals the sampled volume times the fraction of microstates in which a depletant overlaps both colloids.

Two limitations of our modeling approach are important to note. First, the coarse-grained model of polymer coils necessarily neglects coil shapes that deviate from ellipsoidal, which may affect how nonadsorbing polymers interact with hard colloidal surfaces. Second, although neighboring coils in a *θ* solvent do not influence each other’s shape distribution, since polymer segments are effectively noninteracting, excluded-volume interactions between segments in a good solvent can lead to correlations between shapes of neighboring coils. Such correlations may affect the strength and range of polymer depletion-induced interactions, especially in semi-dilute or concentrated polymer solutions. The present modeling approach, which describes only the PMF induced by independent ellipsoidal polymer coils, neglects such effects. In Sec. IV, we discuss implications and potential remedies of these limitations of our approach.

## III. MONTE CARLO SIMULATIONS

To compute the potential of mean force between colloids induced by depletion of nonadsorbing polymers that fluctuate in size and shape according to either RW chain statistics (*θ* solvent) or SAW statistics (good solvent), we used Monte Carlo (MC) simulation methods.^{123} Applying Eq. (33), we determined the average overlap volume $Vov(r)$ by placing two hard-sphere colloids in a rectangular parallelepiped simulation box at center-to-center separation *r*, inserting a polymer ellipsoid at a random position with random orientation and shape governed by the appropriate gyration tensor eigenvalue probability distribution [Eqs. (6) or (9)], and counting the fraction of double overlaps, i.e., insertions leading to an overlap of the ellipsoid with both spheres.

As noted in Sec. II, since our model constrains polymer coils to have only ellipsoidal shapes, our approach, although it captures the gross shapes of polymers, neglects any influence of non-ellipsoidal conformations on polymer–colloid interactions. Furthermore, since we insert polymer coils only one at a time, our approach, when applied to SAW polymers in good solvents, is limited to dilute polymer solutions since it neglects possible correlations between shapes of neighboring, interacting coils.

To randomly sample polymer conformations, we implemented a variation of the Metropolis algorithm.^{123} Trial changes in orientation and shape of a polymer ellipsoid were coupled with insertions. To uniformly sample orientations, specified by a unit vector **u** aligned with the longest axis of the ellipsoid, we generated a new (trial) unit vector, **u**′ = (**u** + *τ***v**)/|**u** + *τ***v**|, where **v** is a randomly oriented unit vector and *τ* is a tolerance.^{123} A trial change in the shape from one set of gyration tensor eigenvalues *λ* ≡ {*λ*_{i}} to a new set *λ*′ ≡ {*λ*_{i}′ = *λ*_{i} + Δ*λ*_{i}} with tolerances Δ*λ*_{i} (*i* = 1, …, 3) implies a change Δ*F*_{c} in the coil’s internal free energy,^{64} *F*_{c} = −*k*_{B}*T* ln *P*(*λ*), where *P*(*λ*) is the polymer shape distribution [Eq. (5) with Eqs. (6) or (9)]. A trial conformation was rejected if the inserted polymer ellipsoid overlapped either colloidal sphere. Otherwise, it was accepted with probability

If the trial conformation was accepted, the ellipsoid’s orientation and shape were updated and the double-overlap counter was incremented. Limiting our study to dilute solutions, we inserted polymers one at a time, thus neglecting polymer–polymer interactions. To diagnose the overlap of a colloid and a polymer, we computed the shortest distance between a point (sphere center) and the ellipsoid surface, requiring evaluating the roots of a sixth-order polynomial.^{124} This sampling method yields the average volume of the polymer depletion region surrounding two colloidal spheres and hence, from Eq. (33), the PMF.

## IV. RESULTS AND DISCUSSION

To compare potentials of mean force between colloidal hard spheres induced by depletion of nonadsorbing polymers in *θ* and good solvents, we implemented the ellipsoidal polymer model described in Sec. II and performed a series of Monte Carlo simulations. Using the polymer trial insertion method outlined in Sec. III, we computed the PMF over a range of colloid separations. The side lengths of the rectangular parallelepiped simulation box were set small enough to maximize the acceptance ratio, while large enough to avoid the interaction of a polymer with periodic images of the colloids. Tolerances for polymer trial moves, optimized by trial and error, were fixed at *τ* = 0.001 for rotations and Δ*λ*_{1} = 0.01, Δ*λ*_{2} = 0.003, and Δ*λ*_{3} = 0.001 for shape changes. For a given colloid pair separation, we performed 10^{6} independent trial polymer insertions and then computed statistical uncertainties (error bars) as standard deviations from five independent runs for a total of 5 × 10^{6} trial insertions.

To validate our methods, we first implemented the original AO model of spherical polymers of fixed size and confirmed that our algorithm reproduces the exact PMF predicted by Eqs. (33) and (34). We then proceeded to simulate the ellipsoidal polymer model for polymers whose sizes and shapes are governed by RW and SAW chain statistics, corresponding to polymers in *θ* and good solvents, respectively. To compare depletion of polymers of equal segment number in different solvents, we converted polymer-to-colloid size ratios between RW and SAW statistics using Eq. (14). This conversion requires specifying the ratio of the colloid radius *R*_{c} to the polymer segment length *l*. To make potential contact with experiments, we chose typical values of *R*_{c} = 100 nm and *l* = 0.76 nm, corresponding to polyethylene glycol (PEG) in water.^{125}

The results of our calculations for the PMF are presented in Fig. 4 over a range of size ratios. Note that the vertical axis is scaled by the polymer osmotic pressure Π_{p}, rendering the plotted PMF independent of polymer concentration. As discussed above in Secs. II and III, however, the results shown for SAW polymers are physically meaningful only for dilute polymer solutions. For each bare size ratio of RW polymers in the series *q*_{RW} = 0.1, 0.2, 0.3, 0.4, we calculated the corresponding bare size ratio of SAW polymers *q*_{SAW} of equal *N* from Eq. (14). We then ran simulations for the effective size ratios of ellipsoidal RW and SAW polymers, *q*_{eff,RW} and *q*_{eff,SAW}, calculated from Eqs. (29) and (31), respectively. For comparison, we also simulated spherical SAW polymers (AO model) with the effective size ratio *q*_{eff,SAW-AO}, calculated from Eq. (30). The system parameters are tabulated in Table III.

q_{RW}
. | q_{eff,RW}
. | q_{SAW}
. | q_{eff,SAW}
. | q_{eff,SAW-AO}
. | N
. |
---|---|---|---|---|---|

0.1 | 0.118 21 | 0.199 09 | 0.219 89 | 0.203 25 | 1040 |

0.2 | 0.231 38 | 0.449 84 | 0.471 86 | 0.436 14 | 4155 |

0.3 | 0.340 18 | 0.724 67 | 0.724 16 | 0.669 35 | 9350 |

0.4 | 0.445 18 | 1.016 41 | 0.971 16 | 0.897 65 | 16620 |

q_{RW}
. | q_{eff,RW}
. | q_{SAW}
. | q_{eff,SAW}
. | q_{eff,SAW-AO}
. | N
. |
---|---|---|---|---|---|

0.1 | 0.118 21 | 0.199 09 | 0.219 89 | 0.203 25 | 1040 |

0.2 | 0.231 38 | 0.449 84 | 0.471 86 | 0.436 14 | 4155 |

0.3 | 0.340 18 | 0.724 67 | 0.724 16 | 0.669 35 | 9350 |

0.4 | 0.445 18 | 1.016 41 | 0.971 16 | 0.897 65 | 16620 |

In Fig. 4, the solid symbols represent PMF data for polymers in a good solvent, modeled as fluctuating ellipsoids obeying SAW chain statistics [Eq. (6)]. The lightly shaded symbols represent PMF data for polymers in a *θ* solvent, modeled as fluctuating ellipsoids obeying RW chain statistics [Eq. (9)]. For comparison, the open symbols represent PMF data for polymers in a good solvent, modeled as spheres of fixed size (AO model). Since in the dilute limit, the PMF is proportional to the polymer osmotic pressure, we plot the dimensionless quantity of *v*_{mf}(*r*) scaled by Π_{p} times the cube of the colloid diameter *σ*_{c} = 2*R*_{c}. Note that negative values of *v*_{mf}(*r*) imply an attractive pair interaction between colloids. To within statistical uncertainty, our data for the AO model are perfectly fit by the analytical expressions of Eqs. (33) and (34) (curves).

Having validated our methods, we now compare our PMF data for polymer coils modeled as fluctuating ellipsoids (solid symbols in Fig. 4) with corresponding data for coils of the equal segment number modeled as spheres of fixed size (open symbols in Fig. 4), both obeying SAW chain statistics. Evidently, fluctuating ellipsoid polymers induce generally weaker PMFs than fixed-sphere polymers (AO model). For ellipsoidal polymers, the contact value *v*_{mf}(*σ*_{c}) is consistently lesser in magnitude, while the range of *v*_{mf}(*r*) is consistently longer than for spherical polymers. Furthermore, deviations between the PMFs from the ellipsoidal and spherical polymer models grow with the increasing size ratio.

These results may appear surprising, considering that depletion-induced attraction in the AO model strengthens with increasing size ratio and given that the effective size ratio in the ellipsoidal polymer model [Eq. (31)] exceeds that in the spherical polymer model [Eq. (30)] (the integrated mean curvature being smaller for an ellipsoid than for a sphere of equal radius of gyration). Nevertheless, these trends are quite consistent with the extra conformational freedom of fluctuating ellipsoids to elongate to lengths beyond their mean diameter and to deform to avoid hard surfaces. In previous work,^{117} we showed that the fluctuating ellipsoid polymer model, when implemented with the appropriate effective size ratio, nearly exactly reproduces the PMF computed from “lattice–polymer” simulations of RW polymers whose segments are confined to the sites of a cubic lattice.^{96,97}

Next, we compare our results for PMFs induced by depletion of polymer coils, modeled as fluctuating ellipsoids, that obey either SAW or RW chain statistics. From Fig. 4, we see that SAW polymers in a good solvent (solid symbols) induce PMFs that are significantly stronger—both greater in magnitude and longer in range—than RW polymers of the equal segment number in a *θ* solvent (lightly shaded symbols). Qualitatively, this trend is consistent with the more extended conformations and correspondingly larger effective radius of gyration of SAW polymers compared with RW polymers of the same contour length. Quantitatively, it is interesting that, at least for the system parameters considered here, the attractive well of the PMF induced by SAW polymers is roughly three times deeper than that of the PMF induced by RW polymers.

Figure 5 collects data from Fig. 4 in one plot to summarize the dependence on the polymer-to-colloid size ratio of the PMF induced by SAW polymers. With the increasing effective size ratio, the depth and range of *v*_{mf}(*r*) both steadily grow. As noted above, the fluctuating ellipsoid polymer model predicts a generally weaker, but longer-ranged, PMF compared with the AO model. Although we are not aware of lattice–polymer simulations of colloids dispersed in SAW polymer solutions, our predictions could be tested against such molecular-scale simulations.

As mentioned in Secs. II and III, our modeling approach is based on two main approximations. First, it neglects non-ellipsoidal shapes of nonadsorbing polymer coils interacting with colloidal surfaces, and second, it neglects interactions and associated shape correlations between neighboring coils. The accuracy of the first approximation can be quantified by comparing predictions of the coarse-grained polymer model with those of molecular-scale polymer models. In our previous study of depletion interactions induced by polymers in a *θ* solvent,^{117} we compared our PMF results against simulations of a model of RW polymers on a lattice, demonstrating remarkable accuracy of the coarse-grained polymer model. A similar comparison for SAW polymers in a good solvent could further assess the accuracy of the coarse-grained model.

The second approximation is well justified for polymers in *θ* solvents, except to the extent that the *θ* temperature may vary with polymer concentration. For SAW polymers in good solvents, however, neglecting interactions between neighboring coils strictly limits application of our approach to dilute polymer solutions. Consequently, the model may not accurately describe the influence of depletion forces on thermodynamic properties, including bulk phase separation in concentrated mixtures of colloids and polymers in good solvents. Extending the model beyond the dilute regime would require incorporating the influence on polymer conformations of interactions and correlations between segments of different coils. This extension could be achieved, for example, by simulating a concentrated solution of a few explicit SAW interacting polymers, computing their gyration tensors, and fitting the eigenvalue probability distributions over a range of polymer concentrations.

## V. SUMMARY AND CONCLUSIONS

In summary, we have implemented Monte Carlo simulation methods for computing the potential of mean force between hard-sphere colloids induced by nonadsorbing polymer coils dispersed in good and *θ* solvents. For computational efficiency, we modeled the polymer coils as general ellipsoids whose conformations (size and shape) fluctuate according to statistics of either self-avoiding walks (good solvent) or non-self-avoiding random walks (*θ* solvent). The principal radii of the equivalent ellipsoid representing a polymer coil were determined from accurate fits to probability distributions for the eigenvalues of the gyration tensor of a SAW or RW.

In the colloid limit in which the polymer radius of gyration is smaller than the colloid radius, we determined the PMF by computing the average volume of the polymer depletion region surrounding a pair of colloids using a polymer insertion algorithm. Because polymer conformational distributions vary with solvent quality, the average depletion volume and therefore the PMF differ between good and *θ* solvents. Our results demonstrate that the dependence of the PMF on polymer shape and solvent quality can be quite significant.

Comparing the ellipsoidal and spherical (AO) models of polymers in good solvents, we showed that the former model yields a generally weaker PMF, which we attribute to the freedom allowed by the fluctuating ellipsoid model for a polymer coil to adapt its shape to a crowded environment. This finding is consistent with conclusions from our previous studies of depletion interactions in colloid–polymer mixtures dispersed in *θ* solvents.^{116,117}

Comparing the ellipsoidal model of SAW and RW polymers in good and *θ* solvents, respectively, we showed that polymers in good solvents induce considerably stronger PMFs than polymers of the same number of segments in *θ* solvents. This trend is explained by the more extended conformations of SAW polymers, which tend to enlarge the depletion region around colloids. The coarse-grained model of depletion interactions induced by nonadsorbing polymers could be further tested against experiments and molecular-scale models of polymers, especially for SAW coils dispersed in good solvents.

Our work suggests the possibility of tuning effective depletion-induced interactions between colloids, and thereby thermodynamic phase behavior of colloid–polymer mixtures, by varying solvent quality, e.g., by selecting particular polymer–solvent combinations or by changing temperature and cosolvent concentration for a given combination. Our modeling approach, which focuses on the geometry of the depletant, may also have relevance for effective interactions induced by other types of soft depletants that can vary in size and shape, such as vesicles^{126} and microgels.^{127,128}

Future work could explore the influence of solvent quality on depletion-induced interactions in the protein limit, in which polymers are large enough to be penetrated by colloids. As shown in our previous studies of polymer–nanoparticle mixtures,^{116,118} accurate modeling would require a reliable expression (from polymer field theory) for the penetration free energy. In this regime, our geometric approach may yield insights complementary to those provided by field theories,^{75–82} integral-equation theories,^{83–87} and density-functional theories.^{88–92} The model could also be extended to concentrated polymer solutions by incorporating interactions and correlations between segments within different coils and correspondingly modifying the polymer shape distributions, accounting for the possible role of polymer density fluctuations near polymer–solvent demixing critical points.^{85} Finally, our approach could be generalized to model depletion interactions between aspherical hard colloids,^{129} such as rods or platelets.

## ACKNOWLEDGMENTS

This work was supported by the National Science Foundation (Grant No. DMR-1928073). We thank Wei Kang Lim for important contributions to the coding of the simulations and Sergio J. Sciutto for helpful correspondence regarding the SAW shape distribution.

## DATA AVAILABILITY

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