Depletion Interactions between Nanoparticles: The Effect of the Polymeric Depletant Stiffness

A Density Functional Theory is employed to study depletion interactions between nanoparticles mediated by semiflexible polymers. The four key parameters are the chain contour length and the persistence length of the polymeric depletant, its radius of gyration, and the nanoparticle radius. In the Density Functional Theory calculation of the depletion interaction between the nanoparticles mediated by semiflexible polymers, the polymer gyration radius is kept constant by varying the contour length and the persistence length simultaneously. This makes it possible to study the effect of the chain stiffness on the depletion potential of mean force between the nanoparticles for a given depletant size. It is found that the depletion attraction becomes stronger for stiffer polymer chains and larger colloids. The depletion potential of mean force is used as input to compute the phase diagram for an effective one-component colloidal system.


Introduction
Spherical colloidal particles dispersed in a solution of nonadsorbing polymers play an important role in various systems of biological and technological interest. Accordingly, a detailed understanding of the phase behavior and stability of such dispersions is of primary importance for the development of various industrial applications. As a result, colloidal dispersions in polymer solutions have been actively studied both experimentally [1,2] and theoretically [3][4][5][6]. Despite these extensive studies, several aspects of this problem have remained relatively unexplored. First, most theoretical studies have been limited to fully flexible polymers. At the same time, all realistic polymeric systems involved in numerous practical applications are characterized by some degree of the chain stiffness. Accordingly, the first important research question to be addressed in the present study concerns the effect of the polymer stiffness on the depletion interaction between the colloids. Second, to the best of our knowledge, a direct link between microscopic depletion interaction and phase diagram has not been established yet for semiflexible depletants. Hence, the second research question deals with establishing such a connection by using the depletion interaction calculated from a microscopic model as input for computing a comprehensive colloid-polymer phase diagram. The importance of these two related research questions stems from the fact that a reliable theoretical approach for studying realistic colloid-polymer mixtures and their phase behavior would be of great utility for analyzing the existing experimental data and making suggestions for future experiments. In order to further highlight the novel aspects of the present work and put it in broader context, we next briefly review the existing literature in this field.
From a microscopic perspective, the polymer-induced depletion attraction between two colloidal particles approaching each other in a polymer solution arises from the fact that the chains are expelled from the region between the two colloidal spheres into the bulk due to the loss of configurational entropy by the chains in the region between the spheres. The resulting unbalanced pressure exerted by the polymer chains on the outward surfaces of the two colloids produces an effective depletion attraction between the spheres, which is quantified by the polymer-mediated potential of mean force (PMF) between the two colloids. The knowledge of this PMF makes it possible to treat the colloidal dispersion as an effective one-component system and to construct the corresponding phase diagram describing the stability of the dispersion against flocculation [7,8].
The early pioneering work of Asakura and Oosawa [9,10] developed a theory for the polymer-mediated PMF in colloid-polymer suspensions by treating colloids as hard spheres and polymers as spheres that are fully penetrable to each other but not to the colloids. The resulting PMF produced by the Asakura-Oosawa model is always attractive, with the strength of attraction increasing monotonically with increasing polymer concentration. Following this early work, a large variety of theoretical approaches have been developed to treat colloidal interactions in polymer solutions, including scaling arguments [11], selfconsistent field theory [12], and the adsorption method [4] based on the superposition approximation of one-particle depletion layers. In addition, several authors have used the polymer reference interaction site model [13] within the framework of the integral equation theory to compute the PMF in colloid-polymer solutions [14][15][16][17][18][19].
As already mentioned, the vast majority of the aforementioned theoretical studies are limited to fully flexible polymer chains. At the same time, numerous polymeric molecules employed in various practical applications are characterized by a certain degree of stiffness, which is quantified by the chain persistence length [20]. The latter parameter strongly affects the thickness of the polymer depletion layer around a colloid [21], and therefore the polymer-mediated PMF between the colloids [21]. However, both experimental [22,23] and theoretical [16,21,24,25] studies of the chain stiffness effects on the polymer-mediated PMF are still relatively scarce. Among theoretical methods applied to this problem, one can mention integral equation theory [16], self-consistent field theory [21], and density functional theory (DFT) [24,25]. While the former two methods were applied to spherical colloids, the latter approach was limited to the studies of polymer-mediated interactions between flat walls. Accordingly, it would be of interest to apply the DFT formalism to compute the PMF between spherical colloids mediated by semiflexible polymers. The first goal of the present study is to perform such a calculation, with a particular focus on the case when the polymer radius of gyration is comparable to the colloid radius.
As mentioned earlier, the fundamental importance of the polymer-mediated PMF is due to the fact that it allows one to map the colloid-polymer binary mixture onto an effective one-component colloidal system. Once this goal is accomplished, the phase diagram of this one-component system can be computed using standard methods [3]. In the earlier self-consistent field theory study [21] of a binary mixture of spherical colloids and semiflexible polymers, the polymer-induced PMF was employed to compute the second virial coefficient for the colloids, which can be used to characterize the stability of the colloidal dispersion. However, the comprehensive phase diagram of the dispersion has not been obtained directly from he PMF, but rather from the Free Volume Theory [26]. Accordingly, the second goal of the present work is to compute the phase diagram directly from the PMF generated by the DFT approach.
The outline of the remainder of the paper is as follows. In Section 2, we specify our microscopic model; in Section 3, we outline the DFT formalism employed to calculate the polymer-mediated PMF between colloids. In Section 4, we present our approach to obtain the phase diagrams. Section 5 presents our results, and Section 6 concludes the paper.

Microscopic Model
We consider hard-sphere colloidal particles with radius R c embedded in a solution of semiflexible polymer chains composed of N tangent hard sphere beads with diameter σ, i.e., all the bond lengths are fixed at l b = σ (σ will be used as the length unit throughout this work). In order to study the effect of chain flexibility on the brush structural properties, we employ a bond-bending potential [27,28]: where θ ijk is the bond angle formed between the two subsequent vectors s i and s i+1 along the bonds connecting monomers i, j = i + 1 and j, k = i + 2, i.e., s i = r i+1 − r i and s i+1 = r i+2 − r i+1 . The energy parameter b then controls the persistence length l p , which is defined as [29] For semiflexible chains with b ≥ 2, one has the persistence length l p /l b ≈ β b = κ, where β = 1/k B T (T is the temperature), and κ is the dimensionless stiffness parameter [27,28].
Both colloid-colloid (v cc (r)) and monomer-monomer (v mm (r)) excluded volume interactions are of the hard-sphere type: v cc (r) = ∞ if r < 2R c and zero otherwise, and v mm (r) = ∞ if r < σ and zero otherwise. The monomer-colloid interaction is modeled as a sum of hard-sphere repulsion at contact (v mc (r) = ∞ if r < R mc = 0.5σ + R c ) and a soft Gaussian repulsion at larger separations: βv where R g is the polymer chain gyration radius. In order to isolate the effect of the chain stiffness parameter κ on the polymer-induced depletion PMF between the two colloids while keeping the depletant size fixed, we simultaneously vary the chain contour length and its persistence length in such a way that R g remains constant.

Density Functional Theory
The major goal of the present work is to compute the phase diagram of an effective one-component colloidal system by tracing out the polymeric component. In order to achieve this goal, one needs to compute the polymer-mediated PMF W(R) between the two colloids separated by distance R. Combining W(R) with the bare colloid-colloid potential v cc (R) gives the total interaction V cc (R) between two colloids: In order to obtain W(R), we define ρ(r, R) as the conditional probability of finding a polymer bead at r given that one colloid is at the origin and the other is located at R (R = |R|). With this definition, the polymer-mediated PMF between the two colloids is given by the following exact relations [14]: where the outwards excess mean force, F(R), is given by: whereR is the unit vector along the line connecting the two colloids, and ρ(r, R) is the anisotropic monomer density profile induced by the two colloids. In the present study, we construct the latter density profile on the basis of the Kirkwood superposition approximation (KSA) [8], whereby ρ(r, R) is approximated by the product of spherically symmetric density profiles ρ(r) around individual colloids: ρ KSA (r, R) ≈ ρ(r)ρ(|r − R|)/ρ, with ρ being the bulk monomer density. The isotropic monomer density profile ρ(r) around a single colloid is obtained from the DFT formalism [30][31][32]. As a starting point of the DFT-based approach, one writes an expression of the grand free energy, Ω, as a functional of the polymer density profile ρ p (R p ), where R p = (r 1 , r 2 , · · · , r N ) is a collective variable with the individual monomer coordinates r i . The average monomer density ρ(r) is related to the molecular density profile, ρ p (R p ), as follows: The minimization of Ω with respect to ρ p (R p ) yields the equilibrium polymer density distribution. The functional Ω is related to the Helmholtz free energy functional, F, via a Legendre transform: where µ is the polymer chemical potential, and V ext (R p ) is the external field, which in the present case is due to the interaction of the polymer beads with the colloid: We employ the following approximation for the Helmholtz free energy functional, which separates it into ideal and excess parts according to [33]: with the ideal functional given by [34,35]: where V bend is given by Equation (1), and V b (R p ) is the binding energy given by [36]: For the excess free energy functional, we adopt the weighted density approximation [37]: with the weighted density given by: In the above, the monomer density ρ(r) is given by Equation (6), and f ex (ρ) is the excess free energy density per site of the polymer solution with site density ρ arising from the short-ranged hard-core repulsive interactions. We compute it from the Wertheim's expression which was obtained on the basis of the first-order thermodynamic perturbation theory [38]: where η = πσ 3 ρ/6 is the monomer packing fraction.
In the present work, we employ the simple square-well form for the weighting function w(r), whose range is given by the diameter σ of the polymer segment [39]: where Θ(r) is the Heaviside step function. While more sophisticated forms of the weighting function are available in the literature (e.g., those used in the Fundamental Measure Theory version of DFT [40]), earlier studies [41] have shown relative insensitivity of DFT results for polymeric systems to the specific choice of the weight function.
The minimization of the grand free energy functional Ω yields the following result for the equilibrium polymer density profile [32]: where λ(r) = β δF ex δρ(r) + βv mc (r).
Substitution of ρ p (R p ) into Equation (6) then yields an integral equation for the monomer density distribution ρ(r) which needs to be solved numerically [32].
Regarding the numerical implementation of the DFT procedure, given that the monomer density distribution around a single colloid is spherically symmetric, the corresponding integral equation for ρ(r) is solved numerically on an equidistant grid along the radial coordinate r with the grid spacing ∆r = 0.02. A simple Picard iteration procedure was employed [42], and tolerance criterion for terminating the iterative procedure was set to 10 −6 .

Phase Diagrams
In order to calculate the phase diagram, we follow the standard approach [3,8] and map the two-component colloid-polymer mixture onto an effective one-component colloidal system by tracing out the polymeric component. This goal is accomplished via Equation (3), whereby the bare hard-sphere colloid-colloid potential v cc (r) is augmented by the polymermediated depletion PMF W(R) obtained from Equation (4) using DFT formalism outlined in Section 3. In calculating the phase diagram, we consider both fluid (vapor and liquid) and solid phases. The phase boundaries are obtained by equating the pressure and the colloid chemical potential in the two coexisting phases [3,8].
The dimensionless pressure of a fluid phase is given by [3]: where V c = 4πR 3 c /3 is the volume of the colloidal sphere, ρ c is the colloid number density, and η c = ρ c V c is the colloid packing fraction. The first term on the right-hand side of Equation (18) originates from the Carnahan-Starling equation of state [43], while the second term comes from the mean-field treatment of the effective colloid-colloid attraction due to the polymer-mediated depletion interaction.
The dimensionless chemical potential of a fluid phase is given by [3]: where Λ c is the de Broglie wavelength of the colloidal particle. The dimensionless pressure of a solid phase is given by [3]: where η cp = π/(3 √ 2) is the value of η c at close packing. Finally, the dimensionless chemical potential of a solid phase is given by [3]:

Results
In order to focus on the effect of the polymer stiffness parameter κ on the colloidcolloid polymer-induced depletion interaction and the corresponding phase behavior of the effective one-component colloidal system, we start by computing the polymer chain radius of gyration R g for a range of values of the chain contour length N and stiffness parameter κ. The radius of gyration for our microscopic model of the semiflexible chain specified in Section 2 is calculated using the DFT methodology described in detail in Ref. [44] (with the Helmholtz free energy functional and the weighted density defined in Section 3). All the results reported below were obtained for the value R g = 10 and the following five pairs of values of N and κ: (N = 40, κ = 25.1), (N = 48, κ = 11.5), (N = 64, κ = 6.2), (N = 80, κ = 4.4), and (N = 96, κ = 3.5). While in the first pair the contour and persistence lengths are comparable, in the last pair l p is nearly two orders of magnitude smaller than contour length. Thus, these selected values span the range from semiflexible to nearly fully flexible polymeric depletants (at the fixed depletant size).
Using the DFT approach, we compute the polymer-induced PMF for five pairs of values of N and κ listed above; the monomer bulk density is fixed at ρ = 0.001. The DFT results for the dimensionless PMF βW(R) are presented in Figure 1 for two values of the colloid radius: R c = 5 in the upper panel, and R c = 20 in the lower panel. Thus, the former case corresponds to the situation R c < R g , while in the latter case has R c > R g . One sees that for both values of the colloid radius the strength of the depletion attraction and its range increase with increasing chain stiffness, in agreement with earlier SCF results [21]. The depletion attraction also becomes stronger with increasing colloid radius, as one would expect [21]. In order to compute the phase diagrams of the effective one-component colloidal systems, we obtain the depletion PMFs for a range of values of the reservoir monomer density ρ r [26] and then use these results to calculate the phase boundaries as described in Section 4. Some representative results in the variables η c -ρ r (reservoir representation [26]) are shown in Figure 2, where the upper panel corresponds to the colloid radius R c = 5, and the lower panel corresponds to the colloid radius R c = 20. In the upper panel, the solid lines correspond to the semiflexible polymeric depletant with N = 40 and κ = 25.1, while the dashed lines correspond to more flexible chains with N = 96 and κ = 3.5. In the lower panel, the solid phase boundary lines correspond to the system (N = 48, κ = 11.5) and the dashed lines to the system (N = 96, κ = 3.5). In both panels, the circles mark the location of the liquid-vapor critical points, while the triangles denote the vapor-liquidsolid triple point coexistence. In the reservoir representation, the triangles marking these three coexisting phases all correspond to the same value of the monomer reservoir density ρ r and therefore lie on a horizontal line (all tie-lines connecting coexisting phases are horizontal in the reservoir representation). As one would expect from the PMF results presented in Figure 1, the phase boundaries move to lower values of ρ r with increasing chain stiffness (solid lines lie below dashed lines in both panels of Figure 2). Likewise, for a given polymeric depletant (N = 96, κ = 3.5), the phase boundaries move to lower values of ρ r with increasing colloid radius, as can be seen by comparing the dashed lines in the upper and lower panels of  Figure 3 replots the data shown in Figure 2 in a system representation [26] in the variables η c -ρ. The transformation from the reservoir to the system representation is achieved [3] via the equation ρ = ρ r (1 − η c ). It is important to note that the range of ρ values in Figure 3, as well as the range of ρ r values in Figure 2, are both sufficiently low, so that no isotropic-nematic transition of semiflexible polymers needs to be taken into account [27,28]. It follows from the results shown in Figure 2 that the monomer reservoir density ρ rc (which corresponds to the liquid-vapor critical point), as well as ρ rt (the monomer reservoir density corresponding to the vapor-liquid-solid triple point), both decrease with increasing chain stiffness κ. This behavior is further illustrated in Figure 4, which plots both ρ rc and ρ rt as functions of the inverse stiffness parameter κ −1 . The upper panel presents the results for the colloid radius R c = 5, while the lower panel shows the results for the colloid radius R c = 20. For both values of the colloid size, ρ rc and ρ rt increase nearly linearly with κ −1 . Furthermore, the slope of the triple point line ρ rt is significantly higher than the slope of the critical point line ρ rc . This behavior is indeed consistent with the phase diagrams shown in Figure 2, where one sees that an increase in chain stiffness (going from dashed to solid phase boundaries) leads to a substantially larger drop in ρ rt (marked by triangles) compared to the drop in ρ rc (marked by circles). By comparing the upper and lower panels of Figure 2, one sees that both ρ rc and ρ rt decrease with increasing colloid radius, which is consistent with W(R) becoming more attractive for larger colloids, as shown in Figure 1. The dependence of the critical and triple reservoir monomer densities on the dimensionless ratio R c /R g is illustrated in Figure 5 for two particular semiflexible chains: (N = 64, κ = 6.2) in the upper panel and (N = 96, κ = 3.5) in the lower panel. One sees that for both depletants ρ rc and ρ rt decrease monotonically with increasing R c /R g , as one would expect based on the results shown in Figures 1 and 2.

Conclusions
In this work, we have developed a Density Functional Theory for the depletion potential of mean force between spherical colloids induced by semiflexible polymer chains. The theory was used to study the effects of the colloid radius (relative to the polymer radius of gyration) and the chain stiffness (for a given value of R g ) on the strength of the depletion attraction. In agreement with earlier self-consistent field theory calculations [21], it was found that depletion attraction becomes stronger for larger colloids and stiffer chains.
The colloid-colloid potential of mean force was calculated for a range of monomer densities, and the results were subsequently used to construct phase diagrams for an effective one-component colloidal system, both in the reservoir and in the system representations of the monomer density. The phase boundaries for vapor-liquid, vapor-solid, and liquid-solid phase coexistence were obtained, as well as vapor-liquid critical points and vapor-liquid-solid triple points. The reservoir monomer densities corresponding to critical and triple point were found to increase nearly linearly with inverse chain stiffness, with the slope of the former being substantially smaller than the latter.
This work can be extended in several important directions. First, it would be of interest to study the morphology of the crystalline phase. To achieve this, the microscopic model must be made more chemically specific in order to enable comparison with existing experimental data obtained by tunneling electron microscopy. Second, given the importance of photocatalytic activity for practical applications of these systems, one could combine the present results with quantum DFT calculations with a goal of estimating the photocatalytic activity and comparing it with previous studies. These directions will be the subjects of future research.