Integral Equation Study of Particle Confinement Effects in a Polymer

J. Phys. Chem. C 2007, 111, 15625-15633

15625

Integral Equation Study of Particle Confinement Effects in a Polymer/Particle Mixture† Douglas Henderson,* Andrij Trokhymchuk,‡ and Yurij V. Kalyuzhnyi‡ Department of Chemistry and Biochemistry, Brigham Young UniVersity, ProVo, Utah 84602

Richard H. Gee and Naida Lacevic Chemistry and Material Science Directorate, Lawrence LiVermore National Laboratory, UniVersity of California, LiVermore, California 94550 ReceiVed: May 10, 2007; In Final Form: August 20, 2007

A product-reactant Ornstein-Zernike (PROZA) extension of the Wertheim’s multidensity integral equation theory techniques is applied to evaluate the structuring of the PDMS polymer melt in the bulk and in the presence of a single and a pair of silica surfaces. By comparing with simulation data, we find that the PROZA theory is not very accurate in reproducing the details of the total (the sum of intramolecular and intermolecular correlations) atom-atom radial distribution function of the united atom PDMS model. At the same time, the inhomogeneous calculations based on the PROZA theory are qualitatively correct in describing the local density distribution of the silicon atoms and a single methyl group near silica surfaces. Relating the flat surface geometry to the spherical geometry by means of the Derjaguin approximation, the results obtained for the two-surface system are used to discuss the two-particle aggregation tendency as well as some general aspects of the origin of the reinforcement phenomenon in a polymer/particle composite.

I. Introduction One of the most pronounced features of polymer melts (for instance, various silicones such as polydimethylsiloxane or PDMS) is its ability to dramatically change its mechanical properties with the addition of an inorganic filler to the pure polymer phase. In particular, the strength of the polymer/particle system increases as the filler concentration increases. The mechanism behind this reinforcement is poorly understood, and consequently, the rational control of the resulting mechanical properties is difficult and usually is based on empirical assumptions. Currently, it is believed that the primary reinforcing mechanism in this polymer/particle system has its origin in some peculiar features of particle-polymer interactions.1 At the same time, there are arguments indicating that taking into account particle-polymer interaction alone is not enough to explain the phenomenon of reinforcement. In particular, the interaction between a single-particle surface and polymer molecules does not depend on how many particles are embedded into a polymer melt. However, on the other hand, the properties of the polymer strongly depend on the fraction of filler particles.2,3 The latter observation suggests that, besides the polymer-filler interaction, there are features of the geometrical restrictions, imposed on the polymer by the presence of particles, which have a major impact on the resulting properties of a composite material. Computer simulations are the most efficient tool to study the polymer/particle systems with complex chemical structures,4-7 performing, if necessary, even a full atomistic level of modeling.1 At the same time, computer simulations for the particular polymer/particle system potentially could be expected to experience problems due to the presence of a large amount of filler particles. In fact, the length scale (diameter) of filler particles †

Part of the “Keith E. Gubbins Festschrift”. Permanent address: Institute for Condensed Matter Physics, Svientsitskoho 1, 79011 Lviv, Ukraine. ‡

usually is much larger than the equivalent polymer length scale (diameter of a polymer bead or radius of gyration or transverse diameter of the polymer chains). Therefore, simulations of a polymer/particle system with realistic particle volume fractions will require that an extremely large number of the polymer molecules be explicitly employed. In contrast, integral equation theory8-10 or density functional theory11 studies do allow one to deal with the system comprising realistic numbers (densities) of both the filler particles and the surrounding polymer melt. However, integral equation theory techniques of the systems of complex chemical structures requires the involvement of some approximations, and the effect of these approximations often is not known in advance. An essential task of any theoretical approach in the case of polymer/particle system is the proper treatment of the polymer subsystem. A natural way to deal with a polymer melt in the integral equation approach to the polymer/ particle mixture is to exploit the methods developed for the pure (or bulk) polymer systems. Although there are two most popular methods within the integral equation theory based either on the reference interaction site model (RISM) approach12 or on the Wertheim’s multidensity approach for associating fluids,13,14 both being successfully applied to the bulk polymers, only the RISM method, so far, has been applied to the polymer/particle problem. The purpose of the present article is to apply the integral equation theory based on the Wertheim’s multidensity formalism to deal with a particular PDMS polymer melt, aiming to obtain an insight into the polymer structuring near a single and molecularly smooth model silica surface and in the slit between two such surfaces. The remainder of the paper is organized as follows. In the following two sections, we describe an implementation of the integral equation theory for the polymer/particle system and pure polymer system. In particular, in section III, we present a rather detailed scheme of the so-called productreactant Ornstein-Zernike approach (PROZA) that will be used

10.1021/jp073582g CCC: $37.00 © 2007 American Chemical Society Published on Web 10/03/2007

15626 J. Phys. Chem. C, Vol. 111, No. 43, 2007

Henderson et al.

in this study to deal with a polymer subsystem. Section IV proceeds with results obtained for a pure and confined PDMS polymer melt. The paper concludes in section V with a summary and discussion. II. Integral Equation Approach and its Application to a Polymer/Particle System A. Initial Polymer/Particle System. Typically, the integral equation technique is based on the Ornstein-Zernike (OZ) equation15

hλµ(r12) ) cλµ(r12) +

∑ Fν ∫ hλν(r13)cνµ(r32)dr3

the polymer molecules as well as any density Fs of the filler particles. However, we restrict ourselves to the case of a oneand two-particle limit. To proceed, let us assume in eq 1 that there are only one or two particles in the mixture. Thus, the particles are present at infinite dilution with zero concentration, Fs ) 0. Further, we assume that the product of the particle density and the particle size D, which is large, vanishes, that is, FsD3 ) 0; otherwise, the particle would affect the entire polymer matrix. In other words, the particles are part of the bulk polymer at a vanishingly small density, that is, Fs f 0. Using this fact, eq 2 results in a system of three equations

(1)

hpp ) cpp + Fpcpp X hpp

(2)

where the Greek subscripts range over all of the system components, that is, polymer (p) as well as particles (s), and rR is the position of particle R. Finally, rλµ ) |rλ - rµ|. The function hλµ(r) ) gλµ(r) - 1 is the total correlation function for a pair of particles of components λ and µ that are separated by the distance r. The function gλµ(r) is the pair or radial distribution function (RDF). The function cλµ(r) is the direct correlation function. The important features of the polymer/particle system of interest are (i) the simultaneous presence of two different mass scales associated with characteristic mass mp of the polymer molecule species and characteristic mass ms . mp of the filler particles and (ii) the simultaneous presence of two different length scales associated with characteristic size d of the polymer molecule species and characteristic size D . d of the filler particles. Issue (i) results in different mobilities of the filler particles and polymer molecules. In particular, this suggests that, to a certain extent, filler particles can be considered as immobile obstacles, while the polymer molecules can be treated as being adsorbed into a disordered environment due to the filler particles. Such a description of the polymer/particle systems can be performed using the so-called replica Ornstein-Zernike (ROZ) formalism originally suggested by Given and Stell16 for a hardsphere fluid/hard-sphere particle system and later extended to a variety of fluid/particle systems, including the polymer/particle systems.17-19 Owing to the loading conditions of the experimentally used polymer/particle systems, the consequence of issue (ii) results in a number density Fp ) Np/V of polymer molecule species that is much larger than the number density Fs ) Ns/V of filler particles. Because of this fact, what we will call the ordinary OZ integral equation theory studies usually are performed for a zero particle density (Fs f 0), which allows for detailed analysis of a polymer structuring in the one- and two-particle limits of the initial polymer/particle system. It is important to emphasize that the replica OZ approach suggests that changes in the properties of the polymer/particle composite system occur due to an effect of the quasi-static disordered particle environment on the polymer melt. In contrast, the ordinary OZ integral equation approach assumes that polymer melt affects the phase behavior of the particle subsystem by inducing the effective interaction between filler particles. Although it is very probable that what occurs in reality is a superposition of these two mechanisms, we felt the replica OZ approach to be beyond the scope of the present study. This is mainly because of an explorative character of the PROZA integral equation study and the fact that similar studies within the PRISM approach have been performed using an ordinary OZ integral equation approach. B. Polymer Melt Next to a Surface. Formally, eq 1 allows one to study the polymer/particle system at any density Fp of

hps ) cps + Fpcsp X hpp

(3)

hoss ) coss + Fpcsp X hps

(4)

ν)s,p

and

For notational convenience, in all of the equations above, we have omitted the dependence on r; the symbol X denotes a convolution in r space. The first equation, eq 2, is just the OZ equation for the homogeneous polymer melt. Equations 3 and 4 involve the presence of the geometrical restriction imposed by a giant filler particle that is equivalent to the effect of a single flat surface. Due to this, eq 3 plays an important role in the statistical mechanics of inhomogeneous systems, that is, systems involving diluted mesoscopic objects. Namely, eq 3 describes the polymer inhomogeneity by means of the local density distribution

Fp(r) ) Fp[hps(r) - 1]

(5)

of the polymer species near the giant particle or flat surface. In turn, eq 4 describes the correlations between the surfaces of two giant particles mediated by the polymer melt and, consequently, yields information that can be used to evaluate the effective interaction between a pair of silica particles. Because the particles are large, the region between the pair of particles can be approximated as a slit of thickness H ) (r - D)/d. An important equilibrium property of the polymer melt confined by a slit-like pore is the pressure acting on the surfaces in the direction perpendicular to the surfaces, or local stress or normal pressure, PN. It can be calculated from the density profiles, and as a function of gap width H, it is given by20-22

βPN(H) ) -β

∑a ∫0

H/2

∂UpW a (z,H) Fa(z,H)dz ∂z

(6)

where z is the distance from the center of the polymer unit to the surface. The normal pressure, PN(H), measured relative to the bulk polymer pressure, PB, that is, the pressure of a pure polymer, will define the solvation or disjoining pressure, Π(H)

Π(H) ) PN(H) - PB

(7)

The effect of surface curvature, if necessary, can be restored to some extent by using the Derjaguin approximation.23 According to Derjaguin, the disjoining pressure between two surfaces can be measured by displacing one of the surfaces a distance of ∆H and the work per unit area to bring the surfaces from infinity to the separation H, or the thermodynamic energy

Integral Equation Study of Particle Confinement

J. Phys. Chem. C, Vol. 111, No. 43, 2007 15627

change of the system (excess film free energy), F, can be calculated as

F(H) ) 2π R

∫H∞ Π(H′)dH′

(8)

where R ) D/2 is the radius of the filler particle and is related to the curvature of the surface. The quantity F(H)/R can be determined directly in a “surface force apparatus” and is called the solvation force. C. Pure Polymer Systems. The integral equation technique that is to be applied to deal with polymer molecules differs from the ordinary IET approach that deals with simple molecular liquids. Currently, several integral equation theories that enable one to predict the structure and thermodynamic properties of the fluid of chain molecules in a bulk phase are available. Among those, the most promising are the theories that are based either on the reference interaction site model (RISM) approach12 or on the multidensity approach for associating fluids.13,14 The theories of the first type are known under the general name polymer RISM (PRISM). Originally, PRISM was developed by Curro and Shweizer24,25 as a simple application of the RISM approach of Chandler and Andersen12 to a fluid of flexible ring molecules. Since then, several modifications have been proposed and used in a large number of applications. Two review papers26,27 give reasonable insight into the older developments of PRISM. Later developments have been reviewed recently.28 Originally, the PRISM approach appeared not to be a selfcontained theory, that is, certain input information on the intramolecular distribution functions was needed. Several different approximations for the intramolecular correlations have been proposed. Later, a self-consistent version of the theory was developed with the intramolecular distribution function obtained from a one-chain computer simulation, that is, the simulation of one molecule with an effective potential between the monomers of the chain that is obtained from the solution of the PRISM equation during the previous iteration cycle. Other self-contained versions of the theory are also available; however, they appear to be less successful. In general, PRISM provides reasonable qualitative, and sometimes quantitative, predictions for the structure of dense polymer systems (melts and alloys), polymer-colloidal mixtures. It is not utilized extensively in the prediction of thermodynamic properties. The version of the multidensity theory that is represented by thermodynamic perturbation theory (TPT)14,29 and its modifications, that is, statistical associating fluid theory (SAFT), proves to be very successful in predicting thermodynamic properties of a fluid of chain molecules and its mixtures. An integral equation version of Wertheim’s multidensity theory for associating fluids13,14 has been also developed and applied to study the properties of the chain molecule fluids.30,31 This extension of the Wertheim’s theory is known as a product-reactant Ornstein-Zernike approach (PROZA). The PROZA is a general statistical mechanical theory of reacting mixtures in which the product molecules and reacting molecules are treated on the same footing. When applied to polymerizing monomers, PROZA yields the mean monomer-monomer pair correlation function, from which the thermodynamics of a fluid of polymerizing molecules can be obtained directly. In the completeassociation limit, the PROZA becomes a theory for the structure and thermodynamics of a fluid of a fully polymerized chain. Several different approximations (PY-like, MSA-like, and HNC-

like) have been proposed and applied to different macromolecular systems, which include hard-sphere chain fluids,32-34 hard-sphere ring fluids,35 and star molecule fluids,36,37 chain fluids with Yukawa38,39 and Coulomb40 interactions between monomers, an so forth. Unlike the PRISM approach, PROZA is a self-contained theory and provides predictions for both intermolecular and intramolecular correlation simultaneously. While PRISM requires the application of numerical methods of solution, analytical solution of the Percus-Yevick (PY) and mean spherical approximation (MSA) versions of PROZA are available for most of the cases. Thus, within PROZA, one can have analytical expressions for both thermodynamics and structural properties. However, PRISM has an advantage of being applicable to a fluid of long-chain molecules, while so far, PROZA-based theories have been used mostly for fluids of relatively short-chain molecules. The application of the latter approach is restricted by the increasing dimensionality of the matrices involved into the corresponding multidensity OrnsteinZernike (OZ) equation, which appears to be proportional to the length of the chain molecules. We will attempt to remove this obstacle and reduce the dimensionality of the matrices by neglecting the chain-end effects, that is, assuming that all intermolecular site-site distribution functions are equal regardless of their position in the chain. This allows us to treat molecules of any length using a corresponding OZ equation of the same dimensionality. In the present study, the PROZA integral equation approach is applied to study the effect of the silica surface on the structuring of the polymer melt near a single surface and in a space between two surfaces. We also are making an attempt to reveal the effect of this polymer structuring on the effective interaction between filler particles by calculating the disjoining pressure exerted by the polymer on two surfaces. We start from PROZA approach for the bulk structure of PDMS melt and then extend it on the inhomogeneous PDMS melts including the presence of a single confining surface and a pair of surfaces, that is, a slit pore of a variable thickness. III. PROZA Multidensity Formalism for Polymer Melts A. Fluid of Linear Homonuclear Chain Molecules. We start with the simplest case of chain molecule fluids, a melt of linear homonuclear chain molecules. Each homonuclear chain molecule consists of m identical repeated monomer units. Each monomer unit is composed of a single atom modeled by a hard sphere with two interacting sites, A and B. Each of these sites bond one site of two different atoms from two neighboring monomers to form the monomer chain. The number m of repeated monomer units serves as a measure of the length of the chain. The multidensity OZ equation for such a model chain fluid reads

hˆ Rβ(k) ) cˆ Rβ(k) + F

∑γ (cˆ Rγ(k) + ∆ˆ Rγ(k))R(hˆ γβ(k) + ∆ ˆ γβ(k)) (9)

where F is the number density of chain molecules and subscripts R, β, and γ stand for the monomers in the system. The notations hˆ Rβ(k) and cˆ Rβ(k) are used for the matrices whose elements are the Fourier transforms of the elements of the matrices hRβ(r) and cRβ(r), respectively. In particular

()

15628 J. Phys. Chem. C, Vol. 111, No. 43, 2007

Henderson et al. them are associated with O and Si atoms, while another two are CH3 groups. This kind of modeling is called the united atom (UA) approach since the methyl group CH3 is treated as a single atom while O and Si atoms are treated explicitly. The full set of parameters of this model have been given in refs 41 and 42. To keep the theory treatable, the UA model has been further simplified in the following way: (i) removing bending and torsional interactions, (ii) adding one more oxygen atom so that the chain from one end starts with, say, a Si atom and ends up with, say, an O atom (for symmetry), and (iii) having, on both ends of the chain, two CH3 groups, instead of three as in the original model. 2. OZ Equation. The multidensity OZ equation for such a model has the following form

cR0β0(r)

cR0βA(r) cR0βB(r)

cRAβ0(r)

cRβ(r) ) cRAβA(r) cRAβB(r) cRBβ0(r)

cRBβA(r) cRBβB(r)

and similarly for hRβ(r). The functions hRmβn(r) and cRiβj(r) are the partial pair and direct correlation functions. The subscripts i and j characterize the bonding state of the site of the monomer R, namely, 0 means an unbound site, A means site A is bonded, and B means site B is bonded. The matrices and ∆ ˆ Rβ(k) are defined as follow

Rij ) 1 - δij + δ0iδ0j and

(

) (

0 0 0 0 0 0 0 ˆ (k) + δR,β-1 0 0 ∆ ˆ Rβ(k) ) δR,β+1 0 0 ∆ ∆ ˆ (k) 0 0 0 0 0

)

(10)

where δRβ is the Kronnekker symbol

∆ ˆ (k) )

1 sin(kL) F kL

with L being the bond length. Taking the sum over R and β from both sides of the OZ eq 9 and assuming that both the direct and total correlation functions are independent of the indices R and β, we have

hˆ (k) ) cˆ (k) + F(cˆ (k) + ∆ ˆ (k))R(hˆ (k) + ∆ ˆ (k))

(11)

where

hˆ (k) )

1 2

m

hˆ Rβ(k) ∑ Rβ

cˆ (k) )

1 m

(

2

0 0 0 ∆ ˆ (k) ∆ ˆ (k) ) 0 0 ˆ (k) 0 0 ∆

cˆ Rβ(k) ∑ Rβ

)

with

1 m - 1 sin(kL) ∆ ˆ (k) ) mF m kL

(

)

(12)

The dimensionality of the set of OZ equations (eq 11) is 3 × 3 regardless of the chain length m. B. Application to the Bulk PDMS Melts. All calculations that will be reported here have been performed for a model polymer system that mimics a polydimethylsiloxane (PDMS) melt. As for the filler prototype, we consider silica particles. 1. United Atom Model. The chemical formula of polydimethylsiloxane is (C2H6OSi)n. This means that the molecules of this polymer melt are modeled by the chains composed of the repeated monomer units that consist of four species; two of

4

ab ab hˆ Rβ (k) ) cˆ Rβ (k) + F

ac ac (cˆ Rγ (k) + ∆ ˆ Rγ (k))Rc(hˆ cb ∑ ∑ γβ(k) + c)1 γ

∆ ˆ cb γβ(k)) (13) Now all of the functions additionally have superscripts a, b, and c that denote the type of atoms comprising the chain monomer unit. Because of this, in contrast to eq 9, there is an extra summation over the species of monomer units where 1 stands for the O atom, 2 stands for the Si atom, and 3 and 4 both stand for the CH3 groups treated as a single atom. The chain length is m, which equals the a number of either oxygen or silicon atoms. Following the scheme developed for the fluid of linear homonuclear chain molecules and taking the sum over R and β from both sides of the OZ eq 13 we have

(

hˆ (k) ) cˆ (k) + F(cˆ (k) + ∆ ˆ (k))R(hˆ (k) + ∆ ˆ (k))

where hˆ (k) ) 11 11 12 12 12 12 h11 00 h0A h0B h00 h0A h0B h0C

13 13 14 14 h12 0D h00 h0A h00 h0A

h11 A0 h11 B0 h21 00 h21 A0 h21 B0 h21 C0 h21 D0 h31 00 h31 A0 h41 00 h41 A0

12 13 14 14 hAD h13 A0 hAA hA0 hAA

11 hAA h11 BA h21 0A 21 hAA h21 BA h21 CA 21 hDA h31 0A 31 hAA h41 0A 41 hAA

h11 AB h11 BB h21 0B h21 AB h21 BB h21 CB h21 DB h31 0B h31 AB h41 0B h41 AB

h12 A0 h12 B0 h22 00 h22 A0 h22 B0 h22 C0 h22 D0 h32 00 h32 A0 h42 00 h42 A0

12 hAA h12 BA h22 0A 22 hAA h22 BA h22 CA 22 hDA h32 0A 32 hAA h42 0A 42 hAA

h12 AB h12 BB h22 0B h22 AB h22 BB h22 CB h22 DB h32 0B h32 AB h42 0B h42 AB

h12 AC h12 BC h22 0C h22 AC h22 BC h22 CC h22 DC h32 0C h32 AC h42 0C h42 AC

)

(14)

13 13 14 14 h12 BD hB0 hBA hB0 hBA 23 23 24 24 h22 0D h00 h0A h00 h0A 22 23 24 24 hAD h23 A0 hAA hA0 hAA 23 23 24 24 h22 BD hB0 hBA hB0 hBA 23 23 24 24 h22 CD hC0 hCA hC0 hCA 22 23 24 24 hDD h23 D0 hDA hD0 hDA 33 33 34 34 h32 0D h00 h0A h00 h0A 32 33 34 34 hAD h33 A0 hAA hA0 hAA 43 43 44 44 h42 0D h00 h0A h00 h0A 42 43 44 44 hAD h43 A0 hAA hA0 hAA

(15)

The matrix for ∆ ˆ (k) has the same 12 × 12 dimensionality with all elements equal to zero except

ˆ 21 ˆ 12 ˆ 21 ∆ ˆ 12 AB(k) ) ∆ BA(k) ) ∆ BA(k) ) ∆ AB(k) )

∆ ˆ 23 ˆ 32 CA(k) ) ∆ AC(k) ) and

m - 1 sin(kL12) m mFkL12 (16)

sin(kL23) mFkL23

(17)

Integral Equation Study of Particle Confinement

∆ ˆ 24 ˆ 42 DA(k) ) ∆ AD(k) )

sin(kL24) mFkL24

J. Phys. Chem. C, Vol. 111, No. 43, 2007 15629

(18)

The chain length, m, is measured by the number of either oxygen or silicon atoms. The partial pair correlation functions hRβ(r) that result from the solution of eq 11 do not have an immediate physical meaning. The usual physical meaning has the atom-atom (the same as that in the case of the homonuclear chain) total pair correlation function or atom-atom total radial distribution function

gtotal ab (r) )

∑ij habij (r) + 1

(19)

where the sum is over the bonding states of the sites of the atoms a and b. 3. Closure Conditions. To solve OZ eq 14, we have used both the polymer PY (PPY) closure ab cab 00(r) ) fab(r)(t00(r) + 1)

(20)

ab ab ab cab 0K(r) ) fab(r)t0K(r) cK0(r) ) fab(r)tK0(r)

(21)

ab cab MK(r) ) fab(r)tMK(r)

(22)

and the polymer HNC (PHNC) closure ab ab cab 00(r) ) g00(r) - t00(r) - 1

(23)

ab ab ab ab ab cab 0K(r) ) (g00(r) - 1)t0K cK0(r) ) (g00(r) - 1)tK0 (24) ab ab ab ab ab cab MK(r) ) g00(r)(tM0(r)t0K(r) - tMK(r)) - tMK(r)

fab(r) ) exp(-Uab(r)/kT) - 1

(26)

ab ab ab (r) ) hRβ (r) - cRβ (r) tRβ

(27)

) exp(-Uab(r)/kT +

tab 00(r))

(28)

and Uab(r) is the nonbonded Lennard-Jones potential of the UA model of the PDMS melt. C. Application to the Confined PDMS Melts. Within the framework of the PROZA integral equation theory, eq 3 for the polymer partial density profiles has the form 4

yai (z) ) δi0 + 2πF

∑ ∑ ∑k b)1 j

z+(1/2)[σ +σ ] (1/2)[σ +σ ] b rcab ∫z-(1/2)[σ ij (r)dr)Rjkhk (z1)dz1 +σ ] ( ∫|z -z| a

a

b

b

a

b

hai (z) + δi0 ) yai (z) exp[-βUa(z)]

(29)

1

where yai (z) is the partial polymer-surface cavity distribution function. The coefficients Rjk in eq 29 are the elements of the normalized density matrix while the subscripts i, j, and k run over the bonding states of the sites of the O, Si, and CH3 atoms. The most important ingredients of eq 29 are the partial direct correlation functions cab ij that characterize the bulk polymer system at a fixed density F. These functions are the only input for the inhomogeneous system calculations and have been discussed in a previous section, together with the partial pair correlation functions hab ij .

(30)

The function Ua(z) in eq 30 represents the potential of the polymer-surface interaction. There are a few different ways to characterize this interaction. In the present calculations, we have used the polymer-surface interaction that is representative of a surface containing oxygen atoms to characterize an amorphous silica surface. The corresponding potential functions have the form43

[ ( ) ( )] [ ( ) ( )]

Ua(z) ) Aaw ) Aaw

2 σaw 7 3 σaw 7 z 4 z

4

2 σaw 5 z

4

10

-

σaw z

for a ) Si, O for a ) CH3

(31)

where the coefficient Aaw ) 2πFwsσ2awaw. This potential has been obtained by integrating over the surface plane that is assumed to consist of a given density of oxygen atoms, representing the OH groups on the surface. The surface density of OH groups was taken to be Fws ) 0.039 Å-2. This number was obtained in an atomistic simulation of amorphous silica.51 The parameters of the wall-atom interactions are chosen to be the same as those in the hybrid/UA model for the bulk PDMS, assuming that σaw ) σaO and aw ) aO. The corresponding normalized local density profiles ga(z) ) Fa(z)/Fbulk of the PDMS atoms, which are the subject of a interest, can be written as follows

ga(z) ) 1 +

(25)

where

gab 00(r)

By applying the PPY approximation to treat the polymersurface correlations, the relation between yai (z) and partial polymer-surface correlation functions hai (z) reads

∑i hai (z)

(32)

where the sum is over the bonding states of the sites of the atom a. IV. Results The solutions of integral equations were obtained at room temperature T ) 300 K and a number density of PDMS molecules of F ) 0.0004 Å-3. The difference between results obtained within the PPY closure and the PHNC closure is negligible. Therefore, only PROZA/PHNC results are presented and discussed. A. Radial Distribution Functions of a Bulk PDMS Melt. In Figure 1, the total atom-atom radial distribution functions gtotal ab (r), which consist of both intramolecular and intermolecular contributions to the atom-atom correlations, are displayed. The symbols correspond to our MD simulation data for an original UA model. According to the parameters of the UA model of the PDMS melt, the bonding distances between methyl groups and silicon atom, L23 and L24, and especially between oxygen and silicon atoms, L12, are rather short. With this choice of the model parameters, the oxygen atom is completely buried inside of the repulsive core of the two neighboring silica atoms. This creates some problems for PROZA theory since its predictions (as well as predictions of the so-called “proper” theories44,45) for the rigid chains and site-site molecules with short interatomic bonding distances become less reliable.46-50 As a result, the structural properties of a four-atom UA model of the PDMS melt, which follows from the PROZA theory, seems to be not very accurate when compared against simulation data. In contrast, the PRISM theory, which uses the intra-

15630 J. Phys. Chem. C, Vol. 111, No. 43, 2007

Henderson et al.

Figure 1. Total radial distribution functions gtotal ab (r) of the PDMS melt that results from MD simulations (open circles) and are calculated from PROZA/PHNC theory (solid lines), both using the united atom model.

molecular part gintra ab (r) of the atom-atom radial distribution function as an input, seems to provide a better description of the PDMS structure, as follows from the results reported by Curro et al.24-28 However, the subject of our interest in the description of a bulk polymer is the intermolecular correlation functions ginter ab (r) that describe correlations between atoms belonging to different chain molecules. To exclude the source of the abovediscussed problems for the PROZA theory, we simplified further the original UA model of the PDMS molecule by removing the oxygen site and considering the chain with neighboring Si atoms connected directly. Figure 2 shows the results for intermolecular radial distribution functions that have been calculated from the PROZA theory using this, as we called the simplified UA model, when the presence of oxygen atoms was counted in an effective way. As can be expected, the simplified version of the UA model becomes more flexible, which immediately results in the atom-atom intermolecular radial distribution functions ginter ab (r) that are less structured. However, we still can see that theoretical curves reproduce general features of the “experimental” radial distribution functions obtained from MD simulations of the original hybrid/UA model but fail to reproduce the details. B. Local Density of a PDMS Melt near a Single Silica Surface. The normalized local density distributions ga(z) of silicon and carbon atoms next to a single silica surface that follow from the PROZA theory calculations using the simplified UA model are presented on the bottom part of Figure 3. The upper part of the same Figure 3 shows results that have been

obtained from MD simulations of the full hybrid/UA model of the PDMS melt next to an atomistic silica surface. Both simulations and theory indicate that carbon atoms show a rather strong tendency to form layers next to the surface; at the same time, layering within the pendant methyl groups is weaker while these groups seems to be closer to the silica surface than the silicon atoms on the backbone. This might be expected since the methyl groups are more exposed and tend to shield the polymer backbone atoms. We can see that interface between the bulk PDMS and silica surface consists of densely packed and partly ordered layers in which the polymer segments tend to run preferentially parallel to the solid surface. This perturbation of the polymer local density and polymer local ordering does not extend into the bulk phase by more than three times the transverse diameter of the PDMS chain, that is, approximately 20 Å. C. PDMS between Two Surfaces. The normalized local density distributions of silicon and carbon atoms across a slit formed by two surfaces each representing a filler particle at different separations are shown in Figure 4. The calculations have been performed in a way that the density of the polymer in the slit is not fixed a priori but is allowed to vary with the slit thickness in such a way that the contents of the slit are always in equilibrium with the bulk PDMS that is kept at fixed conditions. Figure 4a shows the case when the slit thickness allows only for one polymer layer to be accommodated in a slit. Whenever more than one layer is present, the layers are not layers of single PDMS molecules, in the sense that monomers of a given

Integral Equation Study of Particle Confinement

J. Phys. Chem. C, Vol. 111, No. 43, 2007 15631 as shown schematically in Figure 3. The pendant methyl groups of the PDMS chains are closer to the filler surface than the silicon atoms on the backbone. Again, this seems to be a physically correct result since the methyl groups are more exposed and tend to shield the polymer backbone atoms. Once density profiles are known, some other properties of the polymer melt can be calculated. In particular, the pressure Π exerted by a polymer melt on the confining surfaces as a function of the surface separation H can be calculated from the density profiles by using eq 7. To satisfy the self-consistency of the numerical procedure, the pressure of the bulk polymer melt is defined as PB ) PN(H f ∞). Particularly, to calculate bulk pressure, we used a surface-to-surface separation of H ) 50 nm, which, according to Figure 3, provides a well-defined homogeneous region. The pressure Π(H) is shown in Figure 5. We can see that there is a large repulsive barrier at the separation between surfaces around 10 Å. If we look for the density distribution at that separation, we find that this gap allows for a single polymer layer. V. Discussions and Conclusions

Figure 2. Intermolecular part of the total radial distribution functions ginter ab (r) that results from MD simulations (open circles) of the united atom model and are calculated from PROZA/PHNC theory (solid lines) using a simplified united atom model.

TABLE 1: Estimate of the Average Separation between Filler Particles in a Polymer/Particle Mixture particle diameter D, nm

volume fraction φ

weight fraction

separation, eq 33 〈δ〉, nm

10 10 10 10 10

6.1 9.9 15.7 25 35

9.7 15.3 23.5

13 10 6.8 4.4 2.8

polymer molecule participate in several layers and bridge them. It is clearly seen for the slit with two and three layers where the polymer density at the middle of the slit is not zero. From these results, we can see that silicon atoms are responsible for polymer stratification, while carbon groups result in bridging effects. In particular, we can see that the local density distribution indicates a very high probability for the polymer to be adsorbed on the silica surface by means of a single CH3 methyl group while the second CH3 group is placed away from the surface,

One further aspect that should be taken into account concerns the aggregation of silica particles embedded into the PDMS melt. Quite recently,52 it has been shown experimentally that the aggregates of nanometric silica spheres have a strong influence on the mechanical properties of the latex/silica nanocomposite. Usually, it is assumed that the clustering of filler particles occurs because of strong attractive van der Waals forces that are always present between particles at short distances. By applying the Derjaguin approximation (eq 8), we can see that it is very probable that, for the case of the PDMS/silica system, the aggregation (vdW aggregation) will occur only if no polymer molecules will be in the gap. In all other cases, there is a stabilized repulsive barrier due to the adsorbed polymer layers. Another possibility for the vdW aggregation can be expected when the polymer structuring (layering) near a particle surface will be diminished by tuning the polymer-particle interaction (eq 31). Indeed, it was observed52 that the average aggregation number to be tunable via the precursor solution pH may in fact influence the polymer density distribution in the particle’s neighborhood.1 On other hand, the average spacing 〈δ〉 between spherical silica particles of diameter D as a function of volume fraction φ ) (π/6)FsD3 occupied by filler particles can be estimated from

〈δ〉 ) D

[( ) ] 1 πx2 6 φ

1/3

-1

(33)

It follows that, at the same loading (the same volume fraction φ), particles with a larger diameter D will be separated by larger gaps δ. For example, in the case of a particle diameter of D ) 10 nm and a particle volume fraction of φ ) 0.35, the average gap width δ between particles is around 30 Å, while in the case of D ) 100 nm, the gap width increases up to 300 Å. As follows from our local density calculations (Figures 3, 4, and 6), the surface separations around 40 Å are crucial to produce the confinement effect on the surrounding PDMS polymer. Thus, the separation between particles and the conditions of the geometric confinement in the gap between them can be one of the important ingredients for reinforcement to occur. Finally, similar to the monatomic molecular liquids, polymer chains next to the confining surfaces form layers. The layers in molecular liquids are formed by the individual molecules, that is, one molecule could participate in the formation of only one layer. In contrast, the PDMS polymer molecule has an ability to participate in few neighboring layers. Because of this, such layered polymer molecules could be associated with a tightly

15632 J. Phys. Chem. C, Vol. 111, No. 43, 2007

Henderson et al.

Figure 3. Normalized local density of the PDMS melt next to a single silica surface that results from the MD simulations (upper part) of the full hybrid/UA model and are calculated from the PROZA/PHNC theory (lower part) using a simplified UA model. The right-hand side of the lower part shows a schematic interpretation of the density profiles.

Figure 4. Normalized local density of silicon and carbon atoms of the PDMS melt in the four different slits formed by a pair of silica surfaces as calculated from the PROZA/PHNC theory. The notation for the lines is the same as that used in Figure 3.

Integral Equation Study of Particle Confinement

Figure 5. The solvation pressure exerted by a simplified united atom model PDMS melt on silica surfaces as a function of the separation between surfaces.

Figure 6. Normalized local density of silicon and carbon atoms of the PDMS melt in a wide slit formed by the silica surfaces. The coexistence of what can be called a tightly bound (region with density oscillations) and a loosely bound (region without oscillations) polymer can be seen. The notation for the lines is the same as that in Figure 3.

bound polymer. For the same reason, a polymer in the region that follows a layered region still does not represent the bulk polymer since some of the chains from this region are part of the bound polymer; these polymer chains can be associated with a loosely bound polymer (see Figures 2 and 6). The existence of two such regions in confined polymer melts has been suggested by Tsagaropoulos and Eisenberg53 in their model to explain experimental data on two glass transitions in the polymer/particle systems. According to their measurements, a tightly bound polymer has a very low spin-spin relaxation time (T2) that does not change significantly with temperature. This has been explained by the fact that the tightly bound polymer does not participate in either of the two glass transitions. On the other hand, a polymer beyond this layer exhibits a usual glass transition, but the T2 value is notably lower than the T2 of the unfilled polymer. This observation is attributed to the existence of a restricted mobility layer. The estimated thickness of the interfacial polymer layer reported in the literature in experimental studies of PDMS/silica systems is 5 nm using neutron scattering,54 0.8 nm using NMR,55 and 1-25 nm using dielectric relaxation spectroscopy (DRS).56 If the interfacial layer is defined as a region where local polymer density oscillates, then both of our earlier simulation data1 and the present integral equation results are within the range of observations. Acknowledgment. The work of the BYU authors was supported by a grant from DOE through the Lawrence Livermore Laboratory. The authors are grateful for a generous allotment of computer time by the Marylou Fulton Supercomputing Center at BYU. References and Notes (1) Gee, R. H.; Maxwell, R. S.; Balazs, B. Polymer 2004, 45, 3885. (2) Landel, R. F. Trans. Soc. Rheol. 1958, 2, 53.

J. Phys. Chem. C, Vol. 111, No. 43, 2007 15633 (3) Landel, R. F.; Smith, T. L. J. Am. Rocket Soc. 1961, 31, 599. (4) Vacatello, M. Macromolecules 2001, 34, 1946. (5) Vacatello, M. Macromolecules 2002, 35, 8191. (6) Starr, F. W.; Shroeder, T. B.; Glotzer, S. C. Macromolecules 2002, 35, 4481. (7) Starr, F. W.; Douglas, J. F. J. Chem. Phys. 2003, 119, 1777. (8) Yethiraj, A.; Dickman, R.; Hall, C. K. J. Colloid Interface Sci. 1992, 151, 102. (9) Khalatur, P. G.; Zherenkova, L. V.; Khokhlov, A. R. Physica A 1997, 247, 205. (10) Hooper, J. B.; Schweizer, K. S. Macromolecules 2005, 38, 8858. (11) Patel, N.; Egorov, S. N. J. Chem. Phys. 2004, 121, 4987. (12) Chandler, D.; Andersen, H. C. J. Chem. Phys. 1972, 57, 1930. (13) Wertheim, M. J. Stat. Phys. 1984, 35, 35. (14) Wertheim, M. J. Stat. Phys. 1986, 42, 477. (15) Ornstein, L. S.; Zernike, F. Proc. Acad. Sci. Amsterdam 1914, 17, 793. (16) Given, J. A.; Stell, G. Physica 1994, 209, 495. (17) Trokhymchuk, A.; Pizio, O.; Holovko, M. Sokołowski, S. J. Chem. Phys. 1997, 106, 200. (18) Anguiano Orozco, G.; Pizio, O.; Sokołowski, S.; Trokhymchuk, A. Mol. Phys. 1997, 91, 625. (19) Pizio, O.; Trokhymchuk, A.; Henderson, D.; Labik, S. J. Colloid Interface Sci. 1997, 191, 86. (20) Walley, K. P.; Schweizer, K. S.; Peanasky, J.; Granick, S. J. Chem. Phys. 1994, 100, 3361. (21) Henderson, D.; Blum, L.; Lebowitz, J. L. J. Electroanal. Chem. 1979, 102, 315. (22) Rickayzen, G. Mol. Phys. 1985, 55, 161. (23) Derjaguin, B. V. Kolloid Zh. 1934, 69, 155. (24) Curro, J. G.; Schweizer, K. S. Macromolecules 1987, 20, 1928. (25) Curro, J. G.; Schweizer, K. S. J. Chem. Phys. 1987, 87, 1842. (26) Schweizer, K. S.; Curro, J. G. AdV. Polym. Sci. 1994, 116, 321. (27) Schweizer, K. S.; Curro, J. G. AdVances in Chemical Physics; Interscience, New York, 1997; Vol. 98. (28) Curro, J. G. AdV. Polym. Sci. 2005, 173, 209. (29) Wertheim, M. S. J. Chem. Phys. 1987, 87 7323. (30) Chang, J.; Sandler, S. I. J. Chem. Phys. 1995, 102, 437. (31) Yu Kalyuzhnyi, V.; Cummings, P. T. J. Chem. Phys. 1995, 103, 3265. (32) Lin, C.-T.; Kalyuzhnyi, Yu. V.; Stell, G. J. Chem. Phys. 1998, 108, 6513. (33) Kalyuzhnyi, Yu. V.; Lin, C.-T.; Stell, G. J. Chem. Phys. 1998, 108, 6525. (34) Stell, G.; Lin, C.-T.; Kalyuzhnyi, Yu. V. J. Chem. Phys. 1999, 110, 5444. (35) Lin, C.-T.; Stell, G.; Kalyuzhnyi, Yu. V. J. Chem. Phys. 2000, 112, 3071. (36) Kalyuzhnyi, Yu. V. Condens. Matter Phys. 1997, 11, 71. (37) Kalyuzhnyi, Yu. V.; Holovko, M. F. Condens. Matter Phys. 2002, 5, 211. (38) Kalyuzhnyi, Yu. V.; McCabe, C.; Cummings, P. T. Mol. Phys. 2002, 100, 2499. (39) Kalyuzhnyi, Yu. V.; McCabe, C.; Whitebay, E.; Cummings, P. T. J. Chem. Phys. 2004, 121, 8128. (40) Kalyuzhnyi, Yu. V.; Cummings, P. T. J. Chem. Phys. 2001, 115, 540 (41) Sides, S. W.; Curro, J.; Grest, G. S.; Stevens, M.; Soddemann, T.; Habenschuss, A.; Londono, J. D. Macromolecules 2002, 35, 6455. (42) Frischknecht, A. L.; Curro, J. G. Macromolecules 2003, 36, 2122. (43) Nath, S. K.; Frischknecht, A. L.; Curro, J. G.; McCoy, J. D. Macromolecules 2005, 38, 8562. (44) Chandler, D.; Silbey, R. Mol. Phys. 1982, 46, 1335. (45) Rossky, P. J.; Chiles, R. A. Mol. Phys. 1984, 51, 661. (46) Kalyuzhnyi, Yu. V.; Stell, G.; Llano-Restrepo, M. L.; Chapman, W. G.; Holovko, M. F. J. Chem. Phys. 1994, 101, 7939. (47) Yethiraj, A.; Dickman, R.; Szamel, G.; Kierlik, E.; Rosinberg, M. L. Mol. Phys. 1994, 82, 937. (48) Yethiraj, A. Mol. Phys. 1994, 82, 957. (49) Attard, P. Mol. Phys. 1994, 83, 273. (50) Kalyuzhnyi, Yu. V. Mol. Phys. 1999, 96, 1289. (51) Tsige, M.; Soddemann, T.; Rempe, S. B.; Grest, G. S.; Kress, G. D.; Robbins, M. O.; Sides, S. W.; Stevens, M. J. Webb, E. III. J. Chem. Phys. 2003, 118, 5132. (52) Oberdisse, J. Soft Matter 2006, 2, 29. (53) (a) Tsagaropoulos, G.; Eisenberg, A. Macromolecules 1995, 28, 396. (b) Tsagaropoulos, G.; Eisenberg, A. Macromolecules 1995, 28, 6067. (54) Arrighi, V.; Higgins, J.; Burgess, A.; Floudas, G. Polymer 1998, 39, 6369. (55) Litvinov, V.; Spiess, H. Macromol. Chem. Phys. 1991, 192, 3005. (56) Kirst, K.; Kremer, F.; Litviniv, V. Macromolecules 1993, 26, 975.