Repulsion between Colloidal Particles Mediated by Nonadsorbing

Jul 23, 2019 - On the other hand, Louis and co-workers performed Monte Carlo (MC) ... molecules, respectively, at lattice site r, and δ denotes the K...
0 downloads 0 Views 2MB Size
Download PDF
Article Cite This: Macromolecules XXXX, XXX, XXX−XXX

pubs.acs.org/Macromolecules

Repulsion between Colloidal Particles Mediated by Nonadsorbing Polymers: Lattice Monte Carlo Simulations and the Corresponding Self-Consistent Field Calculations Pengfei Zhang†,‡ and Qiang Wang*,‡ †

Downloaded via KEAN UNIV on July 24, 2019 at 09:07:53 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.

Center for Advanced Low-Dimension Materials, State Key Laboratory for Modication of Chemical Fibers and Polymer Materials, Donghua University, Shanghai 201620, China ‡ Department of Chemical and Biological Engineering, Colorado State University, Fort Collins, Colorado 80523-1370, United States ABSTRACT: Using a lattice self-consistent field (SCF) theory and the corresponding lattice Monte Carlo (MC) simulations combined with our recently proposed Z method [Zhang, P.; Wang, Q. Soft Matter 2015, 11, 862], we examined athermal homopolymer solutions confined between two parallel and nonabsorbing surfaces and in equilibrium with a bulk solution and accurately calculated the effective interaction between the two surfaces. By directly comparing our MC results with SCF predictions based on the same model system, we were able to quantitatively and unambiguously distinguish the meanfield and the fluctuation contributions to the effective interaction. We found for the first time the fluctuation-induced repulsion between the two confining surfaces at intermediate separation predicted by Semenov and Obukhov [Obukhov, S. P.; Semenov, A. N. Phys. Rev. Lett. 2005, 95, 038305; Semenov, A. N.; Obukhov, S. P. J. Phys.: Condens. Matter 2005, 17, S1747], which is about one order of magnitude stronger than that due solely to the finite chain length as predicted by the SCF theory.

I. INTRODUCTION Adding polymers into colloidal suspensions provides a simple yet pragmatic means of controlling the effective interaction between colloidal particles, thus the phase behavior of various colloidal systems such as emulsions, paints, and milk.1−5 In the “colloid” limit, where the diameter of colloidal particles is much larger than the radius of gyration of polymer chains, the effective interaction between two spherical particles can be well-approximated through that between two parallel planar surfaces via the Derjaguin approximation.6,7 For the simplest system of an athermal homopolymer solution confined between two flat, hard (impenetrable), homogeneous, and parallel surfaces and in equilibrium with a bulk solution, it is widely accepted that at small intersurface separation D there is an effective attraction (usually termed the “depletion attraction”) between the two surfaces.2−4,8−11 This attraction results from the osmotic pressure imbalance between the confined solution, where chains are strongly depleted, and the bulk solution, where chains adopt coil conformations. On the other hand, when D is sufficiently large such that the two depletion layers near each surface do not overlap, the effective interaction disappears. While quite a few experimental results suggest the possible existence of a repulsive region at intermediate D,12−18 there are theoretically many controversies on the physical mechanisms of this repulsion, and a systematic molecular simulation study of this repulsion is also lacking. In this work, we therefore consider such a confined athermal homopolymer solution and focus on the effective interaction, particularly the repulsion, between the two surfaces induced by © XXXX American Chemical Society

nonadsorbing polymers, which are depleted from the surfaces due to their loss of conformational entropy. The mechanism of repulsion at intermediate D was first discussed by Feigin and Napper using the concept of “depletion stabilization” about 40 years ago.19 They stated that when the surfaces approach each other, chains are gradually depleted from the confined region; this is energetically unfavorable in a good solvent, thus leading to the repulsion.19 To support their qualitative argument, they further constructed a phenomenological theory where ideal-chain conformations are used to obtain the polymer concentration and the Flory−Huggins theory to calculate the intersurface potential W (i.e., the work needed to bring two surfaces from an infinite separation to D).20 They found that W exhibits a positive maximum Wmax at D comparable to the chain size R (thus a repulsion between the two surfaces at D somewhat larger than the location of Wmax) and a negative minimum (the depth of which is denoted by Wmin) at D ≈ 0, and that Wmin is significantly larger than Wmax. Because of the large Wmin at strong confinement, the repulsion can only stabilize the colloidal particles kinetically instead of thermodynamically.20 The above work19,20 unfortunately neglects the interaction between chains and their conformational entropy in the confined, inhomogeneous solution; its reliability is therefore questionable. Another study using the ground-state dominance Received: March 18, 2019 Revised: June 24, 2019

A

DOI: 10.1021/acs.macromol.9b00545 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules

two hard surfaces in an expanded isotension ensemble and found that at high enough polymer volume fractions a shortrange repulsion between the two surfaces appears at a separation on the order of the chain bond length. From the above, we see that there exist two possible physical origins for the repulsion between two nonabsorbing surfaces: chain-end effect and fluctuation effect. While the former is reasonably well understood by the GSDA, GSDE, and full SCF theories, the fluctuation-induced repulsion is only studied theoretically,29,30 and its substantiation by molecular simulations or experiments is still lacking; without an accurate determination of the fluctuation-induced repulsion, a quantitative comparison between the magnitudes of these two effects cannot be performed, and the physical mechanism of the repulsion cannot be unambiguously revealed. In this work, using a lattice SCF theory37,38 and the corresponding lattice MC simulations combined with our recently proposed Z method,39 which is much more efficient and accurate than the repulsive-wall method34 used by Louis and co-workers,31−33 we accurately calculate the effective interaction between two neutral surfaces in the presence of nonadsorbing athermal homopolymer solution. We focus on systems in 2D since they exhibit stronger fluctuations than their 3D counterpart;40 this was also suggested by Semenov and Obukhov.29,30 Extending our study to systems in 3D is straightforward and deferred to a future publication. By directly comparing MC simulation results with SCF predictions based on the same model system, we can quantitatively and unambiguously distinguish the mean-field and the fluctuation contributions.

approximation (GSDA) theory, which takes these effects into account, predicts only attraction between the two surfaces.21 Because the GSDA theory is a mean-field theory valid only in the limit of chain length N → ∞, the repulsion (if any) must result from the finite N and/or the fluctuation effects. Indeed, Scheutjens and Fleer performed lattice selfconsistent field (SCF) calculations for finite N and found similar results to those of Feigin and Napper, except that Wmax is so small that the repulsion has practically no effect on the particle interaction.22,23 A more detailed study by van der Gucht and co-workers showed that W obtained from the lattice SCF calculations22,23 is actually a damped oscillatory function of D for semidilute to concentrated solutions,24 with both the period and the decay length of the oscillation proportional to the coil size and independent of the bulk polymer concentration ϕb. The decay of the oscillation, however, is generally so rapid that only the first maximum of W is observable. Moreover, they found no oscillations in dilute solutions because the depletion correlation length there is larger than the decay length of the oscillations, in analogy with the Fisher−Widom criterion25 in simple fluids.24 In parallel to these numerical lattice SCF calculations, Semenov presented an analytical correction to the GSDA theory due to finite N (termed the GSDE theory hereafter).26 His study suggests a repulsive potential at intermediate D, the magnitude of which is inversely proportional to D and is the largest in the θ-solvent around the polymer overlapping concentration.26 Later, he further developed a universal, asymptotically exact GSDE theory, which elucidates that the depletion stabilization originates from the enrichment of chain ends near the surfaces.27 Physically, this finite-N repulsion can be understood in terms of the formation of two “virtual” opposite polymer brushes due to the enrichment of chain ends near the surfaces. Systematic comparisons with their numerical results of a continuum SCF theory28 showed that the GSDE predictions are good when ξ/Rg ≲ 0.1, where ξ and Rg denote the bulk correlation length and ideal-chain radius of gyration, respectively.28 This repulsion alone, however, is usually not strong enough to provide stabilization for colloidal particles with a diameter 7 so that the finite-size effects in our E

DOI: 10.1021/acs.macromol.9b00545 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules

Figure 2. (a) Reduced grand potential per chain ψ* as a function of the reduced exchange chemical potential μ* for confined solutions at various intersurface separations D, where “Bulk” denotes the corresponding result of the bulk solution. (b) Reduced single-chain partition function of the confined solution in the dilute limit Q*1 as a function of D/Re, where the root-mean-square chain end-to-end distance Re is taken to be R e,0 = N − 1 for the SCF result and that of a single self-avoiding chain in bulk (i.e., 8.1710) obtained via enumeration for the MC result, and the dotted and solid horizontal lines mark, respectively, 0 and −μ1*,ex given in section III.1. N = 20.

monotonically with increasing ϕb in both SCF theory (i.e., eq 32) and MC simulations (using the Z method; see ref 39 for details). For small ϕb, Π*,ex ≈ Bϕb2 is found in both cases with B being the second virial coefficient. While SCF theory predicts B = 1/2 and is independent of N, MC simulations give a smaller value of B ≈ 0.153 ± 0.003 for N = 2039 and B ≈ 0.0786 ± 0.003 for N = 80; that B decreases with increasing N is consistent with earlier perturbation calculations.44 Finally, we list in Table 1 the root-mean-square end-to-end distance Re,b for chains of length N = 20 and 80 in the bulk solutions at various ϕb obtained from our MC simulations (performed in a canonical ensemble), which provides the physical length scale for the confined solutions examined below. The approximate treatment of the excluded-volume effects by the SCF theory can be clearly seen by comparing its corresponding prediction R e,0 = N − 1 (i.e., the ideal-chain result regardless of ϕb) with Re,b and is the key for understanding our results of confined solutions as shown below. In general, the intrachain excluded-volume interaction tends to swell a chain, and the interchain excluded-volume interaction tends to compress it. The mean-field approximation inherent in the SCF theory, however, corresponds to the limit of ρ0 → ∞, where the intra- and interchain excluded-volume contributions exactly cancel each other, resulting in the idealchain conformations in the bulk. At finite ρ0 (e.g., 1 used in MC simulations), such cancellation (i.e., the Flory’s ideality hypothesis45) is only approximately true for polymer melts and concentrated solutions (where ϕb is close to 1); for less concentrated solutions (e.g., those in Table 1), the SCF theory predicts much less expanded chain dimensions than the corresponding MC results since it does not account for the fact that a polymer segment is more likely to interact with another segment on the same chain than on a different chain. III.2. Lattice Self-Consistent Field Results of Confined Solution. III.2.1. Grand Potential per Chain and SingleChain Partition Function. The top three curves in Figure 2a show our SCF results of the reduced grand potential per chain ψ* as a function of μ* for the confined solutions of N = 20 at various intersurface separations D. We see that −ψ* decreases monotonically with decreasing μ* at each D and exhibits the exponential behavior of −ψ* ≈ C1 exp(μ*) at small μ* with C1

> 0 decreasing with decreasing D. Similar results are also found for N = 80 (data not shown). To explain this exponential behavior, we note that the interchain interactions are negligible at small μ* (thus small polymer volume fraction), where the system can be approximated by an ideal gas of chains (note that the confinement effects of the impenetrable surfaces still need to be considered). Starting from eq 33 and assuming ρ̂(r) ≪ ρ0, we have | l o o o o ̂ [ρ0 − ρ (̂ r)]! = expm [ − ]! ρ ρ ln ( ) r } ∑ 0 o o o o r n r ~ | l o o o o ̂ ̂ ̂ ≈ expm [ − ] [ − ] − [ − ] ρ ρ ρ ρ ρ ρ ( ) ln ( ) ( ) r r r } ∑ ∑ 0 0 0 o o o o r n r ~ | l o o o o ≈ expm [ρ0 − ρ (̂ r)][ln ρ0 − ρ (̂ r)/ρ0 ] − (ρ0 V − nN )} ∑ o o o o n r ~

∏

≈ exp(ρ0 V ln ρ0 − nN ln ρ0 − ρ0 V ) ≈ (ρ0 ! )V /ρ0nN

where the Stirling approximation is used to obtain the second and the last lines above, ln(1 − x) ≈ − x is used to obtain the third line above, and the O(ρ̂2(r))-term is neglected to obtain V the fourth line above; and thus Z(n, D) ≈ Zn1(D)ρnN 0 /n! (ρ0!) , where Z1(D) ≡ GQ1* is the canonical partition function of a single chain in the confined solution of D lattice layers in the xdirection, and Q*1 (D) ≡ limμ*→−∞Q* is the corresponding reduced single-chain partition function. Inserting this approximate result into eq 3, where the summation over nS is omitted as the incompressibility constraint is no longer needed, gives Ξ ≈ exp[ρ0VμS + Z1(D) exp (μ)ρN0 ]/(ρ0!)V and finally ψ * ≈ −N

Z1(D) ADzLN − 1

exp(μ*) = −NQ 1*(D) exp(μ*)

(35)

after some algebra. Note that Q*1 monotonically decreases from its bulk value of 1 (i.e., in the limit of D → ∞) with decreasing D as shown by the open circles in Figure 2b, due to the restriction of possible chain configurations by the impenetrable surfaces. This therefore explains the behavior that −ψ* F

DOI: 10.1021/acs.macromol.9b00545 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules

Figure 3. Polymer volume fraction profile in the direction perpendicular to the confining surfaces ϕ(x) in the confined solutions at various D obtained from the SCF calculations, where ϕb = 0.1 in (a) and 0.5 in (b) is marked by the dotted horizontal line. Part (c) shows the average polymer volume fraction ϕ̅ of the confined solutions in equilibrium with bulk solutions at various polymer volume fractions ϕb as a function of D/Re, where the root-mean-square chain end-to-end distance Re is taken to be R e,0 = N − 1 for the SCF results and Re,b given in Table 1 for the MC results. N = 20.

monotonically decreases with decreasing D at given μ* shown in Figure 2a. III.2.2. Polymer Volume Fraction Profile and the “Sandwich” Approximation. Figure 3a shows the polymer volume fraction profile ϕ(x) in the confined solutions of N = 20 at various D and ϕb = 0.1. For small D ≲ 5, ϕ(x) at all x is much smaller than ϕb due to the strong confinement of the two impenetrable surfaces and the excluded-volume interactions between polymer segments. With increasing D, ϕ(x) in the interior of the confined solution increases monotonically toward ϕb, while that in the regions close to a surface increases toward some value much smaller than ϕb; in other words, at this small ϕb, the confined solution can be approximately split into the bulk and surface regions for large D ≳ 10. On the other hand, Figure 3b shows that this “sandwich” approximation is applicable to all D at large ϕb = 0.5. Similar results are found for N = 80 (data not shown). The curves with open symbols in Figure 3c show the SCF results of the average polymer volume fraction in the confined solution ϕ̅ ≡ nN/ρ0V = ∑Dx=1ϕ(x)/D as a function of D for various ϕb. We see that ϕ̅ /ϕb < 1 monotonically decreases with decreasing D at given μ*; moreover, the dependence of ϕ̅ on D is weaker at higher μ* (i.e., higher ϕb); similar results are also found for N = 80 (data not shown). These are well understood based on the depletion of chains from the impenetrable surfaces. III.2.3. Surface Tension and Intersurface Potential. As explained in section II.1 and demonstrated in section III.2.2, ϕ(x) at large D can be well described by the “sandwich”

approximation and the validity of this approximation shifts to smaller D with increasing ϕb. In such cases, the two surface layers are decoupled and the confined solution can be well approximated by two single surface layers in contact with the bulk solution. The bottom two curves in Figure 4a therefore show the reduced (mean-field) surface tension γ* as a function of ϕb for both N = 20 and 80. For small ϕb, γ* ∝ ϕb due to the fact that chains are nearly isolated in dilute solutions, which is used above to explain the behavior of ψ* at small ϕb shown in Figure 2. For large ϕb, γ* also exhibits an approximately linear dependence on ϕb but with a larger slope. Moreover, while γ* is slightly smaller for longer N at small ϕb, the opposite occurs at large ϕb. While the intersurface potential vanishes in the above cases, the two surfaces are no longer decoupled when their surface layers overlap at small D. Figures 4b−d thus show the reduced intersurface potential W* as a function of D/Re,0 at ϕb = 0.1, 0.3, and 0.5, respectively, obtained from our SCF calculations. At small ϕb = 0.1 and 0.3, W* < 0 for small D/Re,0 ≲ 1 and increases with increasing D/Re,0, which is just the well-known “depletion attraction”;8 at ϕb = 0.5, this short-range depletion attraction is not accessible due to the use of a lattice. W* then becomes >0 (thus a repulsion) over some intermediate range * at Dm/Re,0 ≈ 1; in the case DA < D < DR with a maximum Wmax of ϕb = 0.5 and N = 20, W*max and Dm are taken as the leftmost data point again due to the use of a lattice. As D further increases, W* exhibits small damped oscillations around 0.24 Table 2 lists Wmax * , Dm/Re,0, DA/Re,0, and DR/Re,0 for various N and ϕb obtained from our lattice SCF calculations; in G

DOI: 10.1021/acs.macromol.9b00545 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules

Figure 4. (a) Reduced surface tension γ*. (b−d) Reduced intersurface potential W*(D) for the confined solution at various bulk polymer volume fractions ϕb obtained from the SCF calculations.

Table 2. Maximum of the Reduced Intersurface Potential W* (Denoted by Wmax * ) and Its Location (Denoted by Dm/Re), and the Linearly Interpolated First and Second Roots of W* (Denoted by DA/Re and DR/Re, Respectively) Bracketing Dm/Re for the Chain Length of (a) N = 20 and (b) N = 80 at Various Bulk Polymer Volume Fractions ϕb Obtained from Our Lattice SCF Calculations and MC Simulations, Where the Root-Mean-Square Chain End-to-End Distance Re Is Taken to Be R e,0 = N − 1 for the SCF Results and Re,b Given in Table 1 for the MC Results (Number in Parentheses Denotes the Estimated Error on the Last Digit) (a) N = 20 ϕb = 0.3

ϕb = 0.1 SCF * Wmax Dm/Re DA/Re DR/Re

MC −5

8.73 × 10 1.376 1.121 2.255

SCF

7.4(1) × 10 1.111 0.943 1.605

−4

−4

5.70 × 10 0.688 0.661 1.793

* Wmax Dm/Re DA/Re DR/Re

6.92 × 10 0.900 0.757 1.885

−3

6.5(3) × 10 1.00 × 10 0.651 0.688 0.549 0.485 1.301 1.723 (b) N = 80

SCF −4

7.5(5) × 10 1.021 0.805 1.332

2.22 × 10 0.450 0.421 1.675

SCF −2

1.213(5) × 10 0.538 0.427 1.076

1.88 × 10 0.459 − 1.639

ϕb = 0.4 MC

−4

ϕb = 0.5 MC

ϕb = 0.3 MC

−5

SCF −3

ϕb = 0.1 SCF

ϕb = 0.4 MC

SCF −3

7.9(2) × 10 0.412 0.349 0.773

particular, W*max and Dm/Re,0 are given by the actual (i.e., not interpolated) data points, DA/Re,0 and DR/Re,0 are given by the linearly interpolated first and second roots of W*, respectively, and DA < Dm < DR. We see that W*max increases with increasing ϕb and decreasing N, and that Dm/Re,0, DA/Re,0 and DR/Re,0 all decrease with increasing ϕb and N (due to the use of a lattice, our reported Dm/Re,0 is the same for the cases of N = 20 at ϕb = 0.3 and 0.4, as well as for the cases of N = 80 at ϕb = 0.4 and 0.5, respectively). We therefore conclude that this intermedi-

3.07 × 10 0.338 0.327 1.648

1.99(2) × 10−2 0.419 0.352 0.978 ϕb = 0.5

MC −4

MC −3

SCF −3

8.5(1) × 10 0.280 0.262 0.616

MC −4

4.43 × 10 0.338 − 1.624

− − − −

ate-range repulsion at the mean-field level is solely due to the finite N. Our mean-field results are consistent with the previous SCF and GSDE predictions.22−24,26−28 III.2.4. Intersurface Pressure and Accuracy of FiniteDifference Formulas. Figure 5 shows the reduced intersurface pressure P* in the confined solution of N = 80 at various ϕb, where the derivative in calculating the normal pressure Pn is approximated by various finite-difference (FD) formulas given in the Appendix. At small ϕb (=0.1), all the FD formulas give H

DOI: 10.1021/acs.macromol.9b00545 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules

Figure 5. Reduced intersurface pressure P* in the confined solution of N = 80 obtained from the SCF calculations using finite-difference formulas of various orders shown in the legend, where “C”, “F”, “CF1”, and “CF2” denote the centered, forward, and two centered-forward finite-difference formulas, respectively, given in the Appendix.

results close to each other. At larger ϕb (=0.3 and 0.5), however, results of the second-order FD formulas are significantly different from those of the fourth- and sixthorder formulas; in particular, P* calculated using the fourthand sixth-order formulas are much closer to 0. This is clearly due to the use of a lattice, where D can only be integers, and the higher-order FD formulas are expected to give more accurate results. On the other hand, using shorter chains of N = 20 enlarges this lattice discretization effect and various FD formulas give even more different results (data not shown). Finally, we note that in the limit of D/Re,0 → 0, −P* should approach the reduced osmotic pressure Π* ≈ 0.006611, 0.06042 and 0.1994, respectively, for bulk solutions of ϕb = 0.1, 0.3 and 0.5 given by the SCF theory; while this is approximately the case for ϕb = 0.1, for ϕb = 0.3 and 0.5 the polymer volume fraction in the confined space is higher than 0 even at the strongest confinement (i.e., D = 2, the leftmost data points in Figure 5) allowed in the lattice SCF calculations. III.3. Lattice Monte Carlo Results of Confined Solution. III.3.1. Grand Potential per Chain and SingleChain Partition Function. The bottom three curves in Figure 2a show the reduced grand potential per chain ψ* as a function of μ* at various D obtained from our MC simulations. We see that the MC results exhibit qualitatively the same behavior as the SCF results, and the same argument as presented in section III.2.1 explains the exponential behavior of − ψMC * at small μ*. At the same μ* and D (where D/Re in the SCF calculations is larger than that in the MC simulations), however, −ψ*MC is much smaller than − ψ*SCF due to the approximate treatment of the excluded-volume effects by the SCF theory. Similarly, the filled circles in Figure 2b show that the reduced single-chain partition function obtained from our MC

simulations for a single chain in the confined system of D lattice layers in the x-direction with A = N, Q*1 (D) = Z1(D)/ ADzN−1 L , exhibits qualitatively the same behavior as that from the SCF calculations. Quantitatively, Q1,MC * < Q1,SCF * at the same D/Re because the (intrachain) excluded-volume interaction, which is approximately treated by the SCF theory, further limits the allowed chain configurations in addition to the restriction due to the confining surfaces; in the limit of D → ∞, MC simulation gives −ln Q1* = μ1*,ex ≡ μex 1 + (N − 1) × ln ρ0 (= −ln (83,779,155/418) for N = 20) instead of 0 as given by the SCF theory. Finally, we note that with μ*,id ≡ μid + ,ex ln(zN−1 L /ρ0) = ln(ϕ̅ /N), the general result of μ* = −ln Q* in the SCF theory can be easily derived from eq 26. III.3.2. Polymer Volume Fraction Profile. The polymer volume fraction profiles ϕ(x) obtained from the MC simulations in general exhibit qualitatively the same behavior as those from the SCF calculations. Quantitatively, the thickness of the surface layer obtained from the MC simulations is larger than that from the SCF calculations. Taking the system of N = 80 and ϕb = 0.3 as an example, in Figure 6 we compare the polymer volume fraction profiles ϕ(x) in the confined solution at D = 8 obtained from MC simulations and SCF calculations. We see that the MC result lies well below the SCF result and that the surface layer from MC simulations is much thicker than that from SCF calculations. While the “sandwich” approximation works well for the SCF result, it fails for the MC result (for this set of parameter values); this is because the SCF theory predicts much less expanded chain dimensions due to its approximate treatment of the excluded-volume effects. The curves with filled symbols in Figure 3c show the grandcanonical-ensemble average of the polymer volume fraction in I

DOI: 10.1021/acs.macromol.9b00545 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules

layers, consistent with the behavior of ϕ(x) shown in Figure 6. Moreover, γMC * − γSCF * increases with increasing ϕb. Figure 7 compares the reduced intersurface potential W* obtained from our MC simulations and SCF calculations at various ϕb. We see that the MC results exhibit qualitatively the same behaviors as the SCF results. Note that for N = 20 at ϕb = 0.5, while the SCF theory predicts a purely repulsive interaction (within our data range), MC simulations give both a short-range depletion attraction and an intermediate-range repulsion. Table 2 also compares the maximum of the reduced intersurface potential W*max and its location Dm/Re, as well as the range of the repulsion DA/Re and DR/Re, obtained from our MC simulations with those from the SCF calculations quantitatively. We see that Wmax * obtained from the MC simulations is larger than that from the SCF calculations by about one order of magnitude, indicating that the intermediate-range repulsion in the MC simulations is mainly due to the fluctuations neglected by the SCF theory. This is consistent with Semenov and Obukhov’s theoretical prediction.29,30 We also see that the former increases with increasing ϕb and has only a weak dependence on N. On the other hand, Dm/Re, DA/Re and DR/Re from the MC simulations all decrease with increasing ϕb and N, and they are smaller than their counterparts from the SCF calculations (except Dm/Re and DA/Re in the case of ϕb = 0.1 and N = 80). III.3.4. Intersurface Pressure. Figure 8 shows the reduced intersurface pressure P* in the confined solution of N = 80 at various ϕb obtained from our MC simulations. As expected, P* approaches the negative reduced osmotic pressure −Π* at small D/Re,b (due to the complete depletion of polymer chains in the confined space) and 0 at large D/Re,b. Compared to the SCF results shown in Figure 5, we see that the various FD formulas give more consistent results here. Clearly, in all cases P* from the MC simulations exhibits a small positive maximum, with its value increasing and its location decreasing with increasing ϕb. Given the good agreement among the results of various FD formulas, we conclude that this maximum is not due to the numerical errors of the FD formulas but the fluctuations neglected by the SCF theory. In other words, we find the fluctuation-induced repulsion between the two confining surfaces at intermediate separation predicted by Semenov and Obukhov,29,30 with its strength much larger than

Figure 6. Polymer volume fraction profiles in the direction perpendicular to the confining surfaces ϕ(x) in the confined solution at D = 8, where the bulk concentration ϕb = 0.3 is marked by the dotted horizontal line. N = 80.

the confined solution ϕ̅ obtained from our MC simulations as a function of D/Re at various ϕb (i.e., μ*) for N = 20. We again see that the MC results exhibit qualitatively the same behavior as the SCF predictions, but with ϕ̅ MC < ϕ̅ SCF at the same D/Re for ϕb = 0.1 and 0.3 again due to the approximate treatment of the excluded-volume effects by the SCF theory; with increasing ϕb, the intrachain excluded-volume interaction is increasingly screened by the interchain excluded-volume interaction, which leads to ϕ̅ MC ≈ ϕ̅ SCF at ϕb = 0.5. Finally, we note that the partition coefficient ϕ̅ /ϕb and the equivalent free energy of confinement −ln(ϕ̅ /ϕb) have been examined by Cifra and coworkers46,47 and Wang and co-workers48,49 using lattice MC simulations; the latter quantity and the related confinement force50,51 d ln(ϕ̅ /ϕb)/dD, however, are fundamentally different from the intersurface potential and pressure studied in this work. III.3.3. Surface Tension and Intersurface Potential. The top two curves in Figure 4a show that the reduced surface tension γ* obtained from our MC simulations exhibits qualitatively the same behavior as the corresponding SCF result. On the other hand, γMC * > γSCF * at the same ϕb; this is because the approximate treatment of the excluded-volume effects by the SCF theory reduces the loss of the chain conformational entropy due to the formation of the surface

Figure 7. Reduced intersurface potentials W* obtained from the SCF calculations and MC simulations at various ϕb as a function of D/Re, where the root-mean-square chain end-to-end distance Re is taken to be R e,0 = N − 1 for the SCF results and Re,b given in Table 1 for the MC results. N = 20 in (a) and 80 in (b). J

DOI: 10.1021/acs.macromol.9b00545 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules

Figure 8. Reduced intersurface pressure P* in the confined solution of N = 80 obtained from MC simulations using finite-difference formulas of various orders shown in the legend, where “C”, “F”, “CF1”, and “CF2” denote the centered, forward, and two centered-forward finite-difference formulas, respectively, given in the Appendix. In each part, the dotted horizontal line marks the negative reduced osmotic pressure −Π* of the bulk solution.

fluctuation-induced repulsion between the two confining surfaces at intermediate separation predicted by Semenov and Obukhov,29,30 and that SCF and MC results exhibit qualitatively the same behavior with their quantitative differences summarized below. While the MC results of the bulk solution at N = 20 were reported in our previous work,39 our results here further show that for small ϕb, μ*SCF ≈ ln(ϕb/N), which is due to the meanfield approximation leading to ideal-chain conformations in the bulk system; in contrast, μMC * ≈ ln(ϕb/N) + μ1*,ex, where μ1*,ex > 0 denotes the (reduced) excess chemical potential of a single self-avoiding chain. Their difference is due to the excludedvolume effects approximately treated by the SCF theory, which is also the key for understanding the results of confined solutions. Moreover, while the second virial coefficient BSCF = 1/2 and is independent of N, BMC is much smaller and decreases with increasing N. For confined solutions, the negative reduced grand potential per chain −ψ* decreases monotonically with decreasing μ* at each D and exhibits the exponential behavior of −ψ* ≈ C1exp(μ*) at small μ* with C1 > 0. The latter can be understood by approximating the system as an ideal gas of (confined) chains, and the restriction of possible chain configurations by the impenetrable surfaces further explains the monotonic decrease of C1 (thus −ψ*) with decreasing D. At the same μ* and D, however, −ψMC * is much smaller than −ψSCF * due to the approximate treatment of the excludedvolume effects by the SCF theory.

that of the intermediate-range repulsion caused by the finite N as predicted by the SCF theory. Moreover, by comparing the results at N = 20 (data not shown) and 80 at the same ϕb, we find that the strength of the repulsive barrier depends weakly on N; this is consistent with the results of W* presented in the previous section.

IV. SUMMARY Using a lattice self-consistent field (SCF) theory37,38 and the corresponding lattice Monte Carlo (MC) simulations combined with our recently proposed Z method,39 we have examined athermal homopolymer solutions confined between two parallel and nonabsorbing surfaces and in equilibrium with a bulk solution on the square lattice in 2D. In particular, the Z method enables us to accurately calculate the reduced intersurface potential W*. We have used two chain lengths N = 20 and 80, with the total number of polymer segments and solvent molecules on each lattice site ρ0 = 1 (i.e., the self- and mutual-avoiding walk) in our MC simulations and ρ0 → ∞ (i.e., the mean-field approximation) in the SCF calculations, and varied the surface separation D and the reduced exchange chemical potential μ* (or equivalently the bulk polymer volume fraction ϕb). By directly comparing our MC results with SCF predictions based on the same model system, we can quantitatively and unambiguously distinguish the mean-field and the fluctuation contributions to the effective interaction between the two surfaces. We have found for the first time the K

DOI: 10.1021/acs.macromol.9b00545 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules On the other hand, for large D and ϕb, the polymer volume fraction profile in the direction perpendicular to the surfaces ϕ(x) can be described by the “sandwich” approximation, which splits the confined solution into the bulk and surface regions. In such cases, the two surface regions are decoupled and W* ≈ 0, and we have found that the reduced surface tension γ* ∝ ϕb at small ϕb and exhibits an approximately linear dependence on ϕb at large ϕb but with a larger slope; we have also found that, while γ* is slightly smaller for longer N at small ϕb, the opposite occurs at large ϕb. Since the SCF theory predicts much less expanded chain dimensions again due to its approximate treatment of the excluded-volume effects, however, the surface region given by MC simulations is much thicker than the SCF prediction and we have found that γMC * − γSCF * > 0 increases with increasing ϕb. When the two surface regions are no longer decoupled, we have found the well-known short-range “depletion attraction”,8 where W* < 0 for small D/Re < DA/Re (= 0.3−1.1) and decreases monotonically with decreasing D; here the rootmean-square chain end-to-end distance Re is taken to be N − 1 (i.e., that of ideal chains) for the SCF results and that given in Table 1 for the MC results. We have also found an intermediate-range repulsion, where W* > 0 for DA/Re < D/Re < DR/Re (= 0.6−2.3) and exhibits a maximum Wmax * at Dm ∈ (DA, DR). Wmax * increases with increasing ϕb, and Dm/Re, DA/Re and DR/Re all decrease with increasing ϕb and N. Consistent with the previous SCF and GSDE predictions,22−24,26−28 W*max obtained from our SCF calculations decreases with increasing N, indicating that the intermediate-range repulsion at the mean-field level is solely due to the finite N (i.e., the chain ends). In contrast, W*max obtained from our MC simulations is about one order of magnitude larger than that from the SCF calculations and has only weak dependence on N; we therefore attribute the intermediate-range repulsion found in our MC simulations mainly to the system fluctuations neglected by the SCF theory, consistent with Semenov and Obukhov’s theoretical predition.29,30 Furthermore, Dm/Re, DA/Re and DR/Re obtained from our MC simulations are smaller than their counterparts from the SCF calculations (except Dm/Re and DA/Re in the case of ϕb = 0.1 and N = 80). Finally, while our results of the reduced intersurface pressure P* are consistent with those of W*, the lattice discretization

effect requires higher accuracy of P* than given by the commonly used second-order finite-difference formulas in the SCF calculations. In contrast, the various finite-difference formulas give more consistent MC results, where P* exhibits a small positive maximum in all cases due to the fluctuations neglected by the SCF theory, thus confirming the fluctuationinduced repulsion between the two confining surfaces at intermediate separation predicted by Semenov and Obukhov.29,30

■

APPENDIX. FINITE-DIFFERENCE FORMULAS OF THE SECOND-, FOURTH-, AND SIXTH-ORDER ACCURACY FOR CALCULATING THE FIRST-ORDER DERIVATIVE Because a first-order derivative is involved in the calculation of the intersurface pressure (as shown in eq 5), for a lattice model as adopted here, we have to approximate the derivative by the finite difference (FD). To estimate errors introduced by the FD approximation, we use different FD formulas and compare their results. For completeness, we present below the second-, fourth-, and sixth-order FD formulas for estimating the firstorder derivative f i′ of a function f(x) at x = i, where the function value is denoted by f i. With the lattice spacing taken as the unit of length, the second-order FD formulas are l( −3f + 4f − f )/2 for i = 1 o o i i+1 i+2 o f i′ = m o o for i ≥ 2 o(fi + 1 − fi − 1 )/2 n

which are termed the forward (F) and centered (C) formulas, respectively; the fourth-order FD formulas are l(− 25f + 48f − 36f + 16f − 3f )/12 for i = 1 o o i i+1 i+2 i+3 i+4 o o o o o (− 3fi − 1 − 10fi + 18fi + 1 − 6fi + 2 + fi + 3 )/12 for i = 2 m o o o o o o for i ≥ 3 o(fi − 2 − 8fi − 1 + 8fi + 1 − fi + 2 )/12 n

f i′ =

which are termed the forward (F), centered-forward (CF1), and centered (C) formulas, respectively; and the sixth-order FD formulas are

l −(147fi − 360fi + 1 + 450fi + 2 − 400fi + 3 + 225fi + 4 − 72fi + 5 + 10fi + 6 )/60 o o o o o o o o o o−(10fi − 1 + 77fi − 150fi + 1 + 100fi + 2 − 50fi + 3 + 15fi + 4 − 2fi + 5 )/60 f i′ = m o o o o(2fi − 2 − 24fi − 1 − 35fi + 80fi + 1 − 30fi + 2 + 8fi + 3 − fi + 4 )/60 o o o o o o ( −f + 9fi − 2 − 45fi − 1 + 45fi + 1 − 9fi + 2 + fi + 3 )/60 o n i−3

for i = 2 for i = 3 for i ≥ 4

Notes

which are termed the forward (F), centered-forward-1 (CF1), centered-forward-2 (CF2), and centered (C) formulas, respectively. Note that the sixth-order FD formulas were used in our previous work.52,53

■

for i = 1

The authors declare no competing financial interest.

■

AUTHOR INFORMATION

ACKNOWLEDGMENTS

This work was supported by NSF CAREER Award CBET-

Corresponding Author

0847016, which is gratefully acknowledged. P.Z. also acknowl-

*E-mail: [email protected]. ORCID

edges the financial support provided by the National Natural

Qiang Wang: 0000-0001-6456-9288

Science Foundation of China (21803011). L

DOI: 10.1021/acs.macromol.9b00545 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules

■

(26) Semenov, A. N. Theory of Long-Range Interactions in Polymer Systems. J. Phys. II 1996, 6, 1759. (27) Semenov, A. N. Theory of Colloid Stabilization in Semidilute Polymer Solutions. Macromolecules 2008, 41, 2243−2249. (28) Shvets, A. A.; Semenov, A. N. Effective Interactions between Solid Particles Mediated by Free Polymer in Solution. J. Chem. Phys. 2013, 139, 054905. (29) Obukhov, S. P.; Semenov, A. N. Long-Range Interactions in Polymer Melts: The Anti-Casimir Effect. Phys. Rev. Lett. 2005, 95, 038305. (30) Semenov, A. N.; Obukhov, S. P. Fluctuation-Induced LongRange Interactions in Polymer Systems. J. Phys.: Condens. Matter 2005, 17, S1747. (31) Bolhuis, P. G.; Louis, A. A.; Hansen, J. P.; Meijer, E. J. Accurate Effective Pair Potentials for Polymer Solutions. J. Chem. Phys. 2001, 114, 4296−4311. (32) Louis, A. A.; Bolhuis, P. G.; Meijer, E. J.; Hansen, J. P. Polymer Induced Depletion Potentials in Polymer-Colloid Mixtures. J. Chem. Phys. 2002, 117, 1893−1907. (33) Addison, C. I.; Louis, A. A.; Hansen, J. P. Influence of Solvent Quality on Polymer Solutions: A Monte Carlo Study of Bulk and Interfacial Properties. J. Chem. Phys. 2004, 121, 612−620. (34) Dickman, R. New Simulation Method for the Equation of State of Lattice Chains. J. Chem. Phys. 1987, 87, 2246−2248. (35) Stukan, M.; Ivanov, V.; Muller, M.; Paul, W.; Binder, K. Finite size effects in pressure measurements for Monte Carlo simulations of lattice polymer models. J. Chem. Phys. 2002, 117, 9934−9941. (36) Broukhno, A.; Jönsson, B.; Åkesson, T.; Vorontsov-Velyaminov, P. N. Depletion and Bridging Forces in Polymer Systems: Monte Carlo Simulations at Constant Chemical Potential. J. Chem. Phys. 2000, 113, 5493−5501. (37) Scheutjens, J.; Fleer, G. F. Statistical Theory of the Adsorption of Interacting Chain Molecules. 1. Partition Function, Segment Density Distribution, and Adsorption Isotherms. J. Phys. Chem. 1979, 83, 1619−1635. (38) Zhang, P.; Li, B.; Wang, Q. Quantitative Study of Fluctuation Effects by Fast Lattice Monte Carlo Simulations. III. Homopolymer Brushes in an Explicit Solvent. Macromolecules 2012, 45, 2537−2550. (39) Zhang, P.; Wang, Q. Calculating Pressure in Polymer Lattice Simulations. Soft Matter 2015, 11, 862. (40) Zhang, P.; Yang, D.; Wang, Q. Quantitative Study of Fluctuation Effects by Fast Lattice Monte Carlo Simulations. V. Incompressible Homopolymer Melts. J. Phys. Chem. B 2014, 118, 12059−12067. (41) de Gennes, P.-G.; Brochard-Wyart, F.; Quere, D. Capillarity and Wetting Phenomena Drops, Bubbles, Pearls, Waves; Springer: New York, 2004. (42) Fredrickson, G. The Equilibrium Theory of Inhomogeneous Polymers; Clarendon Press: Oxford, 2006. (43) Press, W. H.; Teukolsky, S. A.; Vetterling, W. T.; Flannery, B. P. Numerical Recipes in FORTRAN: the Art of Scientific Computing; Cambridge University Press: 2002. (44) Yamakawa, H. Modern Theory of Polymer Solutions; Harper & Row: New York, 2007. (45) de Gennes, P.-G. Scaling Concepts in Polymer Physics; Cornell University Press: Ithaca, NY, 1979. (46) Bleha, T.; Cifra, P.; Karasz, F. E. The Effects of Concentration on Partitioning of Flexible Chains into Pores. Polymer 1990, 31, 1321. (47) Cifra, P.; Bleha, T. Partition Coefficients and the Free Energy of Confinement from Simulations of Nonideal Polymer Systems. Macromolecules 2001, 34, 605. (48) Wang, Y.; Teraoka, I. Computer Simulation of Semidilute Polymer Solutions in Confined Geometry: Pore as a Microscopic Probe. Macromolecules 1997, 30, 8473. (49) Jiang, W.; Wang, Y. Thermodynamics and Partitioning of Homopolymers into a Slit - A Grand Canonical Monte Carlo Simulation Study. J. Chem. Phys. 2004, 121, 3905.

REFERENCES

(1) Napper, D. H. Polymer Stabilization of Colloidal Dispersions; Academic: London, 1983. (2) Lekkerkerker, H.; Tuiner, R. Colloids and Depletion Interaction; Springer: 2011. (3) Tuinier, R.; Rieger, J.; de Kruif, C. G. Depletion-Induced Phase Separation in Colloid-Polymer Mixtures. Adv. Colloid Interface Sci. 2003, 103, 1−31. (4) Binder, K.; Virnau, P.; Statt, A. Perspective: The Asakura Oosawa model: A Colloid Prototype for Bulk and Interfacial Phase Behavior. J. Chem. Phys. 2014, 141, 140901. (5) Poon, W. The Physics of a Model Colloid-Polymer Mixture. J. Phys.: Condens. Matter 2002, 14, R859−R880. (6) Derjaguin, B. V. Analysis of Friction and Adhesion, IV. The Theory of the Adhesion of Small Particles. Colloid Polym. Sci. 1934, 69, 155−164. (7) Israelachvili, J. Intermolecular and Surface Forces, 3rd ed.; Academic Press: London, 2011. (8) Asakura, S.; Oosawa, F. On Interaction between Two Bodies Immersed in a Solution of Macromolecules. J. Chem. Phys. 1954, 22, 1255−1256. (9) Asakura, S.; Oosawa, F. Interaction between Particles Suspended in Solutions of Macromolecules. J. Polym. Sci. 1958, 33, 183−192. (10) Vrij, A. Polymers at Interfaces and the Interactions in Colloidal Dispersions. Pure Appl. Chem. 1976, 48, 471. (11) Wang, Z.-G. 50th Anniversary Perspective: Polymer ConformationA Pedagogical Review. Macromolecules 2017, 50, 9073− 9114. (12) Ogden, A. L.; Lewis, J. A. Effect of Nonadsorbed Polymer on the Stability of Weakly Flocculated Suspensions. Langmuir 1996, 12, 3413. (13) Zhang, X.; Servos, M. R.; Liu, J. Ultrahigh Nanoparticle Stability against Salt, pH, and Solvent with Retained Surface Accessibility via Depletion Stabilization. J. Am. Chem. Soc. 2012, 134, 9910. (14) Vincent, B.; Luckham, P. F.; Waite, F. A. The Effect of Free Polymer on the Stability of Sterically Stabilized Dispersions. J. Colloid Interface Sci. 1980, 73, 508. (15) Cowell, C.; Li-In-On, R.; Vincent, B. Reversible Flocculation of Sterically-stabilised Dispersions. J. Chem. Soc., Faraday Trans. 1 1978, 74, 337. (16) Clarke, J.; Vincent, B. Nonaqueous Silica Dispersions Stabilized by Terminally Anchored Polystyrene: The Effect of Added Polymer. J. Colloid Interface Sci. 1981, 82, 208. (17) Kozer, N.; Kuttner, Y. Y.; Haran, G.; Schreiber, G. ProteinProtein Association in Polymer Solutions: From Dilute to Semidilute to Concentrated. Biophys. J. 2007, 92, 2139. (18) Kuhl, T. L.; Berman, A. D.; Hui, S. W.; Israelachvili, J. N. Part 1. Direct Measurement of Depletion Attraction and Thin Film Viscosity between Lipid Bilayers in Aqueous Polyethylene Glycol Solutions. Macromolecules 1998, 31, 8250−8258. (19) Feigin, R. I.; Napper, D. H. Stabilization of Colloids by Free Polymer. J. Colloid Interface Sci. 1980, 74, 567−571. (20) Feigin, R. I.; Napper, D. H. Depletion Stabilization and Depletion Flocculation. J. Colloid Interface Sci. 1980, 75, 525−541. (21) Joanny, J. F.; Leibler, L.; de Gennes, P. G. Effects of Polymer Solutions on Colloid Stability. J. Polym. Sci., Polym. Phys. Ed. 1979, 17, 1073−1084. (22) Scheutjens, J.; Fleer, G. F. Effect of Polymer Adsorption and Depletion on the Interaction Between Two Parallel Surfaces. Adv. Colloid Interface Sci. 1982, 16, 361; ibid. 1983, 18, 309. (23) Fleer, G. J.; Scheutjens, J. M. H. M.; Vincent, B. The Stability of Dispersions of Hard Spherical Particles in the Presence of Nonadsorbing Polymer. ACS Symp. Ser. 1984, 240, 245. (24) van der Gucht, J.; Besseling, N. A. M.; van Male, J.; Cohen Stuart, M. A. Coil Size Oscillatory Packing in Polymer Solutions near a Surface. J. Chem. Phys. 2000, 113, 2886−2893. (25) Fisher, M. E.; Wiodm, B. Decay of Correlations in Linear Systems. J. Chem. Phys. 1969, 50, 3756. M

DOI: 10.1021/acs.macromol.9b00545 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules (50) Bleha, T.; Cifra, P. Free Energy and Confinement Force of Macromolecules in a Slit at Full Equilibrium with a Bulk Solution. Polymer 2003, 44, 3745−3752. (51) Bleha, T.; Cifra, P. Polymer-Induced Depletion Interaction between Weakly Attractive Plates. Langmuir 2004, 20, 764−770. (52) Wang, Q. Systematic and Simulation-free Coarse Graining of Multi-component Polymeric Systems: Structure-based Coarse Graining of Binary Polymer Blends. Polymer 2017, 117, 315−330. (53) Qiu, W.; Li, B.; Wang, Q. Lattice Self-Consistent Field Calculations of Ring Polymer Brushes. Soft Matter 2018, 14, 1887− 1896.

N

DOI: 10.1021/acs.macromol.9b00545 Macromolecules XXXX, XXX, XXX−XXX