Statistical Mechanical Theory of Penetrant Diffusion in Polymer Melts

Article pubs.acs.org/Macromolecules

Statistical Mechanical Theory of Penetrant Diffusion in Polymer Melts and Glasses Rui Zhang and Kenneth S. Schweizer* Department of Materials Science and Frederick Seitz Materials Research Laboratory, University of Illinois, 1304 West Green Street, Urbana, Illinois 61801, United States S Supporting Information *

ABSTRACT: We generalize our microscopic, force-level, selfconsistent nonlinear Langevin equation theory of activated diffusion of a spherical particle in a dense hard sphere fluid to treat molecular penetrant diffusion in homopolymer melts and nonaging glasses. A coarse-grained mapping is developed where polymer chains are modeled as disconnected, noninterpenetrating Kuhn scale hard spheres (diameter, σ), and the penetrant is modeled as an effective hard sphere (diameter, d) which can be attracted to the polymer segment. The polymer mapping is a priori carried out by enforcing the effective hard sphere fluid reproduces the specific polymer liquid or glass long wavelength dimensionless collective density fluctuation amplitude. The theory predicts that penetrant diffusivity exhibits supra-Arrhenius temperature dependence in supercooled polymer melts and (near) Arrhenius temperature dependence in quenched nonequilibrium polymer glasses. Polymer−penetrant attraction slows down penetrant diffusivity to a degree that is strongly enhanced as penetrants become smaller. By treating d/σ as the only adjustable material-specific parameter, the theory is in good agreement with experimental diffusivity data spanning more than 10 decades for a wide range of penetrants (from small gas to large organic molecules), amorphous polymers, and temperatures. Optimal d/σ values are consistent with a priori physical estimations of effective spacefilling molecular and Kuhn segment diameters. Through comparative studies, two different a priori choices of penetrant−matrix attraction strength are established for small gas and large organic penetrants. System parameter transferability is examined. The theory represents a microscopic-based statistical mechanical approach for penetrant diffusion in polymers and provides a foundation for treating time-dependent penetrant diffusivity in aging polymer glasses, collective effects induced by finite penetrant loading, and diffusion in heterogeneous polymeric materials.

I. INTRODUCTION The diffusive transport of molecules in polymers is a fundamental physical phenomenon that often plays an important role in polymeric materials. For example, monomer diffusion controls the reaction rates in step-growth polymerization, migration of activated additives in the curing process is a decisive factor for the final morphology of the thermoset resin, and the dynamics of molecular plasticizers/antiplasticizers can greatly influence the material mechanical properties. Penetrant diffusion is also at the heart of many polymerbased applications such as gas diffusion in polymer membranes where selectivity and rapid transport are desired for separation applications,1 barrier materials where gas permeation in coatings is unwanted,2 the long-term sequestration and controlled release of reactive molecules in microcapsules for self-healing applications,3,4 and water purification membranes.5 Penetrant diffusion rates in homopolymer matrices vary over many orders of magnitude and are affected by a wide range of physical and chemical factors. The matrix of interest can be a normal or supercooled liquid, a rubbery network, or a glass. Penetrants can vary greatly in size and shape and their cohesive © XXXX American Chemical Society

interaction with the polymer. Because of this material-specific complexity, and since the penetrant diffusivity can vary over a very wide range, predicting from first-principles the long time activated penetrant diffusion in polymer matrices remains a challenging theoretical problem. Notably, penetrant long time diffusivity is often too small to be directly measured using molecular dynamics or Monte Carlo simulations,6 even for small penetrants such as water.7 While combining simulation and transition-state theory allows penetrant diffusivity several decades slower to be computed for small gas/vapor molecule diffusion in glassy polymers,8 this method is yet to address the even slower penetrant diffusion in barrier and other polymeric materials. Moreover, quantitative predictions for real materials based on simulation are always limited by the quality of force fields. Several phenomenological models of penetrant diffusion in polymers have been developed over the decades. Among them, Received: April 8, 2016 Revised: July 11, 2016

A

DOI: 10.1021/acs.macromol.6b00725 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules

and quantitatively apply our SCNLE approach21 to molecular penetrant diffusion in rubbery and glassy polymer matrices. The remainder of the article is organized as follows. The key elements of SCNLE theory for the primitive hard sphere penetrant system are concisely reviewed in section II.A, and the new generalization of the mapping approach to a thermal polymer matrix is described in section II.B. Section III presents model calculations for penetrant diffusion constant in poly(methyl methacrylate) (PMMA) melts and glasses and discusses general effects of the penetrant−matrix size ratio, temperature, and cohesive attraction. Section IV proposes several a priori methods to estimate the effective penetrant and matrix sphere diameters and penetrant−matrix attraction strength for specific real systems and the mathematical criterion employed for the optimal value of d/σ. Application of our theory to gas/vapor molecule diffusion in various polymers is presented in section V and to the diffusion of larger organic molecules in section VI. The article concludes in section VII with a summary and discussion. The Appendix in the Supporting Information includes additional calculations for 20 penetrant−polymer pairs that serve as predictions which can be tested in future experiments.

the free volume theory is the most widely used. The free volume concept was first proposed by Cohen and Turnbull (CT)9 to treat self-diffusion in hard sphere fluids and molecular liquids and was later generalized by Fujita10 to solvent diffusion in polymers. The modern solvent−polymer free volume theory has been systematically established by Vrentas and Duda (VD).11−14 It postulates that the volume of an amorphous polymer consists of occupied volume, interstitial free volume, and hole free volume,11,14 and only hole free volume is responsible for solvent and polymer diffusion by occasionally opening up voids large enough to allow hopping of elementary units.11,14 By building on the CT9 approach for hard sphere fluids, VD empirically constructed three generalized equations (containing 10 independent parameters) for solvent and polymer self-diffusivity and their mutual diffusion coefficient in solid-state solvent/polymer systems.11 In the penetrant infinite dilution limit considered in this article, the penetrant diffusion constant obeys11 D∞(T ) = D01(T ) exp( −ξV 2*/VFH2(T ))

(1)

Here, four material specific parameters or (thermodynamicstate-dependent) functions enter: the pre-exponential D01(T) that may depend on temperature, for example, in an Arrhenius manner if the penetrant adsorbs on the polymer, V*2 is the critical hole free volume for polymer segment motion (a steric repulsion determined quantity), ξ is the ratio between critical molar volume of solvent and polymer jumping unit, and VFH2 is the temperature-dependent polymer hole free volume. Although the VD free volume model is neither microscopic nor derived from fundamental statistical mechanics, it has been successfully employed to model many experimental observations based on empirical procedures to determine the required (“standardized”) parameters from other material properties (e.g., viscosity, thermal expansion coefficient, matrix glass transition temperature, etc.).15,16 The phenomenological nature of the present understanding of penetrant diffusion in polymers is unsurprising given the underlying physical problem involves activated motion in cold liquids or glasses. Understanding the α relaxation and glass transition of one-component molecular and polymeric liquids remains a frontier scientific area. However, much progress has been recently made for the latter problem, and we believe the time is right to build on these advances to attack the problem of penetrant diffusion at a microscopic force-level using statistical mechanics. We build on the “elastically collective nonlinear Langevin equation” (ECNLE) theory of Mirigian and Schweizer (MS) which has successfully treated the α relaxation of molecular and polymer liquids, with no adjustable parameters, over 10−14 decades in time.17−20 Penetrant diffusion is perhaps an even harder problem due to the asymmetry in structure, interactions, and relevant time scales for a probe molecule in a glassy liquid. The goal of this article is to present, to the best of our knowledge, the first general microscopic theory for penetrant diffusion in polymers in the dilute penetrant limit. It is applicable for a wide range of penetrant sizes and degrees of matrix mobility. In very recent work,21 we heuristically developed the force-level, self-consistent nonlinear Langevin equation (SCNLE) theory for activated diffusion of a single penetrant in viscous small molecule fluids and colloidal suspensions. The basic formulation is for a “primitive sphere system” and allows a unified treatment of the three (predicted) activated hopping regimes of diffusion. In this article, we extend

II. THEORY, MODELS AND MAPPING We very recently formulated the SCNLE penetrant activated diffusion theory21 for a single hard sphere penetrant (diameter, d) dissolved in a dense hard sphere (diameter, σ) fluid of volume fraction, ϕHS = πσ3ρm/6, where ρm is the matrix fluid number density. This approach allows a no-fit parameter, selfconsistent calculation of the penetrant activated barrier hopping time and long-time diffusion constant for diverse penetrant− matrix (p−m) size ratios, p−m attraction strengths, and matrix volume fractions. In section II.A, we briefly review key ingredients of this SCNLE theory; full technical details can be found in ref 21. In section II.B we present the primary new methodological advance in this article which generalizes SCNLE theory to treat molecular penetrant diffusion in homopolymer melts and nonaging glasses. II.A. Hard Sphere SCNLE Theory. Figure 1 illustrates the conceptual problem that the theory aims to address: the timedependent trajectory of a single penetrant dissolved in a dense matrix fluid where penetrant motion is activated. The goal here is to predict the mean activated hopping time and diffusion constant of the penetrant as a function of the species size asymmetry, intermolecular forces, the α relaxation process of the matrix fluid, and thermodynamic state. The starting point is an overdamped, stochastic nonlinear Langevin equation (NLE)21 for the scalar (angularly averaged) displacement of the penetrant from its initial position, rp(t) ζs , p

dr p dt

=−

dFdyn , p(rp , t ) dr p

+ δfs , p (t )

(2)

where rp(t = 0) = 0 and ζs,p is the short time penetrant friction constant calculated using binary collision mean field (BCMF) theory.21−25 The central physical quantity controlling hopping is the time-dependent penetrant “dynamic free energy”21 B

DOI: 10.1021/acs.macromol.6b00725 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules

Physically, the self-consistency condition determines the extent to which the matrix undergoes an α process in order for the penetrant to hop.21 For small enough penetrants, hopping can occur when the matrix only undergoes vibrational motion, per an amorphous solid. For penentrants and matrix particles comparable in size, the matrix must undergo an α relaxation event to allow “enough room” for a penetrant hop. Intermediate sized penetrants require an intermediate amount of matrix relaxation, with a significant time scale separation between penetrant hopping and matrix α relaxation still existing. Given the matrix packing fraction, penetrant−matrix diameter ratio, d/σ, penetrant−matrix pair potential, upm(r), and βK externally chosen,21 the penetrant mean barrier hopping time can be self-consistently determined. The penetrant displacement associated with the localization well and barrier location in the dynamic free energy define rloc,p and rB,p, respectively. Finally, the long time Fickian penetrant hopping diffusion constant is estimated as

Dhop , p = Δrp 2/6τhop , p

where the penetrant jump distance is Δrp = rB,p − rloc,p. II.B. Extension to Molecular Penetrants in a Polymer Matrix. The theory of section II.A can treat real material systems upon introducing a coarse-graining mapping scheme first proposed for one-component thermal liquids composed of nonspherical rigid molecules.17,19 Thermal liquids are replaced by an effective hard sphere fluid based on requiring their dimensionless compressibilities, Sexpt = ρmkBTκT, are identi0 cal.17,19 This idea was successfully implemented in our prior SCNLE study of penetrant diffusion in small molecule liquids, and good agreement with experiment was demonstrated.21 Typical van der Waals and moderately polar molecular liquids have an experimental Sexpt that is well represented as19 0

Figure 1. Schematic of the coarse-graining scheme that maps a real penetrant−polymer system to a model sphere mixture (see the text for details).

Fdyn , p(rp , t ) kBT

1 2π 2

= −3 ln(rp/d) −

∫0

∞

dk

2 2

× k 2[ρm Cpm 2(k)Smm(k)][e−k rp /6] 2

βK

2

× [e−k rloc ,m /6Smm(k)e−(t / τm(k)) ]

(3)

It quantifies the effective confining field (its gradient determines the force) exerted on the penetrant by the dynamically evolving matrix. Here, kB is the Boltzmann constant, T is temperature, Cpm(k) is the wavevector-dependent p−m direct correlation function, and Smm(k) is the matrix collective static structure factor. All structural inputs are calculated from standard integral equation theory based on the Percus−Yevick (PY) closure relations.21,24 The initial dynamical constraints (effective forces) enter via structural pair correlations. These constraints can relax via matrix structural relaxation which enters through the wavevector-dependent matrix collective relaxation time, τm(k), as precisely constructed in refs 21, 26, and 27; its time scale is set by τα,m, the mean α relaxation time of the pure matrix fluid computed using ECNLE theory.18 A distribution of matrix relaxation times enters via the KWW stretching exponent, βK. The transient matrix localization length from ECNLE theory,17−20 rloc,m, quantifies the matrix Debye−Waller factor. In the limit the matrix is a glass or amorphous solid, τm(k) → ∞ in eq 3. The dynamic free energy21 is employed to formally compute a “time-dependent” penetrant hopping rate using Kramers theory.28 To render the approach temporally self-consistent, matrix relaxation is evaluated at a time equal to the (unknown) penetrant relaxation time. This requires the penetrant mean barrier hopping time τhop,p to obey the self-consistent equation:21 τhop , p =

2τs , p d2

∫r

rB, p

dr e Fdyn,p(r , τhop,p)/ kBT

loc , p

× e−Fdyn,p(r ′ , τhop,p)/ kBT

∫r

(5)

1 S0expt

= −A Ns +

B Ns T

(6)

where A and B are site-level entropic (or packing) and liquid cohesion contributions, respectively, for a molecule composed of Ns “interaction sites” (e.g., a site in benzene is a CH group, and hence Ns = 6). For hard spheres, PY integral equation 4 2 theory predicts24 SHS 0 = (1 − ϕHS) /(1 + 2ϕHS) . Equating this expt to S0 obtained from experimental equation-of-state data, the system and temperature-dependent effective packing fraction is ϕHS = 1 + f −

f 2 + 3f

−1 ⎛ B Ns ⎞ ⎟⎟ f = ⎜⎜ −A Ns + T ⎠ ⎝

(7)

The specific chemistry enters via A, B, and Ns. Very recently, the above mapping ideas have been extended to pure polymer melts by disconnecting chains at the level of “Kuhn-sized” hard spheres.20 ECNLE theory captures quite well the polymer α relaxation time over 10−14 decades, Tg values, and emergent glassy elasticity. Despite the global “disconnection” approximation, the influence of chain length still enters in a system-specific manner via the Kuhn length (and hence Ns above) since backbone stiffness varies with molecular weight via the characteristic ratio; quantitative calculations yield excellent predictions for how Tg varies with

r

dr ′

loc , p

(4) C

DOI: 10.1021/acs.macromol.6b00725 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules

Table 1. Key Parameters (Symbols Defined in the Text) for the Five Amorphous Polymers Studied: Poly(vinyl acetate) (PVAC), Poly(methyl methacrylate) (PMMA), Polystyrene (PS), Polycarbonate (PC), and Polysulfone (PSF) (References Refer to Sources of the Numbers) polymer PVAc PMMA PS PC PSF a

Ns 18.4 23.5 26.2 36.4 38.0

Aa 0.762 0.898 0.62 0.674 0.558

Ba (K) 1285 1460 1297 1280 1343

Tg (K) 32

308 37332 37335 41435 45936

αg (10−4/K)b 2.8 2.7 2.5 2.6 2.2

M (Da) 264.1 336.1 341.0 487.7 525.4

σest (Å) 8.9 9.8 10.2 11.0 11.2

σcali (Å) 10.8 11.8 13.5 13.8 13.8c

lk (Å)

ρb (g/cm3)

33

13.6 13.834 14.734 n/a n/a

1.186 1.149 1.023 1.16 1.194

Data from ref 20. bData for polymer glass near Tg from ref 31. cSet as the calibrated value from our O2/PC study.

chain length.20 For molecular penetrants smaller than or comparable to the Kuhn size scale, we physically expect that only local polymer segmental motion is relevant in determining penetrant hopping, and longer range consequences of chain connectivity are minimal. With the above motivation, the penetrant/polymer system is mapped to an effective hard sphere mixture as follows (see Figure 1). For equilibrium polymer melts (T > Tg, where τα,m(Tg) = 100 s), eqs 6 and 7 still hold. The A and B values for five amorphous polymers studied in this article are listed in Table 1. We use a revised scheme to determine Ns compared to the prior pure polymer liquid work20 based on a literal coarse graining scale of the a priori determined Kuhn length. Here, solely to improve quantitative accuracy for the penetrant diffusion problem of interest, we adopt the experimental Tg of the pure polymer matrix as input to calibrate Ns by enforcing ECNLE theory reproduces the experimental Tg. The thusobtained Ns values are listed in Table 1. Importantly, these Ns values are indeed very close to the a priori ones determined directly from polymer conformational statistics as employed previously (Table 1 of ref 20). In this article we do not address the influence of polymer molecular weight on penetrant diffusivity. However, this issue can be studied based on prior work.20,21 Specifically, at fixed chemistry, as polymer chain length grows the typical behavior is Tg and the Kuhn length (and Ns) increase, until saturating in the long chain limit. This implies the ratio of penetrant to Kuhn segment size decreases, and the dimensionless compressibility decreases and hence caging constraints increase.20 These finite chain length effects are present in the dynamical theory and will modify penetrant diffusion. We also study here, for the first time, penetrant diffusion in glassy polymers (T < Tg). For simplicity, we focus on “freshly quenched” or nonaging nonequilibrium polymer glasses, which requies a different mapping scheme. We argue that since below Tg the time scale of polymer structural relaxation exceeds the experimental scale, the corresponding thermal density fluctuations, and hence effective hard sphere volume fraction, only depends on the phonon-like vibration of polymer segments. In this limit, the matrix is effectively “isostructural”, evolving with cooling only via thermal contraction. We thus make the very simple, but no adjustable parameter, assumption that only the density of the polymer glass changes (increases) with cooling via the thermal expansion effect. Employing the experimental polymer thermal expansion coefficient (αg at Tg, Table 1), we obtain for T < Tg ϕHS = ϕg [1 + (Tg − T )αg ]

upm(r ) = −εe−(r − dpm)/ a ,

r ≥ dpm = (d + σ )/2

(9)

where ε and a define the attraction strength and range, respectively. We fix a = 0.4dpm, a reasonable a priori value for nonpolar van der Waals (vdW) molecules.21 The above theoretical ideas allow prediction of the penetrant diffusion constant given d/σ, ε and βK(T). We find that the influence of βK(T) on our diffusivity calculations is minor.21 Given this, and the fact that quantitative experimental knowledge of βK(T) for specific polymers is often limited, we fix βK ≡ 1/3 for all polymers. For poly(vinyl acetate) (PVAc), we sometimes employ βK(T) = 0.51 + 0.0021(T − Tg), as experimentally deduced.29 Concerning d and σ, to be physically consistent their “best values” (if defined in a practical sense that the resulting SCNLE calculation most closely matches experimental data) should be correlated with chemically sensible a priori estimates of the penetrant and polymer Kuhn space-filling diameters. Obviously, a unique first-principles calculation of d and σ is not possible because of the “sphericalization approximation” adopted for real molecules and Kuhn scale segments. But what matters most in the theory is the size ratio, d/σ, which is a steric metric of how a (generally) nonspherical penentrant molecule locally fits into the polymer liquid. Our theory predicts penetrant diffusivity is exponentially sensitive to d/σ for smaller penetrants.21 Since the polymer is mapped to a hard sphere matrix, a “reasonable” value of ε in eq 9 for the effective sphere mixture system should be compared to the εpm value estimated from chemical knowledge of penetrant−segment interactions. We return to these parameter selection issues for real materials in sections IV, V and VI and now present model calculations for the generic influence of d/σ, ε, and Tg/T on penetrant diffusion.

III. MODEL SYSTEM CALCULATIONS All calculations in this section adopt PMMA (Tg = 373 K) as the polymer matrix (A = 0.898, B = 1460 K, Ns = 23.5 in eq 6). All characteristic trends shown are qualitatively the same across different amorphous polymers. Figure 2 presents the penetrant diffusivity for six representative d/σ values (0.3−2) as a function of Tg/T. The p−m attraction strengths (in units of kB) are 0, 500, 900, and 1300 K. These parameters cover essentially the entire range of common material chemistry, from weak to strong p−m attractions and from small to large penetrant sizes. A universal trend is that above and below the polymer Tg the penetrant diffusivity exhibits a supra-Arrhenius and nearly Arrhenius temperature dependence, respectively. An Arrhenius behavior in glasses is a common feature in nonaging liquids and polymers below Tg. To the best of our knowledge, experiments also generally find near-Arrhenius behavior of penetrant

(8)

Finally, nonspherical molecular penetrants are mapped to an effective hard sphere of diameter, d, which interact with the matrix via a hard core plus exponential attraction: D

DOI: 10.1021/acs.macromol.6b00725 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules

range of temperature when d/σ ≥ 1 for weakly attractive penetrants.21 This leads to a modest distinction between regime II (0.5 ≤ d/σ ≤ 1) and regime III (d/σ ≥ 1) for small ε, although such a distinction becomes “smeared” for strongly attractive penetrants.21 The general physical behavior in regimes II and III is that penetrant diffusion requires either a modest or alot, respectively, facilitation via the matrix α relaxation; for sufficiently large d/σ or ε, penetrant motion is essentially slaved to the matrix α relaxation process. Figure 2 shows that Dhop,p curves as a function of Tg/T for large d/σ or ε (e.g., d/σ = 0.5, 1, and 2 at ε = 1300 and d/σ = 2 at ε = 0) are close and parallel to each other. Third, as deducible from the second point above and seen in Figure 2, the penetrant diffusivity strongly decreases with increasing ε when d/σ is small, but the ε dependence is weak for larger penetrants. As an application of the results in Figure 2, Figure 3 presents calculations of the average time (τdiff) for a penetrant to diffuse

Figure 2. Penetrant hopping diffusion constant (in absolute units) in a PMMA matrix as a function of Tg/T for six representative p−m size ratios: (a) 0.3, 0.35 and 0.4; (b) 0.5, 1, and 2. Solid, dashed, dotted, and dotted-dashed curves (if shown) represent results for ε/kB (in Kelvin) = 0, 500, 900, and 1300, respectively. Dhop,p calculations here and in all other figures employ a KWW factor of βK ≡ 1/3, unless otherwise noted.

Figure 3. Average penetrant diffusion time for moving a distance Δl = 100 nm calculated from Ficks law τdiff = Δl2/6Dhop,p as a function of Tg/T for five representative p−m size ratios in PMMA: 0.3, 0.4, 0.5, 1, and 2. Two different p−m attraction strengths are considered: ε/kB (in Kelvin) = 0 (solid) and 1300 (dotted-dashed); the latter represents the case of very strong p−m attraction. Horizontal dotted lines indicate positions of six distinctive time scales (see labels).

diffusivity in a nonaging glassy matrix, although most data are focused on small gas/vapor penetrants because their diffusivities are still large enough to be measured.14,16 The theory can treat arbitrarily slow penetrant diffusion, although most experiments cannot detect a diffusivity below ∼10−18 cm2/s. We note that the physical origin of the Arrhenius behavior in quenched nonaging glasses arises from (by assumption) thermal contraction (densification) of the polymer matrix with cooling. The general influences of d/σ and ε were established in detail in ref 21. Here, we briefly review for context the essentials based on Figure 4 in ref 21 and Figure 2 in this article. First, the penetrant diffusivity decreases the fastest with d/σ in the smallpenetrant “regime I” (d/σ < ∼0.5), where the jump distance universally increases with d/σ for all ε. Comparative calculations show the penetrant diffusion constants calculated using the full SCNLE theory are consistent with a limiting calculation in which the matrix is treated as a vibrating amorphous solid. Physically, this means the activated hopping time scale of small penetrants is much faster than the matrix α relaxation time. A direct consequence is that the matrix KWW exponent βK does not affect the penetrant diffusivity results in regime I. Second, the decrease of penetrant diffusivity with d/σ dramatically slows down as the penetrant becomes larger than about one-half of the matrix sphere size. The jump distance increases with d/σ starting at d/σ ∼ 0.5 (attaining a maximum ∼σ at T = Tg) and is roughly equal to ∼0.5−0.6σ over a wide

a distance of 100 nm, a length scale representative of the thickness of polymer coatings,2 gas/liquid separation membranes,1,5 and polymeric shells of microcapsules used in selfhealing materials.3,4 For any other choice of penetrant diffusion distance Δl, the diffusion time simply scales (by assumption) as τdiff(Δl/100 nm)2. One sees that for the diverse range of p−m size ratios, p−m attraction strengths, and scaled temperatures included in Figure 3, τdiff varies from microseconds to centuries. A finding that may be relevant for using polymers as barrier materials is that temperature must be well controlled since no matter how large and/or attractive the penetrant is, τdiff quickly drops to around seconds if temperature is increased to only ∼1.1Tg. “High Tg” is thus a critical design principle. The influence of pentrant−matrix attraction is presented in a more quantitative fashion in Figure 4, which shows penetrant hopping diffusion constants at T = Tg as a function of ε for the six representative p−m size ratios discussed in Figure 2. For the three smallest penetrants considered (d/σ = 0.3, 0.35, and 0.4), the ε dependence is initially moderate as Dhop,p decreases by only ∼2−3 orders of magnitude when ε/kB grows from 0 to 800 K (∼2.1kBT for PMMA). A much stronger influence of matrix−penetrant attraction emerges as ε/kB increases beyond E

DOI: 10.1021/acs.macromol.6b00725 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules

weak. On the other hand, for small d/σ values such as 0.3, since the penetrant diffusion time scale is still orders of magnitude faster than matrix relaxation at very high ε, the larger local constraining force with increasing ε directly results in a very large slowing down of penetrant diffusion. Before addressing specific real materials, we introduce two physical quantities of interest. The first we call “penetrant diffusion fragility” (PDF): m=−

d(log(Dhop , p(Tg /T )) d(Tg /T )

T = Tg

(10)

This quantity measures the rate of change of penetrant diffusivity with temperature. For the pure PMMA liquid, the dynamic fragility is ∼85 based on SCNLE theory.21 The closer the PDF is to this value, the stronger the “slaving” of penetrant motion is to matrix relaxation. As Figure 6 shows, the PDF values for small, weakly attractive penetrants (representative of common gases) are well

Figure 4. Penetrant hopping diffusion constant as a function of p−m attraction strength in PMMA at T = Tg (373 K) for six representative p−m size ratios: 0.3, 0.35, 0.4, 0.5, 1 and 2.

∼800 K. From ε/kB = 800 to 1500 K (∼2−4kBT for PMMA), Dhop,p decreases by ∼5−7 orders of magnitude for the three small model penetrants. For larger d/σ, the effect of ε weakens signifcantly, and for d/σ = 2 the reduction of Dhop,p is less than a decade when ε/kB is increased from 0 to 1500 K. As a quantitative measure of the local tendency of penetrants to “adsorb” on Kuhn scale segments, Figure 5 presents the

Figure 6. Penetrant diffusion fragility as a function of p−m attraction strength at T = Tg (373 K) in PMMA for six representative p−m size ratios: 0.3, 0.35, 0.4, 0.5, 1, and 2.

below 20, indicating strong decoupling between penetrant motion and matrix relaxation. The PDF universally approaches 85 for large enough d/σ or ε, signifying the “slaved” limit. We also note that the shape of each curve in Figure 6 is similar to its inverted analog in Figure 4, indicating that the magnitudes of PDF and Dhop,p(T = Tg) are well correlated. A second dynamical quantity of experimental interest in the nonaged glassy regime is the penetrant activation energy, or effective barrier, Eact,30 which is extracted by fitting our theoretical results at T ≤ T g to the expression Dhop,p = D0 exp( − Eact /kBT ); we find very high correlation coefficients of >0.999 in all cases. Figure 7 shows the extracted Eact (scaled by kBTg) as a function of ε. One sees that Eact strongly increases with ε for small penetrants in regime I (7fold from ε/kB = 0 to 1500 K), while for larger d/σ the influence of ε is weak. An interesting result is that for sufficiently strong p−m attraction the theory predicts curve crossing and hence a moderately larger Eact for smaller d/σ. This seemingly surprising finding does not conflict with the physically sensible trend (see Figure 2) that at fixed ε and temperature, larger penetrants diffuse slower than smaller penetrants. Rather, it is a consequence of the fact that at sufficiently high p−m attraction, the penetrant to matrix

Figure 5. Contact value of the penetrant−matrix pair correlation function as a function of p−m attraction strength for d/σ = (from bottom to top) 0.3, 0.5, and 1. For each p−m size ratio, the solid and dashed curves represent results at Tg/T = 0.75 and 1.0, respectively.

contact value of the p−m pair correlation functions for selective systems that exhibit the general influence of d/σ, ε, and Tg/T ≤ 1. Since calculations are performed at a fixed value of the latter parameter, the x-axis in the plot corresponds to varying the ratio ε/kBTg. The results clearly demonstrate that the insensitivity of penetrant diffusivity to ε for large d/σ is not of structural origin since the p−m contact values increase strongly (with nearly the same slope) with ε. We interpret this insensitivity as a physical consequence of the following effect. Even if only the steric force contribution is considered (ε = 0), the confining force on a large penetrant is very strong, and its diffusion rate is controlled by the matrix relaxation (“slaved” penetrant diffusion). Thus, although the local force on the penetrant increases with ε, because the matrix structure (which controls the matrix relaxation rate within our theory) is independent of p−m interactions in the dilute penetrant limit, the dynamical effect of ε on large penetrant diffusivity is very F

DOI: 10.1021/acs.macromol.6b00725 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules

careful discussion. Recall the m−m attraction is embedded in an average sense via mapping the polymer matrix to an effective hard sphere fluid that reproduces the experimental S0(T).20 Since the matrix has been mapped to hard spheres, ε should in principle be adjusted accordingly. For example, intuitively if the penetrant chemistry is similar to that of the polymer, and penetrant and segment sizes are close, ε = 0 is a plausible zeroth-order choice as there is no preference for a segment to attract another segment or the penetrant as quantified in gpm(r). On the other hand, for small gas/vapor molecules, a penetrant does not interact with a polymer segment as a literal smooth sphere, but rather will prefer specific location(s) determined by the Angstrom-scale chemical heterogeneity of a polymer monomer and steric shape considerations. Moreover, if d is small, a penetrant will stick to just one or a few matrix spheres because locally the amount of open space is relatively large. Thus, we argue that the relevant zeroth-order choice of ε should be the real (bare gas phase) p−m attraction strength (εpm) for small gas/vapor penetrants. Since do not wish to fit ε and d/σ independently, we focus on two limiting cases: ε = 0 and ε = εpm. We will show in sections V and VI that for each system one of the two choices works well. For the cross-energy parameter, the geometric mean estimate εpm ≈ (εpεm)1/2 is adopted, where εp for gas molecules are taken as their Lennard-Jones (L-J) parameters from ref 30 (Table 2).

Figure 7. Effective penetrant activation energy Eact (scaled by kBTg) obtained by fitting our penetrant hopping diffusivity data in PMMA at T ≤ Tg to an Arrhenius form Dhop,p = D0 exp(−Eact/kBT) plotted as a function of p−m attraction strength for six representative p−m size ratios: 0.3, 0.35, 0.4, 0.5, 1, and 2.

hopping time ratio becomes smaller with cooling due to tighter packing and penetrant dynamics becomes more “slaved” with matrix.

IV. ESTIMATION AND OPTIMIZATION OF PARAMETERS FOR REAL MATERIALS We now consider real world penetrant−polymer systems. A key to map the penetrant diffusion constant to a location in Figure 2 is the identification of the “correct” d/σ and ε. It is first instructive to recall how chemical complexity (molecular shape and attractive interactions) is handled in a predictive way in ECNLE theory for one-component liquids.17−20 Here the molecular size, d, scales out of the problem based on the hard sphere mapping since only the dimensionless packing fraction enters. The effects of molecular shape and cohesive attractions are captured in an average, but no adjustable parameter, manner by requiring the S0(T) deduced from equation-of-state data is exactly reproduced by the mapped hard sphere reference fluid. The above mapping approach is again used for the matrix liquid, but additional information is required to determine the penetrant effective diameter and its interaction with the matrix, specifically d/σ and ε in eq 9. The ratio d/σ reflects the matrix segment−penetrant molecule pair, not two independent parameters. To address this issue, we are guided by the desire to minimize the number of fit parameters to the smallest we believe is possible, one. We first introduce physically motivated a priori approaches to directly estimate the required parameters. Consistent with the “effective sphere” picture proposed by MS,20 the diameter σ of the effective hard sphere fluids can be estimated as the spacefilling diameter of the polymer “Kuhn sphere”:20 σest = (6M /πρmatrix )1/3

Table 2. Key Parameters for Gas and Vapor Penetrants Studied (Symbols Defined in the Text) penetrant

m (Da)

εp/kBa (K)

dKa (Å)

dLJa (Å)

dBHb (Å)

O2 N2 CO2 CH4 C3H8 H2O CH3OH

32 28 44 16 44 18 32

106.7 71.4 195.2 148.6 237.1 809.1 234c

3.46 3.64 3.3 3.8 4.3 2.65 n/a

3.467 3.798 3.941 3.758 5.118 2.641 3.69c

3.38 3.63 3.94 3.71 5.15 2.77 3.71

a Data from ref 30 (except for methanol). bBarker−Henderson diameters calculated using eq 14 at room temperature (T = 300 K). c Data from ref 37.

For matrix segments and large aromatic penetrants, we estimate εm/kB ≈ 90Ns K, where the number 90 is “calibrated” from the benzene L-J parameter (εm/kB ≈ 540 K,37 six rigid sites), which has a chemistry similar to the three aromatic penetrants and five polymers investigated in this article. The reference hard sphere penetrant diameter d should reasonably represent the steric size of a non-spherical molecular penetrant. For large aromatic penetrants, due to their similar size and chemistry compared to the polymer Kuhn segments considered, density calibration is a plausible option: ⎛ ⎞1/3 6m ⎟ ⎜ dest = ⎜ ⎟ ⎝ πρp , RT ⎠

(11)

(12)

where m is penetrant mass and ρp,RT is the penetrant fluid density at the room temperature (a choice made because of easily accessible data). Diameters of three aromatic penetrants estimated this way are listed in Table 3. For small gas/vapor penetrants, one option is based on penetrant L-J parameters and the Barker−Henderson (BH) statistical mechanics method:24

Here, M is the mass of the effective sphere (equals monomer mass times number of monomers per segment), and ρmatrix is polymer mass density. The matrix sphere diameter follows by choosing ρmatrix = ρg (polymer density at Tg); relevant parameters are summarized in Table 1. How to determine the appropriate value of the p−m attraction strength ε in the mapped sphere system requires G

DOI: 10.1021/acs.macromol.6b00725 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules Table 3. Key Parameters for Aromatic Penetrants Studied (Symbols Defined in the Text) penetrant

m (Da)

dest (Å)

ρa (g/cm3)

toluene tetracene rubrene

92.1 228.3 532.7

6.9 8.1 11.3

0.865 1.35 1.18

a

Density of pure penetrant liquids near room temperature from www. wikipedia.org. WCA(r )

⎡⎛ σLJ ⎞12 ⎛ σLJ ⎞6 1⎤ = 4εLJ ⎢⎜ ⎟ − ⎜ ⎟ + ⎥ , r ≤ 21/6σLJ ⎝ r ⎠ 4⎦ ⎣⎝ r ⎠ = 0, r > 21/6σLJ

dBH =

∫0

∞

(13)

dr(1 − e−uWCA(r)/ kBT )

Figure 8. Kuhn scale segment hopping diffusion constant calculated using one-component SCNLE theory21 as a function of Tg/T for the five polymers studied in this article. Dashed curve shows results for PVAc based on βK(T) = 0.51 + 0.0021(T − Tg).29

(14)

where σLJ and εLJ are L-J diameter and attraction strength, respectively. We have also examined another a priori choice of gas/vapor sphere diameters based on “kinetic diameters” dK summarized by Breck.30,38 “Kinetic diameters” are typically thought of as collision diameters, or distances of closest approach, of a pair of small molecules or atoms modeled at the simple L-J level with pair potential parameters chosen by fitting the second virial coefficient or shear viscosity data. Per ref 30, dK and σLJ are the two most widely used estimates of penetrant size for gas diffusion studies. All relevant gas/vapor parameters employed in this article are listed in Table 2. We will show in section V that calculations which employ dK consistently yield good agreement with experiment. In section V (gas/vapor penetrants) and section VI (large aromatic penetrants), we treat d/σ as an adjustable parameter and a priori set all other input parameters. The mathematical criterion for the best value of d/σ is as follows. Let the experimental penetrant diffusivity data (in units of cm2/s) at different temperatures in a specific penetrant−polymer system be {D̃ p,exp(T1), D̃ p,exp(T2), ..., D̃ p,exp(Tn)}, where Ti’s represent n temperatures at which experimental data are available. Let the corresponding n theoretically computed diffusion constants obtained using a specific value of d/σ be denoted as {D̃ p,th(T1,d/σ), D̃ p,th(T2,d/σ), ..., D̃ p,th(Tn,d/σ)}. Then define a so-called theory−experiment deviation (TED) number λ as follows:

constants in these polymers in Figures 9−13. We note that the five curves (all use βK ≡ 1/3) in Figure 8 almost collapse in the rubbery regime since ECNLE theory predicts little α relaxation time fragility change for these polymers.20 The two curves for PVAc employ βK ≡ 1/3 and βK(T) = 0.51 + 0.0021(T − Tg)29 and do slightly differ, indicating the minor effect of the KWW factor. In the glassy regime, curves for different polymers splay apart somewhat as a result of their different polymer thermal expansion coefficients (listed in Table 1).

V. DIFFUSION OF SMALL GAS/VAPOR PENETRANTS Our prior work considered atomic and molecular penetrant diffusion in small molecule liquids. Specifically, we studied the rare gas tracers argon, krypton, and xenon in supercooled methanol and toluene21 and aromatic tracers (tetracene and rubrene) in supercooled o-terphenyl.21 We found these experimental systems rarely fall into regime I because the size ratio between penetrant and matrix molecule is not small enough. SCNLE penetrant diffusion theory thus has not really been tested against material systems in regime I. For penetrant diffusion in polymers, however, many gas/vapor penetrant− polymer systems are in regime I. This section presents our results for these systems. V.A. Temperature-Dependent Oxygen Diffusivity. To best test our theory, a systematic set of diffusivity data for a fixed gaseous penetrant in multiple amorphous polymers over a wide range of Tg/T is desired. On the basis of a literature search, we find the oxygen penetrant best serves this purpose. Figure 9 presents the optimized theory (only one parameter, d/σ, is adjusted) and available experimental data of oxygen diffusivity in four chemically distinct polymers (PVAc, PMMA, PS, and PC) as a function of Tg/T. As we a priori expect, the optimized theory results based on ε = εpm (solid curves) generally agree better with experiments than those based on ε = 0 (dashed curves). Table 4 lists the optimal d/σ values and associated devth‑exp numbers. In PMMA, PS, and PC, the optimal values of devth‑exp for ε = 0 are 1.7−2.9 times as large as the corresponding values for ε = εpm. Only in PVAc are the two devth‑exp values close. We believe that the theory−experiment comparison of oxygen diffusivity in PS and PC matrices is more significant than in PVAc and PMMA since the former examines a wider temperature window and hence a larger variation of oxygen diffusion constants (∼1.5 orders of magnitude). The

n

λ (d / σ ) =

∑i = 1 {log[Dp̃ ,th (Ti , d /σ )] − log[Dp̃ ,exp(Ti )]}2 n (15)

The optimal d/σ value is taken to be the one that yields the smallest λ based on SCNLE theory calculations. To quantitatively compare the extent of agreement between theory and experiment, a more reasonable measure is devth‑exp = λoptimal/Q, where Q is the number of decades of experimental data spanned from the highest to lowest penetrant diffusivity: log[D̃ p,exp(max)/D̃ p,exp(min)]. We will show that the optimal d/σ value shows good physical consistency and transferability with the a priori estimation. As foundational background, we first apply the theory to calculate the segmental hopping diffusion constant of the five pure polymer liquids studied in this article: PVAc, PMMA, polystyrene (PS), polycarbonate (PC), and polysulfone (PSF). The results are shown in Figure 8, which serves as a reference when we present in the next section penetrant diffusion H

DOI: 10.1021/acs.macromol.6b00725 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules Table 4. Summary of Key Parameters and Results for Systems Discussed in Figures 9−13 penetrant/ matrix

decades of exp dataa

a priori d/σb

ε/kB (K)

optimal d/σ

devth‑exp (%)

0.350 0.322 0.330 0.292 0.294 0.254 0.304 0.252 0.334 0.208 0.412

7.4 7.1 17.9 10.1 17.5 10.5 18.8 6.6 9.1 4.6 25.6

O2/PVAc

0.7

0.389

O2/PMMA

0.5

0.353

O2/PS

1.5

0.339

O2/PC

1.7

0.315

H2O/PVAc

1.7

0.298 (0.245)

methanol/ PVAc

2.3

0.417 (0.343)

0 420 0 475 0 502 0 592 0 1157 0

0.775 (0.639)

623 0c 623c 0

0.334 0.392c 0.330c 0.982

14.6 12.0c 7.7c 6.7

0.388 0.798c 0.350c 0.434 0.282 0.610

10.2 6.9c 8.0c 8.7 6.4 2.8

toluene/ PVAc

4.1

toluene/PS

4.4

0.676 (0.511)

tetracene/ PS

6.9

0.794 (0.600)

976 0c 976c 0 1164 0

rubrene/PS tetracene/ PSF rubrene/ PSF

6.6 2.8

1.108 (0.837) 0.723 (0.587)

1954 0 0

0.232 1.102 0.464

3.0 2.9 2.2

2.9

1.009 (0.819)

0

0.628

2.0

Figure 9. Theory (curves) and experimental (square: PVAc;32 triangle: PMMA;32 circle: PS;35 diamond: PC35) results for the oxygen diffusion constant as a function of Tg/T in four different polymers. Solid and dashed curves represent the best fit d/σ (Table 4) results for ε = a priori real p−m attraction strength (εpm) and ε = 0, respectively.

set τm(k) → ∞ in eq 3) for oxygen and several other small gas/ vapor penetrants studied below in this section (N2, CO2, CH4, and H2O). Thus, there is rigorously no influence of the matrix stretching exponent, βK(T), on our predicted diffusion constants for these penetrants. The full theory results are moderately larger than frozen-matrix limiting calculations for propane and methanol penetrants, indicating they are likely just barely in regime II. V.B. Water vs Oxygen Diffusion in PVAc. Gas molecules are typically dilute in polymers and can thus be reasonably well studied using our single penetrant theory. In contrast, applying the theory to water diffusion in polymers requires a special technique to keep water loading very dilute. Experimental data of this kind are limited, in particular those addressing the temperature dependence. We are only aware of water penetrant diffusion experimental data in PVAc16 among the five polymer matrices studied in this article. In Figure 10, we present our optimized theory results and extended predictions. The oxygen diffusivity results in PVAc employing ε = εpm are included for a comparison purpose. We again perform both ε = εpm and ε = 0 calculations and search for the optimal d/σ. The two optimal d/σ values and

Calculated as log[D̃ p,exp(max)/D̃ p,exp(min)]. bCalculated as the ratio between the a priori estimation of d(dK for oxygen and water, dBH for methanol, and dest for aromatic penetrants; see main text) and the a priori estimation of σ (set as the space-filling diameter of the “Kuhn sphere”, σest). Numbers in parentheses represent the “transferred size ratio” if one sets σ = σcali listed in Table 1. cResults based on using βK(T) = 0.51 + 0.0021(T − Tg), as proposed in ref 29 for PVAc. a

devth‑exp values of oxygen diffusivity in PC are 6.6% and 18.8% for ε = εpm and ε = 0, respectively. A nonattractive penetrant approximation is thus less accurate than the ε = εpm based calculation, a chemically sensible conclusion. On the basis of optimized d/σ values, Figure 9 shows theoretical predictions for a wide range of temperatures beyond what was measured experimentally that covers the rubbery to deep glassy regimes (Tg/T = 0.8−1.6). The difference between results based on the ε = εpm and ε = 0 choices is minor in the rubbery regime. The most noticeable contrast between the two approaches lies in their predictions for the rate of penetrant diffusivity decrease in the glassy regime. One can see that the results based on ε = 0 are well removed from the experimental trend, while the calculations based on ε = εpm capture the data better. For oxygen in PS and PC, the ε = εpm calculations predict Eact = 20.5 and 24.9 kJ/mol, respectively, while the corresponding experimental values are 29.7 and 40.2 kJ/mol.35 Oxygen diffusion in polymers follows regime I physics, as evidenced by the huge separation between oxygen and polymer segmental diffusivity scales (Figure 8). Indeed, we have confirmed that full SCNLE theory results are the same as computed from the amorphous solid matrix limit (effectively

Figure 10. Best fit d/σ (see legend, Table 4) theory results (curves) for water (solid orange: ε = εpm; dashed orange: ε = 0) and oxygen (black) diffusion constants in PVAc as a function of Tg/T. Symbols represent the experimental data (circle: water;16 square: oxygen32). I

DOI: 10.1021/acs.macromol.6b00725 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules associated devth‑exp numbers are listed in Table 4. For the ε = εpm calculations, devth‑exp is 4.6% (optimal d/σ = 0.208), smaller than the analogue (9.1%) based on ε = 0 calculations (optimal d/σ = 0.334). We believe if additional experimental data in glassy PVAc become available, the ε = 0 based calculations will be increasingly less accurate than their ε = εpm analogues, per what we found in the investigation of oxygen diffusion (Figure 9 and Table 4). The solid orange curve in Figure 10 represents our prediction of water penetrant diffusion in PVAc over a wide range of temperature. We hope new experiments/simulations for this system at Tg/T > 1 will be performed in the near future to test our calculations which cover almost 6 decades of variation of the water diffusivity. Concerning the chemical sensibility of the water size derived from the optimal d/σ values, we note that the ε = εpm and ε = 0 treatments lead to effective hard-core diameters of dH2O ≈ 2.3 and 3.6 Å, respectively; here, σPVAc is set to be 10.8 Å, a sensible value calibrated by oxygen diffusion data, per section V.C. It seems reasonable to us that the real water molecule mean size lies in between these two values; for example, its van der Waals diameter is quoted as 2.82 Å.39 Since PVAc is a hydroscopic polymer and the attraction between water and PVAc is strong (εpm/kBTg,PVAc ∼ 4; for comparison, εpm/kBTg,PVAc ∼ 1.3 for oxygen penetrant), water would rather strongly adsorb onto PVAc. It is thus chemically sensible to expect the effective hardcore diameter of water in the water−PVAc system is smaller, reflecting the tight packing between a water molecule and PVAc segment. Figure 10 also illustrates an interesting competition between the penetrant size and attraction effect. Oxygen is larger than water and therefore for this reason should diffuse slower. On the other hand, because water has a stronger attraction with PVAc than oxygen, attractive interactions favor slower diffusion of water. When temperature is lower, the attraction effect is more dominant since ε/kBT becomes larger but d/σ does not change. The above analysis qualitatively explains our prediction that water diffusivity drops below oxygen diffusivity in the deep glassy regime of PVAc. V.C. Physical Consistency and Transferability. To test how well the optimal d/σ values deduced from the oxygen diffusion study can be transferred to other gas molecules, and which of dK and dBH better serves as an a priori estimation of d, we carry out comparative calculations (Tables 5 and 6). In this subsection, only ε = εpm results are discussed, as validated in the

Table 6. Two Types of No-Fit-Parameter SCNLE Calculations (ε/kB = [(εp/kB)·90 K]1/2Ns1/2, See the Text) of Dp (in Units of 10−8 cm2/s) for Four Gas Penetrants in Glassy Polysulfone at T = 308 K and the Corresponding Experimental Data (References Refer to the Source) penetrant O2 N2 CO2 CH4

Dp,exp

d/σb

Dap,th

d/σc

Dbp,th

N2 CO2 CH4 C3H8

40

0.267 0.244 0.279 0.318

24.6 7.9 2.1 0.0005

0.272 0.296 0.279 0.387

19.3 0.05 2.1 3 × 10−9

6.0 5.141 1.342 0.00243

44

3.4 0.8544 2.044 0.2744

d/σa

Dap,th

d/σb

Dbp,th

0.252 0.264 0.240 0.276

6.0 7.1 0.7 0.11

0.252 0.271 0.294 0.276

6.0 4.4 4 × 10−4 0.11

a

The d/σ value used is transferred from the optimal size ratio of the O2/PSF system (assumes it equals the optimal size ratio of the O2/PC system (= 0.252, see Table 4)) and dK of oxygen and the penetrant (Table 2) as d/σ = 0.252dK,p/dK,oxygen. bBased on dBH (Table 2) to determine the transferred d/σ value as 0.252dBH,p/dBH,oxygen.

prior two subsections for gas/vapor penetrants. First, on the basis of the four optimal values (d/σ)optimal derived from oxygen diffusion in PVAc, PMMA, PS, and PC (Figure 9 and Table 4), we estimate σ for each polymer in two ways (called type a and type b), either calibrated by dK as σa = dK,O2/(d/σ)optimal or calibrated by dBH as σb = dBH,O2/(d/σ)optimal. For PSF, since its chemical structure is close to PC, we approximate the segment size of these two polymers as the same. Second, we apply the above two estimates of polymer segment size to carry out SCNLE calculations for two groups of experimental systems. The first is nitrogen, carbon dioxide, methane, and propane in PS (Table 5), and the second is oxygen, nitrogen, carbon dioxide, and methane in PSF (Table 6). Without fitting experimental data, we directly employ an a priori estimation of d/σ for each penetrant−polymer pair based on the two choices: d/σ = dK,p/σa (type a) or d/σ = dBH,p/σb (type b). The resultant no-f it-parameter theory results for penetrant diffusivity are denoted as Dap,th and Dbp,th, respectively. We find that the type-a calculations consistently yield predictions in good accord with experiment, while type-b calculations work poorly for CO2 and propane (Tables 5 and 6). Thus, the kinetic diameters summarized in ref 30 serve as a better source of a priori gas/vapor penetrant sphere diameters. The polymer segment size calibrated in the “type-a” manner is denoted as σcali (=σa) and listed in Table 1. All of them are within 20% or less of σest. As Table 4 shows, dK/σcali = 0.245 for water in PVAc, which is closer to the 0.208 value obtained from the fitted experimental data (Figure 10) compared to the a priori estimation (0.298). We do not wish to overemphasize the accuracy of our type-a no-adjustable-parameter results in Tables 5 and 6. On one hand, it is encouraging that the deviation between the theoretical and experimental penetrant diffusivities is always within 1 order of magnitude. On the other hand, the theory and mapping fail to capture some qualitative trends. For example, while the theory correctly reproduces the diffusion constant ordering for penetrant−PS systems (Table 5) as N2 > CO2 > CH4 > C3H8, it predicts the diffusivity of N2 > O2 in PSF at 308 K (Table 6), which is not in agreement with the experimental trend. Figure 11 presents results for methanol in PVAc for two different p−m attraction strengths and two different βK(T) functions. Comparing the four devth‑exp values in Table 4, we find that ε = εpm (as expected) yields smaller theory− experiment deviations than calculations that employ ε = 0. Another sensible trend is that by employing a presumably more

Table 5. Two Types of No-Fit-Parameter SCNLE Calculations (ε/kB = [(εp/kB)·90 K]1/2Ns1/2; See the Text) of the Penetrant Long-Time Diffusivity Dp for Four Gas Penetrants in Glassy Polystyrene at T = 300 K and the Corresponding Experimental Data (References Refer to Source)a penetrant

Dp,exp

All Dp’s are in units of 10−8 cm2/s. bThe d/σ value used is transferred from the optimal size ratio for the O2/PS system (0.254, see Table 4) and dK of oxygen and the penetrant (Table 2) as d/σ = 0.254dK,p/ dK,oxygen. cBased on dBH (Table 2) to obtain the transferred d/σ value as 0.254dBH,p/dBH,oxygen. a

J

DOI: 10.1021/acs.macromol.6b00725 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules

Figure 12 presents optimized theory results for the three aromatic penetrants in a PS matrix. For toluene, experimental

Figure 11. Best fit d/σ (see legend, Table 4) theory results (curves) for methanol (above) and toluene (below) diffusion constants in PVAc as a function of Tg/T using different βK and ε values indicated in the legend (ε = εpm if nonzero). βK(T) = 0.51 + 0.0021(T − Tg) is adapted from experiment.29 Symbols represent the corresponding experimental data.16

Figure 12. Best fit d/σ (see legend, Table 4) theory results (curves) for the diffusion constant for three aromatic penetrants (from top to down: toluene, tetracene, and rubrene) in PS as a function of Tg/T using ε values indicated in the legend (ε = εpm if nonzero). Symbols represent the corresponding experimental data.16,45

accurate βK(T), better predictions are found compared to the βK = 1/3 analogues. Encouragingly, we find dBH/σcali = 0.343 is very close to the best p−m size ratio deduced by fitting experimental data to our theory. We use dBH since dK of methanol is not available, but our findings do suggest that dBH of methanol is a good approximation of its dK value. Overall, the physical consistency and transferability of the microscopic parameters used in our theory for small gas/vapor penetrants seem reasonable. In the Supporting Information, additional theoretical predictions of temperature-dependent diffusivity are presented for four gases (O2, N2, CO2, and CH4) and five polymers (PVAc, PMMA, PS, PC, and PSF); results for all 20 possible gas−polymer combinations are given. We hope they will be soon tested by experiment/simulation.

data in the glassy regime are available,16 and our calculations extend to the glassy regime. The ε = 0 versus ε = εpm comparative studies are performed for toluene and tetracene diffusion. Again, the ε = εpm calculations always result in unphysically small p−m size ratios, though similar devth‑exp numbers are found as for the ε = 0 calculations. Calculations based on ε = 0 consistently lead to physically reasonable d/σ values that are close to a priori estimations. We note that unlike the small gas/vapor penetrant case, the a priori p−m size ratios dest/σest for some systems are closer to the best value obtained by fitting experimental data than dest/σcali. We think this is unsurprising since σcali is calibrated based on oxygen diffusion data, and the large aromatic penetrants interact with polymer segments in a very different manner. The experimental data included in Figure 12 span nearly 10 decades (the most among Figures 9−13), which is well captured by the theory. As a cautionary comment, we note that employing ε = εpm for large organic penetrants can still empirically capture well the experimental data even if an unrealistically small d/σ is used. As

VI. DIFFUSION OF LARGE ORGANIC PENETRANTS We now apply the theory to study the diffusion of larger organic molecules in polymers. Due mainly to the availability of experimental data, three specific aromatic penetrants (toluene, tetracene, and rubrene) are examined. In addition to the methanol/PVAc system, Figure 11 also shows results for toluene diffusion in PVAc based on using two different ε values and two different βK(T) functions. In contrast to small gas/vapor penetrants, but consistent with our analysis and physical arguments in section V, ε = 0 now yields much better agreement with the experiment than calculations based on the choice ε = εpm. Here, we do not mean the “optimized diffusivity curve” is much closer to experimental data for ε = 0 than ε = εpm; in fact, all “optimized diffusivity curves” largely overlap (weak sensitivity to ε) and capture well the experimental data. Rather, we mean the optimal p−m size ratios derived from ε = 0 calculations are chemically reasonable, while those from ε = εpm are unphysically small (