Article pubs.acs.org/JPCB
Entropic Interactions in Semiflexible Polymer Nanocomposite Melts Yangwei Jiang,† Dong Zhang,† Linli He,‡ and Linxi Zhang*,† †
Department of Physics, Zhejiang University, Hangzhou 310027, China Department of Physics, Wenzhou University, Wenzhou 325035, China
‡
S Supporting Information *
ABSTRACT: By employing molecular dynamics simulations, we explored the effective depletion zone for nanoparticles (NP) immersed in semiflexible polymer melts and calculated the entropic depletion interactions between a pair of NPs in semiflexible polymer nanocomposite melts. The average depletion zone volumes rely mainly on polymer chain stiffness and increase with chain stiffness increasing. In the semiflexible polymer nanocomposite melts, the entropic depletion interactions are attractive and anisotropic, and increase with chain stiffness increasing. Meanwhile, the attractive interactions between NPs and polymers can also affect strongly the entropic depletion interactions. For the semiflexible polymer nanocomposite melts in the athermal system, the entropic depletion interactions change from anisotropic to isotropic when the NP/polymer interactions increase. For NPs in the rodlike polymer melts, a mixture structure of contact/“bridging” aggregations for NPs is formed at a strong attractive NP/polymer interaction. Our calculations can provide an effective framework to predict the morphology of NPs immersed in semiflexible polymer melts.
1. INTRODUCTION Materials that consist of mixtures of polymers and organic/ inorganic particles, ranging from nanometers to microns, are used in a wide variety of applications, and the inherent macroscopic properties (such as mechanic, optical, electronic and so on) of polymer nanocomposites depend strongly on the microscopic morphology of constituent nanoparticles (NPs) in the polymer matrix.1 Interparticle effective potentials, including the entropic interactions associated with the mismatch in species of different size and geometry and the direct enthalpic interactions such as van der Waals attractions and Coulomb forces, play a key role in understanding in this concerned issue. Over the past few decades, there have been tremendous researches from experiments, theories, and computer simulations focused on the entropic interactions in flexible polymer nanocomposites1−17 or NP/rod mixtures.18−22 However, to the best of our knowledge the study on entropic interactions between NPs in semiflexible polymer nanocomposite melts is rare. As an important property, polymer stiffness, undoubtedly, has a crucial impact on the entropic interactions between NPs in semiflexible polymer nanocomposites. For instance, the change of the persistence length will result in an isotropicnematic transition for semiflexible polymer melts23 and this transition is expected to affect the entropic interaction. A study on the entropic interaction between NPs in semiflexible polymer nanocomposite can connect the gap between the flexible polymer nanocomposite and NP/rod mixture and provide a comprehensive view to understand the behavior of NPs in polymer matrix with various polymer stiffnesses, such as flexible, semiflexible, and the rodlike limit. For dilute or semidilute flexible polymer concentrations, the entropic interaction between adventive NPs is strongly related © 2015 American Chemical Society
to a dimensionless length scale ratio, namely, R/Rg for the dilute concentration2 and R/ξ for the semidilute concentration,3 where R is the radius of NP, Rg is the radius of gyration of flexible polymer chains, and ξ is the collective mesh or density−density screening length or correlation length.2,3 For the NP/rod mixture with the rod concentration below the Onsager concentration, this key dimensionless length scale ratio becomes R/L with L indicating the length of rod.18−20 However, for polymer melts or the concentrated solutions, the physical picture for the entropic interaction is very different and more complex. For the flexible polymers, some previous studies suggested that the entropic interaction plays an important role in high-polymer concentrations associated with local monomer packing correlations and the ratio D/σ (D and σ are the diameters of NPs and polymer monomers, respectively), while the range of entropic interaction may be controlled by the density fluctuation correlation length.14 That is to say, in the athermal limit there is a polymer density fluctuation around an NP because polymer configurational freedom is restricted in the vicinity of an NP. Thus, each NP is surrounded by a “depletion zone” in which the local polymer density is less or greater than the bulk value. When the depletion zones of two NPs overlap, entropic interactions will be observed. If the local polymer density is less than that in bulk, the entropic interaction is attractive; otherwise, it is repulsive. Many polymer molecules have internal stiffness and even for polymers considered flexible, like polyethylene, the Rouse model fails as soon as the local chemical structure can no longer Received: September 30, 2015 Revised: November 29, 2015 Published: December 31, 2015 572
DOI: 10.1021/acs.jpcb.5b09551 J. Phys. Chem. B 2016, 120, 572−582
Article
The Journal of Physical Chemistry B be neglected.24 One of the effects of this structure is a certain stiffness in the polymer chain due to the valence angles of the polymer backbone.24 Thus, it is interesting to investigate the aggregation structures and effective interactions of NPs in semiflexible polymer nanocomposite melts. Because semiflexible polymers are considerably more difficult to treat theoretically as they require the fulfillment of additional constraints such as keeping the total chain contour length fixed, in this paper we employed MD simulations to study the entropic interactions between a pair of NPs in semiflexible polymer nanocomposite melts. We first determined the effective depletion zone for NPs in semiflexible polymer nanocomposite melts and calculated the entropic potentials between a pair of NPs with various polymer stiffnesses and polymer concentrations. Unlike in flexible polymer melts or in dilute concentrations, the effective depletion zones for NPs in semiflexible polymer nanocomposite melts are anisotropic. According to our simulations, there are at least three different predictable behaviors on the profiles of entropic interactions according to different polymer stiffness. Meanwhile, we also investigated the effects of NP/polymer interactions on the entropic interactions and aggregation structures of NPs in semiflexible polymer nanocomposite melts. In the rodlike nanocomposites, the entropic interactions transform from anisotropic to isotropic when the NP/polymer interactions increase, which rely mainly on the competition between entropy and enthalpy of the system.
NPs are modeled as Lennard-Jones spheres of diameter D. Mass densities of polymer beads and the NPs are the same, and D = 5σ. The NP/NP and NP/polymer interactions are represented by truncated and shifted LJ potential of the form ⎧ ∞ r ≤ REV ⎪ 12 6⎤ ⎡ ⎪ ⎛ σ ⎞ ⎪ ⎛ σ ⎞ ULJ(r ) = ⎨ 4ε⎢⎜ ⎟ −⎜ ⎟ ⎥ + ULJ(rc) REV < r < rc ⎢ r − R ⎝ ⎠ ⎝ r − REV ⎠ ⎥⎦ ⎪ ⎣ EV ⎪ ⎪ 0 r ≥ rc ⎩
(4)
The purely repulsive NP/NP interactions are described using the Weeks−Chandler−Anderson (WCA) potential given by ε = kBT, rc = REV + 21/6σ, and REV = 4σ.27 NP/polymer interactions are varied from repulsive to strong attractive. The repulsive NP/polymer interactions are described by the WCA potential with ε = kBT, rc = REV + 21/6σ, and REV = 2σ. We refer to this system as ULJ(r) = WCA. Meanwhile, the attractive NP/ polymer interactions are also considered, and we refer to this system as ULJ(r) = LJ. Similar to previous literature,16 the attractive NP/polymer potentials are given by a shifted and truncated LJ potential with rc = REV + 2.5σ and REV = 2σ, and the interaction strengths are varied from very weak attractive εnp = 0.1kBT to a strong attractive interaction εnp = 9.0kBT. MD simulations were carried out in the NVT ensemble by using the open source software, LAMMPS,28 at the reduced temperature T* = kBT/ε = 1.0, and the simulation box is a periodic cubic cell with the size of (30σ)3. For the maximum polymer density, ρ0σ3 = 0.8, the total number of chains is 720. Following the work of Hooper et al.,29 the polymer density ranges from ρ0σ3 = 0.4 to 0.8 by reducing the number of polymer chains in the simulation boxes of the same box size, and the selection of densities is in the melt regime.30 Average bond orientations are chosen to be almost parallel to the z-axis direction in the rodlike polymers. Each simulation ran long time enough to reach the equilibrium state and all data were sampled far from the polymer chain relaxation time. In general, the polymermediated potential of mean force (PMF) U(r) between a pair of NPs can be calculated directly from simulations of the semiflexible polymer nanocomposites containing two NPs and is given by16
2. MODEL AND SIMULATION METHOD Molecular dynamics (MD) method was performed on semiflexible polymer/NP composite melts. A standard beadspring model is used to model polymer chains,25 and each polymer chain consists of 30 monomers with a bead diameter of σ and a mass of m. The nonbonded interaction between all the polymer monomers is represented by a purely repulsive truncated and shifted Lennard-Jones (LJ) potential ⎧ ⎡⎛ ⎞12 ⎛ ⎞6 ⎤ σ σ 1 ⎪ r < 21/6σ ⎪ 4ε⎢⎜ ⎟ − ⎜ ⎟ ⎥ + ⎝ ⎠ ⎝ ⎠ r ⎦ 4 Ubead(r ) = ⎨ ⎣ r ⎪ ⎪ ⎩ 0 r ≥ 21/6σ
(1)
U (r ) = −kBT ln p0 (r )
where r is the distance between the ith and jth beads of polymer chains, and ε = kBT. All neighboring beads in polymer chains interact with the well-known finitely extensible nonlinear elastic (FENE) potential26 ⎡ ⎛ r ⎞2 ⎤ kR 02 ⎢ UFENE(r ) = − ln 1 − ⎜ ⎟ ⎥ , ⎢⎣ 2 ⎝ R 0 ⎠ ⎥⎦
where p0(r) is the probability of finding NPs separated by distance r during simulations. In fact, p0(r) is quite small for separations of interest, precluding sufficiently accurate sampling during a standard MD simulation run. Here we used the combining umbrella sampling31 with weighted histogram analysis method32 to obtain p0(r). In the umbrella sampling, the biasing potential was a simple harmonic spring
r < R0 (2)
Uib =
where r is the distance between two neighboring beads, k = 30ε/σ2, and R0 = 1.5σ.26 The average length between the bonded monomers is 0.98σ in our model. The stiffness of a polymer chain is described by the angle bending potential Ubending = b(1 + cos θ )
(5)
1 k(r − R i)2 2
(6)
Here k = 4kBT/σ2. For each semiflexible polymer nanocomposite, eight umbrella sampling simulations were run with Ri = 5.0σ, 6.0σ, ..., 12.0σ.
3. RESULTS AND DISCUSSION 3.1. Entropic Interactions between NPs in Semiflexible Polymer Nanocomposite Melts. In this section, the purely repulsive NP/polymer interaction of WCA is adopted. The number bead density of polymers ranges from
(3)
where θ is the angle between two consecutive bonds and the stiffness of polymer chains is controlled by varying the value of b. 573
DOI: 10.1021/acs.jpcb.5b09551 J. Phys. Chem. B 2016, 120, 572−582
Article
The Journal of Physical Chemistry B
NPs in semiflexible polymer nanocomposites with different polymer stiffness at ρ0σ3 = 0.6 are shown in Figure 1. In the case of flexible polymer chains, our results are in agreement with previous studies9,16−18 as shown in Figure 1a. The entropic interactions in the flexible polymer nanocomposite melts are oscillatory with attractive and repulsive components. The loss of the configurational freedom in the vicinity of the NP reaches the minimum, and the contact attraction (U(r)) of entropic interaction reaches the minimum at r ≈ D. For semiflexible polymer chains with the moderate stiffness, as shown in Figure 1b, the contact attraction and the range of entropic interactions increase monotonically with the polymer stiffness increasing. For the cases of b = 30−80kBT, a theoretical treatment of the entropic interaction is very difficult. Empirically, our simulation data can be well fitted by the relation U(x) = −A(1 − x/B)3 with x = (r − D)/σ (labeled by the solid curves in Figure 1b). This relation shares the same form with the well-known Derjaguin model, which can well describe the entropic interaction between two spheres of radius R immersed in a dilute suspension of thin rods of length Le when R ≫ Le and be written as18,19 3 ⎛ π x⎞ UD(x) = − kBTnrRLe2⎜1 − ⎟ 6 Le ⎠ ⎝
(7)
where nr is the equivalent number density of the rods, and x = (r − D)/σ. For R ∼ Le, there are fewer excluded orientations for a rod near a sphere than for a rod near a plate, and the Derjaguin model may overestimate the entropy gain from overlapped excluded volume. In general, for semiflexible polymers, the polymer chains within the depletion zone can be regarded as rods. According to the curve fittings in the Derjauin model, the equivalent number density of the rods nr is less than that of the polymer chains. When the polymer stiffness or the polymer bead number density increases, the scope of the depletion zone becomes large correspondingly and then the equivalent length Le of these rods also increases, see Table 1. Both the polymer stiffness and polymer bead number density can strengthen the contact attractions of NPs in semiflexible polymer nanocomposite melts. Meanwhile, the entropic interactions in rodlike polymer nanocomposite melts are shown in Figure 1c. A previous study of depletion interaction in a nematic of a hard, rodlike particle suggested that both the strength and range of the interaction crucially depend on the configuration of the spheres relative to the nematic director.33 For the case of parallel alignment, it gives the depletion potential as
Figure 1. Entropic depletion potentials U(r) between a pair of NPs in flexible (a), semiflexible (b), and rodlike (c) polymer nanocomposites with different chain stiffness at a polymer bead number density ρ0σ3 = 0.6. The open circles represent simulation data, and the solid lines are fits to the Derjaguin model for semiflexible polymer nanocomposites (b) and to the nematic model for rodlike polymer nanocomposites (c).
ρ0σ3 = 0.4 to ρ0σ3 = 0.8, and the polymer stiffness (i.e., bending energy) ranges from b = 0 to b = 400kBT. Considering the contour length of polymer chain is about Lc = 29σ, we can roughly divide the polymer stiffness into three different regions: flexible (b < 10kBT), semiflexible (b = 30−80kBT), and rodlike (b ≥ 150kBT), see Supporting Information. Correspondingly, the behaviors of entropic interactions in these three regions are different, and the entropic interactions U(r) between a pair of
⎛ x⎞ U (x) = −3πkBTnrR2Le⎜1 − ⎟ Le ⎠ ⎝
(8)
where x = (r − D)/σ and we refer to this form as the nematic model. For the case of perpendicular alignment, it declares that the depletion potential is still monotonic and its strength is weaker than that of parallel alignment. As shown in Figure 1c,
Table 1. Number Density nr and the Length Le of Rods for the Derjaguin Model in Figures 1b and 2a ρ0σ3 = 0.6
nr
Le
b = 30 kBT
b b b b
± ± ± ±
± ± ± ±
ρ0σ ρ0σ3 ρ0σ3 ρ0σ3
= = = =
30 50 70 80
kBT kBT kBT kBT
0.34 0.35 0.38 0.38
0.02 0.03 0.01 0.03
4.89 5.54 5.86 6.32
3
0.11 0.12 0.14 0.09 574
= = = =
0.4 0.5 0.6 0.7
nr 0.28 0.31 0.34 0.35
± ± ± ±
Le 0.03 0.02 0.02 0.01
2.88 4.29 4.89 5.59
± ± ± ±
0.08 0.11 0.11 0.13
DOI: 10.1021/acs.jpcb.5b09551 J. Phys. Chem. B 2016, 120, 572−582
Article
The Journal of Physical Chemistry B Table 2. Number Density nr and the Length Le of Rods for the Nematic Model in Figures 1c and 2b ρ0σ3 = 0.6 b b b b
= = = =
150 200 300 400
kBT kBT kBT kBT
nr 0.16 0.16 0.15 0.14
± ± ± ±
Le 0.01 0.01 0.01 0.01
2.74 3.39 4.14 4.72
± ± ± ±
b = 30 kBT ρ0σ3 ρ0σ3 ρ0σ3 ρ0σ3
0.12 0.05 0.06 0.05
= = = =
0.4 0.5 0.6 0.7
nr 0.05 0.09 0.15 0.22
± ± ± ±
Le 0.01 0.01 0.01 0.02
3.76 3.96 4.14 4.65
± ± ± ±
0.02 0.05 0.06 0.09
Figure 2. Entropic depletion potentials U(r) between a pair of NPs in polymer nanocomposites with various polymer bead number densities ρ0σ3 for different bending energies b = 30kBT (a) and b = 300kBT (b). The open circles represent simulation data, and the solid lines are fits to the Derjaguin model for the bending energy b = 30kBT (a), and to the nematic model for the bending energy b = 300kBT (b).
our data agrees well with the nematic model in the short-range, and the values of nr and Le are given in Table 2. However, at a large NP separation our simulation data is smaller than that in the nematic model and the ranges of entropic interactions are also larger than that in the nematic model. It is suggested that this deviation may be caused by the polymer stiffness or rod flexibility. Rodlike chains undulate in long-range and possess more configurational entropy than straight rods. Meanwhile, the mean end-to-end length of flexible rods is smaller than their contour length. Because of these two possible factors, the nematic model may underestimate the loss of configurational entropy of flexible rods in the vicinity of the depletion zone. On the other hand, the influence of polymer number bead density in semiflexible and rodlike polymer nanocomposite melts is shown in Figure 2. When the polymer density increases from ρ0σ3 = 0.4 to 0.8, the contact attractions and the range of entropic interactions also increase monotonically for NPs in semiflexible polymer nanocomposite melts, and the corresponding fitting data is also given in Figure 2. For the rodlike polymers with polymer bead number densities ρ0σ3 = 0.4−0.8, the polymer bead number density brings significant influences on the curves of entropic interaction. As shown in Figure 2b,
Figure 3. Reduced density fluctuation map of polymer beads ρ* around a NP in the ρ−z plane for various bending energies b = 0 (a), b = 30kBT (b), and b = 300kBT (c) with a polymer bead number density ρ0σ3 = 0.6.
with the increases of polymer number bead density the difference between our simulation data and nematic model becomes more obvious. As the polymer chains possess more configurational entropy in the higher density melts, there are stronger contact attractions for NPs in the higher density melts. In order to investigate the origin of entropic interactions in semiflexible polymer melts, we studied the effective depletion zone around a NP in various nanocomposite melts using MD 575
DOI: 10.1021/acs.jpcb.5b09551 J. Phys. Chem. B 2016, 120, 572−582
Article
The Journal of Physical Chemistry B
Figure 4. Depletion zone around a NP immersed in flexible (a), semiflexible (b), and rodlike (c) polymer melts.
Figure 5. The 3D shapes of low (upper panels, ρ* ≤ 0.8) and high (lower panels, ρ* ≥ 1.2) reduced polymer densities around a NP in polymer melts with various bending energies: left, b = 0 (a1, b1); middle, b = 30kBT (a2, b2); and right, b = 300kBT (a3, b3).
fluctuation becomes anisotropic. As shown in Figure 3b, in the lateral direction there exists an obvious oscillation. However, in the longitudinal direction the reduced density varies form ρ* = 0 to ρ* = 1.0 monotonically. Therefore, the effective depletion zone becomes anisotropic and extends along the longitudinal direction. In the case of rodlike polymer chains, as shown in Figure 3c, a more obvious anisotropy can be observed. The reduced density oscillates more severely in the lateral direction; while in the z-direction, the area with cyan marked is much larger than that for the semiflexible polymers. Therefore, the entropic interactions are anisotropic, and the optimal entropic interactions tend to be parallel to the z-axis direction. Correspondingly, snapshots of a single NP immersed in polymer nanocomposites are shown in Figure 4 for a clear look of effective depletion zone, where the depletion zone around each NP is visualized. There is a great extension in the longitudinal direction and a small narrowing in lateral direction for semiflexible polymers. In other words, an isotropic− anisotropic transition for the depletion zone near NPs in nanocomposite melts occurs, following the similar predicted
simulation. Figure 3 shows the density distribution of polymer chains around a single NP with a reduced density fluctuation map of polymer beads ρ* (ρ* = ρlocal/ρ0 with ρlocal denoting the local bead density) in the ρ−z plane for three typical bending energies b = 0, b = 30kBT, and b = 300kBT. Here, the z- axis is parallel to the nematic direction of a rodlike chain, ρ = ± x 2 + y 2 , and the sign of ρ is just used to distinguish between one-half (ρ < 0) and the other (ρ > 0). For flexible polymers, fluctuations in lateral and longitudinal directions share a similar trend. According to the local bead density fluctuation map, we can determine the effective depletion zones. The area filled with cyan corresponds to the reduced density of 0.4
