Collective and Single-Chain Correlations in Disordered Melts of

Article pubs.acs.org/Macromolecules

Collective and Single-Chain Correlations in Disordered Melts of Symmetric Diblock Copolymers: Quantitative Comparison of Simulations and Theory Jens Glaser, Jian Qin, Pavani Medapuram, and David C. Morse* Department of Chemical Engineering and Materials Science, University of Minnesota, 421 Washington Ave SE, Minneapolis, Minnesota 55455, United States ABSTRACT: We present a detailed comparison of simulations of disordered melts of symmetric AB diblock copolymers to predictions of the renormalized one-loop (ROL) theory. The behaviors of the structure factor S(q) and of single-chain correlations are studied over a range of chain lengths (N = 16, ..., 128) for two models: one with harshly repulsive pair interactions and another with very soft interactions. The ROL theory is shown to provide an excellent description of the dependence of S(q) on chain length and thermodynamic conditions for both models, even for very short chains, if we allow for the existence of a nonlinear dependence of the effective interaction parameter χe upon the strength of the AB repulsion. The decrease in peak wavenumber q* with increasing χe is shown to be unrelated to changes in single-chain correlations. Results for all quantities are consistent with the hypothesis that the ROL theory gives an exact description of the dominant 6(N̅ −1/2) corrections to RPA and random-walk predictions in the limit of infinite chain length N. the ratio cR3/N of the volume R3 pervaded by a random walk of size R = N1/2b to the volume N/c actually occupied by the polymer. The Fredrickson−Helfand theory and its immediate descendants were designed to describe only the dominant effects of strong composition fluctuations very near the order− disorder transition (ODT) of long, nearly symmetric diblock copolymers. After Fredrickson and Helfand, several authors who tried to construct more complete theories of corrections to the SCFT/RPA theory reported serious technical difficulties11,12 that delayed further progress. The most serious conceptual problems arose from the fact that predictions obtained from highly coarse-grained models were found to be very sensitive to the treatment of short wavelength correlations, at length scales that the underlying models were never designed to describe. A more recent generation of “renormalized oneloop” (ROL) theories13−21 have overcome this problem by introducing a renormalization scheme in which all contributions that depend sensitively on monomer scale correlations are explicitly absorbed into corrections to the values of the χ parameter and statistical segment lengths. This and other technical refinements have allowed the development of much more rigorous and complete theories of corrections to the random walk model of single-chain statistics in polymer melts,14−16 and of corrections to the random walk and random phase approximations in polymer blends13,18,19 and in diblock copolymer melts.18,21

1. INTRODUCTION For many years, theoretical study of block copolymers has been dominated by the self-consistent field (SCF) theory of phase behavior1 and by the random phase approximation (RPA) theory of composition fluctuations and scattering.2 These theories are intimately related: The RPA for the structure factor S(q) in a homogeneous phase is derived by using a SCF approximation to estimate the free energy cost of infinitesimal composition fluctuations. The physical approximations underlying the SCF and RPA theories rely upon the existence of a separation of length scales in a polymer liquid, between the size of a monomer and that of an entire polymer coil. It has long been believed, partially on intuitive grounds, that the SCF and RPA theories should become more accurate with increasing chain length. Until recently, however, there was very little clear evidence for this from simulations or experimental data. Understanding of the limitations of SCFT and the RPA for block copolymers has long relied primarily on the work of Fredrickson and Helfand.3 These authors adapted earlier work by Brazovskiı̆4 on weakly first-order crystallization transitions to describe the order−disorder transition (ODT) of nearly symmetric diblock copolymers. The Fredrickson−Helfand theory3 and subsequent closely related work by several authors5−10 predict corrections to the SCFT and RPA theories whose magnitude is controlled by a dimensionless number N̅ = N(cb3)2, sometimes referred to as the invariant degree of polymerization. Here, N is the degree of polymerization (number of monomers per chain), c is the monomer concentration, and b is a statistical segment length. Predicted corrections decrease with increasing N̅ , thus recovering SCF and RPA predictions in the limit N̅ → ∞. The parameter N̅ is a measure of chain overlap, given by the square N̅ = (cR3/N)2 of © 2014 American Chemical Society

Received: August 12, 2013 Revised: December 2, 2013 Published: January 13, 2014 851

dx.doi.org/10.1021/ma401694u | Macromolecules 2014, 47, 851−869

Macromolecules

Article

In this paper, we compare results of extensive simulations of disordered symmetric diblock copolymers to predictions of the ROL theory.18,21 We compare simulations of two coarsegrained models to ROL predictions and to each other. One model (model H) uses a harshly repulsive Weeks−Chandler− Anderson pair interaction. This model is similar to that of Kremer, Grest, and co-workers.22,23 The other (model S) uses a softer interaction and higher monomer concentration, and is almost identical to the model introduced by Groot and Warren24 for use in dissipative particle dynamics (DPD) simulations. A preliminary comparison of results of model H to ROL predictions for S(q) was presented in ref 25. Some data for S(q) from both of these models were also presented in ref 26 for a narrower range of chain lengths than discussed here. The analysis presented in ref 26, however, focused exclusively on testing the consistency among results obtained from simulations of different simulation models with matched values of N̅ , rather than on testing the ROL theory. The present article is a comprehensive report of tests of the ROL theory for disordered symmetric diblock copolymer melts. Here, we combine data from models H and S to cover a range of N̅ up to 2000, to test the consistency of results for two models, while also testing predictions for a wider range of physical properties than in our previous work. This includes the first tests of predictions of the ROL theory for single-chain properties in block copolymer melts. The article is organized as follows: Section 2 presents notation and theoretical background. Section 3 presents the models and simulation methods used here. Section 4 presents the analysis of homopolymer simulations that we use to obtain initial estimates for the statistical segment and the Flory− Huggins interaction parameter χ. Section 5 focuses on the behavior of S(q) in the limit of thermodynamically “ideal” copolymers, with χ = 0. Section 6 presents an analysis of results for the parameter dependence of the peak intensity S(q*) in systems with nonzero χ. The analysis differs from that given for model H in ref 25 in that we allow for the existence of a nonlinear dependence of χ on the strength of repulsion between A and B monomers. This difference is found to dramatically improve the level of agreement between the ROL theory and simulations of short chains. Section 7 tests ROL predictions for the peak wavenumber q*. Section 8 presents a more detailed analysis of results for S(q) based on the use of the Ornstein−Zernike equation to explicitly separate S−1(q) into intra- and intermolecular contributions. One important conclusion of this analysis is to confirm an earlier prediction21 that the decrease in q* with increasing χN is essentially unrelated to changes in the coil dimensions. Section 9 presents results for radii of gyration of the chain and of individual blocks. Conclusions are summarized in section 10.

Sij(q) ≡

∫ dr eiq·r⟨δci(r)δcj(0)⟩

(1)

where δci(r) is the deviation of the concentration ci(r) of i monomers from its spatial average. In a dense liquid, long-wavelength fluctuations in the total monomer density c(r) ≡ cA(r) + cB(r) are much smaller than corresponding fluctuations in the composition field ψ (r) ≡

1 [δcA(r) − δc B(r)] 2

(2)

It is thus convenient to define a scalar structure function S(q) ≡

∫ dr eiq·r⟨ψ (r)ψ (0)⟩

(3)

to characterize pure composition fluctuations. The scalar structure factor S(q) for a disordered diblock copolymer exhibits a maximum value S(q*) at a peak wavenumber q* of order the inverse polymer coil size. No liquid is completely incompressible. In a slightly compressible liquid, we expect the matrix Sij(q) at q of order the inverse coil size to have one large eigenvalue, with an eigenvector ε− = (1,−1) that corresponds to fluctuations of ψ(r) at fixed concentration and a second much smaller eigenvalue with an orthogonal eigenvector ε+ = (1,1) that corresponds to concentrations in total concentration c(r). The function S(q) gives the variance of the composition mode. The vectors ε+ and ε− defined above are also exact eigenvectors of Sij(q) in the special case considered in this paper, as a result of symmetry. The simulations presented here all use symmetric models of diblock copolymers, with equal number of A and B beads of identical structure. In this case, the system is symmetric under exchange of the labels of A and B beads. As a result, eigenvectors of Sij(q) must be either even or odd under this symmetry, and so must be given exactly by ε+ and ε− for all q, even in a compressible liquid. Let Ωij(q) denote the intramolecular contribution to Sij(q), i.e., the contribution to Sij(q) for q ≠ 0 that arises from correlations between pairs of monomers of types i and j within the same chain. This can also be written as a single-chain average c Ωij(q) = ∑ ∑ ⟨eiq·[R(s) − R(s ′)]⟩1‐chain N s ∈ i s ′∈ j (4) in which c/N is the number concentration of polymers, s and s′ are monomer indices that are restricted to lie within the i and j blocks of the copolymer, respectively, R(s) is the position of monomer s, and ⟨···⟩1‑chain denotes an average over conformations of a single chain. Prediction of coarse-grained theories that are based upon a model of polymer as continuous random walks, such as the RPA and ROL theories, must be invariant under a rescaling of N and c that corresponds to a redefinition of what we mean by a “monomer”. RPA and ROL predictions for S(q) are thus given as predictions for a normalized inverse structure factor

2. THEORETICAL BACKGROUND In this section, we define notation for quantities that are measured in our simulations and review theoretical predictions for these quantities. We consider a homogeneous liquid containing M symmetric AB diblock copolymers in a volume V. Each chain contains N monomers, N/2 of each type. The total monomer concentration is c ≡ MN/V. Correlation Functions. In a liquid containing two types of monomer, A and B, indexed by a type indices i and j, we may define a 2 × 2 correlation function matrix with elements

H(q) ≡ cNS −1(q)

(5)

which is invariant under this rescaling. Corresponding predictions for Ωij(q) may be given as predictions for Ωij(q)/ cN. Random Phase Approximation (RPA). The RPA for S(q) yields2 852


Macromolecules H0(q) = F0(q) − 2χe N

Article

diblock copolymers by Leibler,2 generalized to allow for nonGaussian single-chain statistics. The function χa(q), which we refer to as an “apparent” interaction parameter, is a generalization of the RPA interaction parameter χe. The key difference between the RPA and incompressible OZ equations is that the RPA is a predictive theory, whereas the incompressible OZ equation is merely descriptive. The RPA can be recovered from eq 9 by assuming random-walk statistics, which yields an explicit expression for F(q) as a function of q and N, and assuming that χa(q) is independent of q and N. Equations 9−11 provide definitions of F(q) and χa(q) but do not predict anything about how these quantities depend on q and N. Given measurements of S(q) and Ωij(q), one can always use eq 10 to calculate F(q) and then use eq 9 to extract χa(q). One-Loop Theories. In what follows, we compare simulations results to the ROL theory for diblock copolymers,18,21 which predicts corrections to both the RPA for S(q) and to the random walk model of single-chain statistics. The ROL theory considered here is the most rigorous variant of a family of coarse-grained field theories that all descended from the Fredrickson−Helfand theory,3 as discussed below. The original Fredrickson−Helfand theory3 and closely related work by Dobrynin and Erukhimovich 8,9 gave predictions for the dominant singular behavior of the peak intensity S(q*) in the disordered phase near the ODT but did not allow for the possibility of a shift in q*. Subsequent work by Olvera de la Cruz5,6 and Barrat and Fredrickson7 generalized the theory so as to allow for a shift in q*. These early theories for S(q) all relied on a somewhat ad hoc approximation in which the SCFT free energy is used as an effective Hamiltonian within a functional integral of a composition order-parameter field. We refer to this as a mean-field effective Hamiltonian (MFEH) approximation. A recent analysis20 compared predictions of the MFEH approximation to those of an auxiliary-field approach in which the partition function is instead expressed rigorously as a functional integral of a pair of complex chemical potential fields, as done in both the ROL theory13,18 and in recent field-theoretic numerical simulations by Fredrickson and co-workers.32−34 It was shown that the MFEH approximation, in the absence of further mathematical approximations, introduces a q-dependent error of order 6(N̅ −1/2) in predictions for cNS−1(q) but has no effect on predictions of the dominant singular dependence of cNS−1(q) on χN in the strong-fluctuation regime near the ODT. Barrat and Fredrickson7 and Erukhimovich and Dobrynin10 have also long ago presented predictions for changes in the size of individual polymers and individual blocks in a diblock copolymer melt near the ODT. These calculations of single chain dimensions, and more recent ROL calculations of singlechain properties,16,18,19,21 are all based on an simple theory that was first considered by Edwards for solutions.35,36 In this approach, each chain is taken to interact with itself via a screened interaction,35,37,38 and the effect upon single-chain statistics is calculated to first order in the strength of the screened interaction. We note in passing that the expression for the screened interaction reported by Barrat and Fredrickson (eq 27 in ref 7) is different from that obtained by other workers10,18,38 and appears to us to be incorrect. The calculations of Erukhimovich and Dobrynin10 for chain and block dimensions in both symmetric and asymmetric copolymers appear, however, to start from an unrenormalized

(6)

in which F0(q) is a function that reflects the statistics of random walk polymers, which was first calculated for block copolymers by Leibler,2 and χe is an effective interaction parameter. Here and hereafter, we use a subscript 0 to denote quantities such as H0(q) ≡ cNS0−1(q), S0(q), and F0(q) that are calculated using the RPA or random-walk approximations. In its most general form, the RPA allows χe to exhibit an arbitrary dependence on the temperature, concentration, and parameters of the Hamiltonian but requires that χe be independent of q and N. The analytic expression for F0(q) for symmetric diblock copolymers (which we do not display) can be expressed as a dimensionless function F0(x) of a quantity x ≡ (qRg0)2, where Rg0 = (Nb2/6)1/2 is the radius of gyration of a Gaussian chain with statistical segment length b. For symmetric copolymers, this function has a minimum value of F0(x*) = 20.99 at x* = 3.785.3 This gives a peak in S(q) at q*0 = 1.946/Rg0 with H0(q*0 ) = 2[(χN)s − χeN], where (χN)s = F0(x*)/2 = 10.495 is the predicted spinodal value of χeN. In what follows, as in previous work,19,21,26 we sometimes characterize the maximum peak intensity S(q*) measured in a simulation by giving a value for an “apparent” interaction parameter χ*a . We define χ*a to be value of the RPA interaction parameter χe that would be inferred by fitting the actual peak intensity to the RPA. This quantity is defined by the equation cNS −1(q∗) ≡ 2[(χN )s − χa* N ]

(7)

in which (χN)s = 10.495 is the RPA prediction for the spinodal value of χN. Ornstein−Zernike (OZ) Equations. Our starting point for analyzing deviations from the RPA is an Ornstein−Zernike equation for polymer fluids27−31 Sij−1(q) = Ωij−1(q) − Cij(q)

(8)

−1

in which Ωij (q) is the matrix inverse of the 2 × 2 matrix Ωij(q), and Cij(q) is the so-called direct correlation function. The direct correlation function Cij can be interpreted, somewhat heuristically, as an effective interaction between monomers of type i and j. In an incompressible system, or an AB system that is symmetric under relabeling of A and B monomers, the OZ equation for Sij(q) may be reduced to an equation for the inverse scalar correlation function S−1(q). In these cases, S−1(q) = ∑ijε−i ε−j Sij−1(q), where ε− = (1,−1). Applying this to eq 8, and multiplying by cN, we obtain an expression H(q) = F(q) − 2χa (q)N

(9)

closely analogous to the RPA, in which F(q) ≡ cN ∑ εi−εj−Ωij−1(q) ij

= cN

Ω11(q) + Ω12(q) + Ω 21(q) + Ω 22(q) Ω11(q)Ω 22(q) − Ω12 2(q)

(10)

and χa (q) ≡

c 2

∑ εi−εj−Cij(q) ij

(11)

We will refer to eq 9 as the incompressible Ornstein−Zernike (OZ) equation. The function F(q) defined in eq 10 is a generalization of the function that was calculated for Gaussian 853


Macromolecules

Article

in which N̅ = N(cb3)2. Here, F1, X1, and H1 are dimensionless functions of qRg0 and χ*a N that satisfy H1 = F1 − 2X1. These functions are defined by Fourier integrals,18 that we evaluate by numerical integration. Because the ROL theory predicts corrections to the RPA that decrease as N̅ −1/2, the RPA is recovered in the limit N̅ → ∞. ROL predictions for F(q) are obtained from corresponding predictions for Ωij(q)/cN. ROL predictions for the radii of gyration of polymers and individual blocks, which are discussed in section 9, have been inferred from the low-q behavior of the elements of Ωij(q)/cN. The notation used in eq 13 is intentionally ambiguous, since it indicates that the one-loop corrections given by the functions F1, X1, and H1 depend on an invariant interaction parameter “χN” but does not specify which of several possible definitions of χ should be used to calculate this parameter. This ambiguity reflects the fact that there are actually two versions of the theory, which we will refer to as the “perturbative” and “selfconsistent” variants. In the “perturbative” theory, which is the version discussed in the original discussion of renormalization,18 the input parameter to the Fourier integrals that define these one-loop corrections is taken to be the RPA interaction parameter χe, which is independent of N. In the “selfconsistent” version of the theory, which was introduced in subsequent work on blends19 and diblock copolymers,21 this parameter was taken to be the “apparent” interaction χ*a defined by eq 7, which is determined by the self-consistently determined peak intensity. In the self-consistent ROL theory, H(q) is thus given by an equation of the form

theory equivalent to that underlying the ROL theory considered here. These early theoretical predictions of corrections to the RPA for S(q)3,5,7 and to the random-walk model of single-chain statistics7,10 all relied on a set of physical and mathematical approximations that were expected to capture only the dominant corrections to the RPA in the region of strong fluctuations near the ODT, where universal 6(N̅ −1/2) corrections are dominated by composition fluctuations with wavenumbers very near q*. These theories were thus intended to be valid only in this limit. The ROL theory is instead expected to remain valid for sufficiently long chains (N̅ ≫ 1) arbitrarily far from the ODT, i.e. at arbitrary values of χN in the disordered phase. We have argued elsewhere18,20 that the ROL yields an exact expression for the universal 6(N̅ −1/2) corrections to the RPA within an expansion in powers of N̅ −1/2. One goal of this work is to test more thoroughly whether simulation results are consistent with this claim. In general, predictions of the ROL theory reduce to those of corresponding earlier theories in the strong-fluctuation, large χN limit for which the earlier theories were designed.20 Numerical comparisons of predictions of the ROL and FH theories for S(q*) to each other21 and to simulation data25 have been presented previously. A comparison of both ROL and Fredrickson−Helfand predictions to simulations of model H for symmetric diblock copolymers of N = 64 beads25 showed that both theories agree well with simulation results at very large values of χN but that Fredrickson−Helfand theory fails further from the ODT, giving the wrong sign for corrections to the RPA at low values of χN. The ROL theory instead remains quite accurate down to χN = 0. The accuracy of the ROL theory for symmetric diblock copolymers in the limit χN = 0 is a result of the fact that its predictions for single-chain statistics become equivalent in this limit to those of the very accurate ROL theory for corrections to the random statistics for homopolymers in a dense homopolymer melt by Beckrich et al.16 and the fact that the corrections to the RPA for S(q) are controlled in this limit by corrections to single chain statistics.21 ROL predictions for the dependence of q* upon χN are also expected on theoretical grounds20 to reduce near the ODT to those given by the generalization of the FH theory that allows for a shift in q*.5,7 Differences between the ROL theory and the FH theory become less important with increasing N̅ but remain crucial for the modest values of N̅ used in many of our simulations.21 Because the full ROL theory includes the older theories discussed above as special cases, but has a much wider range of validity, particularly for modest values of N̅ , we have chosen to only present comparisons of simulation data to predictions of the full ROL theory in this paper. ROL predictions for S(q) are most conveniently expressed as expressions for the invariant quantities H(q) ≡ cNS−1(q), F(q), and X(q) ≡ χa (q)N

H(q) = H0(q) + N̅ −1/2H1(qR g0 , χa* N )

in which χ*a is related to the minimum value H(q*) = minq[H(q)] by eq 7, yielding an implicit equation for H(q*). The self-consistent version of the ROL theory is designed to mimic the behavior of the Brazovskiı̆−Fredrickson−Helfand theory near ODT. The use of a self-consistently determined value of χ*a N is neither required nor prohibited by the logic of the original ROL theory.18 It is a physically motivated modification that was introduced21 in an attempt to combine the virtues of the ROL theory, which provides a systematic expansion that is expected to be accurate where corrections to the RPA are small, with those of Brazovskiı̆−Fredrickson− Helfand theory, which was designed to capture the dominant behavior near the ODT, where corrections to the RPA become very large. The way we introduce self-consistency is also closely analogous to the “Hartree” approximation used by Wang13 in his ROL theory of fluctuations in near-critical homopolymer blends. In most of this paper, we will compare simulation results to the self-consistent ROL theory. The only exception is section 5, where we use a perturbative version of the theory with χe = 0 to describe “ideal” copolymers with χe = 0. Correlations in Real Space. For some purposes, it is useful to have an approximate expression for the real-space correlation function S(r) = ⟨ψ(r)ψ(0)⟩ as well as for its Fourier transform S(q). Near the ODT, the function S(r) exhibits exponentially decaying oscillations with an oscillation period d* ≡ 2π/q*.28 In this limit, we may approximate the RPA prediction for S(r) (where r ≡ |r|) by expanding the RPA prediction for H(q) = cNS−1(q) to harmonic order in q − q*, giving H(q) ≃ H(q*)[1 + (q − q*)2ξ2]. By using this approximation for H(q) to

(12)

with H(q) = F(q) − 2X(q). ROL predictions are of the form F(q) = F0(q) + N̅ −1/2F1(qR g0 , χN ) X(q) = χe N + N̅ −1/2X1(qR g0 , χN ) H(q) = H0(q) + N̅ −1/2H1(qR g0 , χN )

(14)

(13) 854


Macromolecules

Article

For each model and each value of N, we carry out simulations at fixed values of kBT, εAA = εBB, σ, κ, and l0 over a range of values of a parameter α ≡ εAB−εAA (19)

evaluate the inverse Fourier transform of S(q), and retaining only the dominant contribution for ξ ≫ d*, one obtains28 S(r ) ≃

Ac 2 sin(q*r ) −r / ξ e (N̅ τ )1/2 q*r

(15)

which we use to control the thermodynamics of mixing. All data shown here are taken within the disordered phase. We have identified values of α at the order−disorder transition for all of the models presented here and will discuss this phase transition elsewhere. The value of κ used in model S was chosen as to obtain a value for the quantity cb3 in model S that is larger than that for model H by exactly a factor 2. This choice yields matched values of N̅ ≡ N(cb3)2 for simulations of models H and S for values of N that differ by a factor of 4. Models H64 and S16 thus have matched values of N̅ ≃ 240, while H128 and S32 both have N̅ ≃ 480. Consistency of results obtained with models H and S can be tested by directly comparing results obtained from these two pairs of simulations for various physical properties, as done previously26 for q* and S(q*). All simulations were carried out in NVT ensemble in periodic cubic simulation cells of volume V = L3. In all of the new simulations presented here, the number of molecules M was chosen so as to yield a cell length L ≃ 3d0* approximately 3 times the RPA layer spacing of d*0 ≡ 2π/q*0 , or L ≃ 9.67Rg0, in systems that have exactly the chosen monomer concentration for each model. (The length L and the monomer concentration c = MN/L3 cannot both be set to exactly specified values simply because M must be an integer.) A list of system sizes used for different chain lengths is given in Table 1.

in which the ratio of the decay length ξ to the layer spacing d* = 2π/q* is given for symmetric copolymers by3 ⎛ 2F ″(x*)(x*)2 ⎞1/2 ⎛ 0.70 ⎞1/2 ξ ⎟ =⎜ =⎜ ⎟ d* ⎝ (2π )2 H(q*) ⎠ ⎝ H(q*) ⎠

(16)

where F″(x*) = d2F(x)/dx2|x* = 0.9624 and A = [27x*/ (π2F″(x*)]1/2 = 3.28. A similar analysis can be applied to the true correlation function, simply by using the actual minimum value of H(q) and the second derivative of this function with respect to q rather than RPA predictions for these quantities in the Taylor expansion of H(q). Corrections to the RPA are known to have a large effect on the minimum value of H(q) near the ODT, but smaller effects on the curvature of this function and on the value of q*. (Recall that the original Fredrickson−Helfand theory3 predicted no change in either of these two quantities.) We thus expect eqs 15 and 16 to provide useful approximations for ξ and S(r) even near the ODT, as long as we define H(q*) using the measured peak intensity, rather than the RPA prediction for this quantity. When used in this way, eq 16 provides a useful starting point for understanding finite size effects in simulations.

3. SIMULATION MODELS AND METHODS We have conducted NVT simulations of two coarse-grained simulation models that we refer to as model H (for “hard”) and model S (for “soft”). Simulations of each model have been carried out using four chain lengths, N = 16, 32, 64, and 128. Throughout our discussion of results, we use codes consisting of the letter H or S followed by a value of N to identify specific simulations. Thus, for example, H64 refers to a simulation of model H using chains of 64 beads. In model H, all nonbonded beads interact via a purely repulsive Weeks−Chandler−Anderson (WCA) pair potential. The nonbonded potential between i and j particles is Vij(r ) = 4εij[(r /σ )−12 − (r /σ )−6 + 1/4]

Table 1. Model Dimensions, Chain Length N, Number of Molecules M and Cubic Box Length La model S N̅

M

L/σ

N

N̅

M

L/σ

16 32 64 128

238.9 477.7 955.4 1910.9

952 1347 1904 2694

17.19 24.31 34.37 48.62

16 32 64 128

59.8 119.6 239.2 478.4

470 665 940 1330

22.07 31.21 44.13 62.42

The monomer concentration is exactly c = 0.7σ−3 for model H and c = 3.0σ−3 for model S, while M is chosen to so as to give a length L of approximately 3 times the RPA layer spacing 2π/q*0 . a

(17)

All of the new simulations presented here are GPUaccelerated molecular dynamics simulations conducted using the HOOMD-blue39 code. All MD simulations presented here use a time step of Δt = 0.005 (in LJ units) for both models, using a Nosé-Hoover thermostat. We have confirmed that this yields results equivalent to those obtained using smaller time steps and to earlier hybrid Monte Carlo simulations.25,26 Figure 1 shows an example of the behavior of S(q) for a simulation of model S32 at several values of α. Each discrete point corresponds to data from a family of symmetry related wavevectors that are commensurate with the simulation unit cell. The lines through the points in this graph are fits to an empirical model of the form

for separations r less than a cutoff distance rc = 2 σ, and Vij(r) = 0 for r > rc. Consecutive beads within each chain interact via a harmonic bond potential Vbond(r) = κ(r − l0)2/2. In this model, the nonbonded pair interaction is suppressed for bonded particles. All simulations of this model use parameters ε = εAA = εBB = kBT, κ = 400kBT/σ2, l0 = σ, and a monomer concentration c = MN/V = 0.7σ−3, as in our previous work.25,26 In model S, all particles (including bonded particles) interact via a nonbonded pair interaction 1/6

Vij(r ) = εij(σ − r )2 /(2σ 2)

model H

N

(18)

for r < σ and Vij(r) = 0 for r > σ. Bonded particles interact via an additional harmonic bond potential with zero rest length (l0 = 0), giving Vbond(r) = κr2/2. All simulations of this model use parameters ε = εAA = εBB = 25kBT and κ = 3.406kBTσ−2, with a monomer concentration c = 3.0σ−3.

cNS −1(q) ≃ F0(qR g 0) + a + bq2 + cq 4

(20)

in which F0(qRg0) is the RPA prediction for F(q), and in which a, b, and c are adjustable parameters. We have used this fit to estimate values of the peak wavenumber q* and the maximum 855


Macromolecules

Article

Figure 1. Fit of data for S(q) of model S, N = 32 and α = 0.0, 0.05, 1.52, and 1.9 (from bottom to top), to eq 20 (error bars included). The fits yield (bottom to top) q*Rg0 = 1.911, 1.85, 1.65, 1.75 and cS(q*) = 1.60, 4.86, 17.43, 55.80, or, equivalently, χ*a N = 0.52, 7.20, 9.59, 10.20.

S(q*) for each simulation but attach no physical meaning to the fit parameters a, b, and c. The fit was applied to structure factor values S(q) with 0 < qRg0 ≤ 6.78 in all models, corresponding to wave vectors with Miller indices i, j, k ≤ 6. Finite size effects are obviously a source of concern for simulations with such small unit cells. In earlier simulations of model H,25 results for q* and S(q*) were reported for at least two different cell sizes with L/d*0 = 2−4 for each chain length, but results obtained using different systems sizes were found to be indistinguishable except for the longest chains very near the ODT. Equation 16 for ξ provides a partial explanation of this somewhat surprising lack of finite size effects. Finite size effects are expected to be small when the correlation length ξ remains much less than box size L. Equation 16 predicts a correlation length ξ greater than one layer spacing d* only for H(q) < 0.7. Because we study the entire range from χe = 0 (or α = 0) up to the ODT, most of our data automatically fall outside the region very near the ODT in which ξ > d*. Equation 16 also, however, makes it clear that caution is required for quantitative analysis of data very near the ODT. The peak intensity S(q*) for the largest value of α shown in Figure (1) yields a value of H(q*) ≃ 0.6, or χa*N = 10.2, corresponding to ξ ≃ d*. This example is near the limit beyond which we do not believe it is possible to accurately fit the peak using data from a cell of size L ≃ 3d*. To more carefully quantify finite size effects, we have compared simulations of model S64 with M = 1904 molecules, the system size used in most of simulations of this model, to simulations conducted with 4 times as many molecules (M = 7616) at the same monomer concentration. Figure 2 shows an example of the raw data and fits for S(q) obtained from these pairs of simulations at two different values of α = 0.75 (upper panel) and 0.95 (lower panel). These two values of α are somewhat below and above the value of α at which we begin to see significant finite size effects, respectively. The results illustrate that our fitting procedure yields different values for S(q*) from different system sizes only above a critical value of α. Figure 3 shows fitted values of S(q*)/c obtained from these simulations of different size systems, plotted vs α. The horizontal dotted line shows the value of S(q*)/c = 91.4 for which H(q*) = 0.70, corresponding to ξ ∼ d*. Significantly different values for S(q*) are obtained only above this bound. This simple criterion seems to adequately identify the onset of finite size effects for systems of size L ≃ 3d* and so is used in

Figure 2. Simulated S(q) data and resulting fits for two different system sizes M = 1904 and M = 4 × 1904 (M: number of molecules) of model S64 at α = 0.75 (upper panel) and 0.95 (lower panel). The value of α in the lower panel falls in the range for which we expect finite-size effects, based on an analysis of the peak scattering intensity S(q*)/c (Figure 3 and main text).

Figure 3. Maximum structure factor S(q*) vs α for different system sizes of model S64. The dashed line represents the limiting value cNS−1(q*) ≲ 0.35, beyond which we expect to observe finite-size effects.

our analysis to identify data that may be corrupted by finite size effects.

4. PARAMETER CALIBRATION The RPA for S(q) contains two phenomenological parameters, b and χe. The ROL theory is built on the RPA and contains the same two parameters. The ROL theory yields random walk statistics only in the limit N̅ → ∞, and reduces to the RPA for S(q) in the same limit. When comparing simulation results to ROL predictions, the RPA parameters b and χe should thus be understood to be properties of a hypothetical system of 856


Macromolecules

Article

infinitely long chains. Values of these parameters can thus be determined by extrapolating results of simulations of finite chains to infinite chain length. We obtain estimates for both b and χe at small values of α by extrapolating results obtained from simulations of homopolymer (α = 0) reference systems, as discussed below. 4.1. Statistical Segment Length. In what follows, we use b to denote the statistical segment length of a hypothetical system of infinite chains. This is given for homopolymers (or diblocks with α = 0) by the limit b2 ≡ lim 6R g 2(N )/N N →∞

(21)

in which Rg2(N) is the measured radius of gyration in a dense melt of chains of length N. The ROL theory predicts19 that Rg2(N) should vary with N as R g 2(N ) =

⎤ γ Nb2 ⎡ 1.42 + ...⎥ ⎢⎣1 − 1/2 + ⎦ 6 N N̅

(22)

The dominant 6(N̅ −1/2) correction is a universal correction that arises from the incompressibility of the liquid, as first discussed by Wittmer, Semenov, Beckrich, and co-workers 14−16 , who also gave an analogous expression for the root-meansquare end-to-end length of a finite chain segment. The predicted prefactor of 1.42 for the radius of gyration of monodisperse chains was calculated by Qin and Morse.19 The subdominant 6(1/N ) contribution is nonuniversal because it includes corrections that arise from the discrete structure of the chain, and thus has a coefficient γ that must be treated as a model-dependent fitting parameter. The upper and lower panels of Figure 4 show results for the ratio 6Rg2(N)/N vs N̅ −1/2 from homopolymers of length N = 16, 32, 64, and 128, for models S and H, respectively. The solid line in each figure is a fit to eq 22 in which the coefficients b and γ have been adjusted to fit the data, but in which the theoretically predicted coefficient has been used for the N̅ −1/2 term. The dashed line is the asymptotic prediction obtained by setting γ = 0 after determining b. The predicted asymptotic slope agrees very well with the data, confirming that the ROL theory accurately describes the dominant corrections to the random walk model in homopolymer melts. These fits yield values of b = 1.404σ for model H and b = 1.088σ for model S. As already noted, this yields values of cb3 = 1.937 for model H and cb3 = 3.864 for model S that (by design) differ by a factor of 2. In the ROL theory, b is a phenomenological parameter whose value is understood to be sensitive to details of local fluid structure. When analyzing simulations of each model over a range of values of α, it would be logical to allow b to depend on α. The above analysis of homopolymer simulations only gives the value of b at α = 0. The dependence of b on α is expected to be weak within the disordered phase of systems of long chains, however, because α is constrained to remain below a value 6(1/N ) to avoid microphase separation. The simplest way to model this would be to assume a function b2(α) ≃ b2(0)[1 + dα] and to treat the coefficient d as another phenomenological parameter. For α of 6(1/N ), this could produce 6(1/N ) fractional changes in b and corresponding nonuniversal 6(1/N ) corrections to predictions for q* and Rg. Such corrections would be smaller than the dominant N̅ −1/2 corrections for long chains but are potentially measurable.

Figure 4. Fit of the squared and normalized radius of gyration Rg2(N)/ N at α = 0 as a function of N̅ −1/2 for models S (upper panel) and H (lower panel). Data for each model is shown for N = 16, 32, 64, and 128. The solid curve is a fit to eq 22, in which b and γ are treated as fitting parameters. The dashed curve is the asymptote obtained by setting γ = 0.

In a previous study,26 we compared results of models H and S for q*Rg0 while using the homopolymer (α = 0) value for b to calculate Rg0 ≡ (Nb2/6)1/2 and found (perhaps surprisingly) almost perfect consistency between results obtained with the two models in systems with matched values of N̅ . This indicates that the α dependence of b in these two models must be either very weak or very similar. The level of agreement of results for q*Rg0 for model H with the ROL theory25 suggests, moreover, that any such dependence must be quite weak. Any attempt to estimate the intrinsic α dependence of b by measuring Rg in blends with α > 0 would require a more subtle extrapolation procedure than used above to estimate b at α = 0 because a series of blends with any fixed positive value of α would phase separate above a critical value of N, complicating extrapolation to infinite N at a fixed value of α. For simplicity, in this paper, we have chosen to ignore the possibility of a slight dependence of b on α and again use the α = 0 value of b throughout our analysis, while keeping an eye out for discrepancies that might be attributed to this approximation. 4.2. Perturbation Theory for χe. Attempts to make precise comparisons of simulations of block copolymers and polymer blends have often been constrained by the lack of a rigorous procedure for estimating the RPA interaction parameter χe. In ref 40, Morse and Chung, refining an earlier proposal by Müller and Binder,41 provided a partial solution to this by using thermodynamic perturbation theory to estimate χe(α). Reference 40 analyzed a simple perturbation theory for the free energy of a symmetric blend of A and B homopolymers with a pair potential of the form Vij(r) = εiju(r) used in both of the models studied here, where u(r) is a model-dependent 857


Macromolecules

Article

dimensionless function, using α = εAB − εAA as an expansion parameter. The excess free energy per monomer fex for such a blend is given to first order in α by fex(α) = αz(N)ϕAϕB, where z(N) is a dimensionless measure of the degree of contact between A and B monomers in the homopolymer (α = 0) reference state. This coefficient z(N), known as the effective coordination number, is given by ratio z(N) ≡ einter/ε, where einter is the average intermolecular pair energy per monomer in a homopolymer melt with V(r) = εu(r). The coefficient z(N) approaches a finite limit z∞ as N → ∞ but is expected to exhibit a weak dependence on chain length for finite chains. Morse and Chung40 predicted a dependence ⎡ (6/π )3/2 δ⎤ + ⎥ z(N ) ≃ z∞⎢1 + 1/2 N⎦ ⎣ N̅

The dashed line show the asymptotic behavior obtained by setting δ = 0, which shows that the data is in excellent agreement with the theoretically predicted asymptotic slope. From these fits, we obtain values of z∞ = 0.2965 for model H and z∞ = 0.237 for model S. By comparing the results of this perturbation theory to ROL predictions for polymer blends, Morse and Chung40 concluded that the RPA interaction parameter χe(α) is given to first order in an expansion in powers of α by χe(α) ≃ z∞α/(kBT), plus corrections of 6(α 2). In what follows, we use the notation z α χe(1) (α) ≡ ∞ kBT (24) to denote the resulting first order approximation for χe(α). This approximation was used in ref 25 to analyze peak intensity data for model H.

(23)

The dominant 6(N̅ −1/2) correction to the limiting value has a universal coefficient of (6/π)3/2 that is predicted by the ROL theory. The 6(1/N ) correction, however, has a coefficient δ that is nonuniversal because it includes model-dependent contributions that arise from the chain ends and so must be treated as an additional adjustable parameter. Figure 5 shows simulation results for values of z(N) vs N̅ −1/2 from simulations of homopolymer melts of different chain

5. IDEAL DIBLOCKS (χe = 0) In this section, we compare simulation results and ROL predictions for S(q) in the case of thermodynamically “ideal” diblocks, with α = 0. In this limit, A and B monomers are physically indistinguishable but are labeled for the purpose of defining S(q). This is an idealization of a small-angle neutron scattering (SANS) experiment that measures the correlation hole in a melt of diblock copolymer containing blocks of hydrogenated and deuterated versions of the same monomer, idealized by ignoring the small but nonzero χ between hydrogenated and deuterated monomers. It has been shown analytically, in Appendix B of ref 21, that in the special case of perfectly symmetric diblock copolymers with α = 0, χa(q) is exactly zero for all q. This implies that, in this case, cNS −1(q) = F(q)

(25)

where F(q) is defined using the exact single-chain correlation functions in a dense homopolymer liquid. Equation 25 is an exact result that is a consequence of the symmetry of the system under interchange of the labels of A and B monomers and thus applies only to this special case. In this case, however, it implies that S(q) is completely determined by single-chain statistics and that deviations from RPA predictions are a result of non-Gaussian chain statistics. Each panel in Figure 6 shows a comparison of simulation results to ROL and RPA predictions for a pair of simulations of models H and S with matched values of N̅ . The top panel shows results for simulations H64 and S16, which both have N̅ ≃ 240. The lower panel shows a comparison for model H128 and model S32, for which N̅ ≃ 480. For both values of N̅ , the results obtained with the two models agree very well with one another and with ROL predictions at wavenumbers near the peak in S(q). It is also clear in the top panel of Figure 6, however, that results from model S16 shows a systematic deviation from results of model H64 and from ROL predictions at higher values of q. These high-q deviations are results of the discrete chain structure, which become noticeable as q approaches the inverse monomer size 1/σ. This deviation from the ROL prediction is more pronounced for model S16 than for model H64 simply because H64 has 4 times as many monomers per chain. Corrections to random-walk behavior that arise from discrete chain structure can be distinguished from the universal

Figure 5. Fits of the intermolecular coordination number z(N) in a homopolymer melt (α = 0) vs N̅ −1/2, for models S (upper panel) and H (lower panel). The solid line in each graph shows the result of a fit to eq 23, using z∞ and γ as fitting parameters. The dashed line shows a truncated expansion, with ε = 0, showing the asymptotic slope.

lengths, for models S (upper panel) and H (lower). In each panel, the solid line is a fit to eq 23, in which the coefficients z∞ and δ have been adjusted to fit the data but in which the predicted coefficient of (6/π)3/2 has been used for the N̅ −1/2 term. The intercept on the left axis is the desired value of z∞. 858


Macromolecules

Article

Figure 7. Difference between measured values of S−1(q)/c in homopolymer melts and ROL predictions for this quantity vs qσ. Data are shown for four different chain lengths N = 16, 32, 64, and 128 for model S (filled symbols) and model H (open symbols). Solid curves are fits to an empirical function cSlocal−1(qσ) = a(qσ)2 + b(qσ)4. The fit parameters are a = 0.0227, b = −0.0140 (model S) and a = −0.0111, b = −0.0376 (model H).

∝ N−1/2 but can be important for very short chains. Figure 8 shows an analysis of our results for S(q) for model H16, for

Figure 6. Comparison of the simulated S(q) to one-loop theory (solid curve) and RPA (dashed curve), at α = 0, for simulations of model S and model H with matched values of N̅ ≃ 240 (upper panel, S16 and H64) and N̅ ≃ 480 (lower panel, S32 and H128).

corrections predicted by the one-loop theory because they are expected to exhibit a different dependence on q and N. Universal corrections should exhibit a universal (i.e., modelindependent) dependence on qRg0 and on N̅ . Corrections that arise from discrete chain structure should instead depend on the absolute wavenumber q, independent of Rg0 or N. If the universal and nonuniversal correction to S−1(q) are approximately additive, and if universal contributions are adequately approximated by the ROL theory, then S−1(q) should be well approximated as sum S −1(q) ≃ SROL−1(q) + S local −1(q)

Figure 8. Comparison of the simulation results for S(q) (open squares), simulation results that have been corrected by subtracting Slocal(q) (open circles) to the RPA (dashed curve) and ROL (solid line), for model H16 at α = 0. The corrected results were obtained by subtracting the N-independent function that was obtained by fitting the data shown in Figure 6 from the raw data.

(26)

−1

in which SROL (q) = FROL(qRg0,N̅ )/cN is the ROL prediction, which exhibits a universal dependence on qRg0 and N̅ , and Slocal(q) is a “local” contribution arising from discrete chain structure. The local contribution may be entirely different for different models, but should be independent of N. In Figure 7, we test eq 26 by plotting the difference c[S−1(q) − SROL−1(q)]/c vs qσ for multiple chain lengths of each model, in an attempt to isolate Slocal−1(q) . As proposed, results for this difference from different chain lengths collapse onto a common curve for each model when plotted vs qσ (rather than vs qRg0) but yield different curves for models H and S. Consistent with our interpretation, these deviations are large only for qσ ≲ 1. The solid lines through the data for each model are purely empirical fits of the function Slocal−1(q) for each model, as discussed in the figure caption. This identification of a nonuniversal “local” contribution to S−1(q) can be used to refine our comparison of simulation data for S(q) to the ROL theory. The local contribution is expected to have a negligible effect for q ∼ q* for long chains because q*

which the effects of discrete chain structure are strongest. The open squares represent the raw simulation results for S(q)/cN. The open circles represent the results for a function S′(q) obtained by subtracting our estimate of SROL−1(q) from each measured value of S−1(q). The ROL theory, indicated by the solid line, shows a somewhat astonishing level of agreement with the corrected data. Agreement is equally good for all remaining data sets for both models. Explicitly correcting for the effects of discrete chain structure thus brings the ROL into nearly perfect agreement with our simulation results in the limit α = 0, even for very short chains.

6. PEAK SCATTERING INTENSITY Some predictions of the ROL theory can be tested without estimating the interaction parameter χe. No estimate of χe is required to test predictions for symmetric copolymers with α = 0 because χe is known21 to be exactly zero in this limit. One can also test predictions for relationships among observable properties, such as the relationship between the peak 859


Macromolecules

Article

been plotted using eq 24 for χ(1) e (α). Results for model S and model H are presented in the upper and lower panels, respectively. The RPA predicts a linear function cNS−1(q*)/2 = 10.5 − χeN, which is shown by the dashed line in each plot. The horizontal dotted line in each plot shows the value of cNS−1(q*)/2 = 0.35 at which eq 16 yields ξ ≃ d*. Data from simulations very near the ODT, for which cNS−1(q*)/2 falls below this bound, have been excluded from both the plots and the fit to avoid data that may be corrupted by finite size effects. In the lower panel, data from older simulations of model H25 (open symbols) using two different box sizes for each chain length are plotted along with the results of new MD simulations (closed symbols) to show that the results are independent of L within this range of values of the ordinate cNS−1(q*)/2. When assessing the successes and failures of this approach, it is important to remember that the value of χe at a given value of χ(1) e N is inversely proportional to N, and that α is roughly proportional to χe. The range of values of α accessed in simulations of different chain lengths thus decreases by approximately a factor of 2 when N is doubled. Upon inspection of these plots, a few trends become apparent. Agreement between theory and simulation results is reasonably good for the longest chain lengths, particularly for model H. Agreement is also very good for all chain lengths, in both models, at small values of χ(1) e N, corresponding to small values of α. Agreement thus appears to be best wherever α is smallest. Conversely, agreement is worst where α is largest, at large values of χeN and small values of N. Interestingly, the nature of the discrepancies between simulation results and theory is qualitatively different for the two models. In model H (lower plot), simulations of H16 and H32 yield values of cNS−1(q*)/2 that are significantly greater than theoretical predictions. In model S (upper plot), simulation results for cNS−1(q*)/2 are instead all less than the predictions. The fact that the largest errors are consistently found where α is largest, independent of χe(1)N, suggests that these discrepancies may be primarily the result of the breakdown of our linear approximation for χe(α) at large values of α, rather than a failure of the ROL theory. This hypothesis is also suggested by the fact that the errors are very different for these two models, since the nature of nonlinearities in χe(α) may be completely different for different models. Finally, this hypothesis is consistent with the conclusion of ref 26, in which we showed that the results of three different models are consistent with the existence of underlying universal scaling relation H(q*) = H*(χe(α)N,N̅ ) in which χe(α) is a modeldependent nonlinear function. The comparison of data from different models in ref 26 did not, however, test whether the universal function H* is adequately approximated by the ROL theory. We answer this question below. 6.2. Nonlinear χe(α). In order to allow for possibility of a nonlinear dependence of χe(α) upon α, we have expressed the function χe(α) for each model as a nonlinear function with 2−3 free parameters and adjusted the values of these parameters to fit the data for four chain lengths. Despite the introduction of several adjustable parameters, this procedure provides a nontrivial test of the accuracy of the ROL theory because we attempt to simultaneously fit data from simulations of several chain lengths using a single function χe(α), without allowing χe(α) to depend on N. This is conceptually very similar to a hypothetical procedure of fitting the function χe(α) for each model so as to fit the data for one chain length, using as many

wavenumber q* and the peak intensity S(q*) or χ*a N, which was examined in ref 25. Predictions of how the peak intensity S(q*) depends upon χe cannot, however, be tested without somehow estimating how χe depends upon the control parameter α used in our simulations. The problem is completely analogous to the problem of estimating the dependence of χe(T) upon temperature T when analyzing experimental data. In this section we present comparisons of simulation results to ROL predictions for the peak intensity S(q*) using two different schemes to estimate χe(α). First, in subsection 6.1, we use the linear approximation χ(1) e (α) provided by first-order perturbation theory, as done previously in ref 25 for model H. The virtue of this approach is that it involves no free parameters, since the required coefficient of proportionality z∞ was obtained from independent simulations of homopolymer melts. This approach is found to yield reasonable agreement between theory and simulation for small values of α, and thus for long chains, but to fail at the larger values of α accessed in simulations of shorter chains. Then, in subsection 6.2, we take an approach closer to that normally used to analyze experiments, in which χ(T) must be fitted to experimental data. There, we approximate χe(α) for each model by a nonlinear function that we estimate by simultaneously fitting ROL predictions to results of simulations of several different chain lengths over the entire accessible range of values of α. 6.1. Linear χe(α). Figure 9 shows results for the normalized inverse scattering intensity cNS−1(q*)/2 vs χ(1) e (α)N. In this figure, the abscissa for each of the simulation data points has

Figure 9. Inverse peak scattering intensity cNS−1(q*)/2 vs χ(1) e N = z∞αN. Data are shown for model S (upper panel) and model H (lower panel). Solid curves are predictions by ROL theory. Data for model H corresponds to the parameters given in Table 1 (filled symbols) and to data previously published in ref 25 (open symbols). 860


Macromolecules

Article

parameters as required to fit that data, and then using the resulting function χe(α) within the ROL theory to predict scattering intensities for the other three chain lengths. Our choice of what functional form to use for χe(α) for each model is somewhat arbitrary. The only rigorous constraints are that, for symmetric diblocks, χe(α) must vanish for α = 0, for reasons of symmetry,21 and that dχe(α)/dα must approach the known value of z∞/kBT at α = 0. Different functional forms were chosen for models H and S, motivated by physical considerations. Our choice of a functional form for χe(α) in model H was motivated by the realization that, for this model, χe(α) must approach a finite limit as α → ∞. The WCA pair potential used in this model is a cutoff and shifted Lennard-Jones potential for which Vij(r) = 0 at a fixed cutoff rc = 21/6σ. In the limit εij → ∞, Vij(r) approaches a hard-sphere interaction with a hard-sphere diameter rc. In the limit α → ∞, pairs of A and B monomers thus experience a hard-core interaction with a hard-core diameter rc ≃ 1.12σ, while pairs of like (AA or BB) monomers experience a WCA potential with ε = kBT, corresponding to a slightly smaller effective hard-sphere diameter of approximately σ. In this limit, the model is thus expected to yield behavior similar to that of hard spheres with slightly unequal hard-core diameters, for which we expect a finite value of χe. Some of our simulations of model H16 use values of εAB up to 7kBT, with εAA = kBT, and thus may approach this limiting behavior. To allow for the resulting saturation, we chose to fit kBTχe(α) to a ratio of polynomials kBTχe (α) =

Figure 10. Fit of the inverse scattering intensity (cN/2)S−1(q*) using a function χe(α) for four different chain lengths N = 16, 32, 64, and 128 in model S (upper panel) and H (lower panel). The function χe(α) has been chosen as χe(α) = (z∞,Sα + aα2)/(1 + bα) for model S and as χe(α) = (z∞,Hα + cα2)/(1 + dα + eα2) for model H. The optimal fit parameters are a = 0.101, b = 0.285, c = 1.175, d = 3.515, and e = 0.750. The dashed line indicates the estimated onset of finitesize effects for S−1(q*), as discussed in the main text, and below which data points have been excluded from the fit.

z∞α + aα 2 1 + bα + cα 2

(27)

that is designed to yield a finite value kBTχe → a/c as α → ∞, while giving χe ≃ χ(1) for small α, where a, b, and c are e adjustable parameters. Simulations of model S use a much softer potential and access a range of values (α < 5 for S16) in which α = εAB − εAA remains significantly smaller than εAA = 25kBT. The function χe(α) for this model is thus expected to exhibit a nearly linear dependence on α over the relevant range of values. We have thus chosen to fit χe(α) for model S to a modified form of eq 27 in which we set c = 0 and treat a and b as adjustable parameters, giving a function that increases linearly for large α. Figure 10 shows a comparison of simulation results and ROL predictions for cNS−1(q*)/2 vs χe(α)N that are obtained by this fitting procedure. Here, values of the abscissa χe(α)N for each simulation data point have been calculated using a function χe(α) that was obtained from a simultaneous least-squares fit of data from all four chain lengths to the ROL theory. Agreement between the ROL theory and the simulation data in Figure 10 is excellent for both models. The use of a nonlinear function for χe(α) dramatically improves the agreement between simulation and theory in regions of large α, without degrading agreement in regions of small α in which the linear approximation is adequate. The extent of this improvement in two qualitatively different models strongly supports our earlier conjecture25 that the largest source of error in our comparisons of cNS−1(q*) vs χ(1) e N was the failure of the linear approximation for χe(α). Figure 11 shows the resulting fits for χe(α) for each model. Results for model S and H are shown in the upper and lower panels, respectively. Solid circles along these curves show values of α at which χe(α)N = 15 for each chain length, in order to give a sense of the range of values of α accessed for each chain

length. The resulting estimate of χe(α) for model H (lower panel) shows that simulations of models H64 and H128 accessed a range of values of α in which χe(α) is approximately linear but that values of χe(α) in simulations with N = 16 approached the limiting value.

7. PEAK WAVENUMBER Previous experiments42,44 and simulations45,47 have shown that the peak wavenumber q* in symmetric diblock copolymer melts decreases monotonically with increasing χeN, with a particularly fast drop near the ODT. The RPA instead predicts a constant peak wavenumber q*0 . Figure 12 shows a comparison of simulations results and ROL predictions for the fractional deviation (q* − q0*)/q0* from the RPA prediction plotted vs χa*N for both model S (upper panel) and model H (lower panel). In this and several subsequent plots, we plot a quantity of interest vs the apparent interaction χa*N rather than vs an estimate of χeN. Because χa*N is a measurable quantity, which is related to the inverse peak intensity by eq 7, this approach allows us to test many theoretical predictions in a manner that does not rely on any estimate of χe(α). Predicted values of q* have been calculated by numerically identifying the maximum of S(q*) predicted by the self-consistent ROL theory at each value of χa*N. Agreement between theoretical predictions and simulation is excellent for all four chain lengths for model S, though some statistically significant discrepancies are visible near the ODT, 861


Macromolecules

Article

particularly for the shortest chain, S16. Agreement is also good for models H64 and H128, which correspond to the same values of N̅ as for S16 and S32, but deteriorates noticeably for H32 and (particularly) H16. In the limit of very long chains, the ROL theory predicts a fractional shift in q* that decreases asymptotically as N̅ −1/2 with increasing N̅ . To show this, and as an introduction to the analysis of this shift given in section 8, it is useful to analyze this asymptotic behavior. The ROL prediction for q* is given by the minimum of eq 13 for H(q) with respect to normalized wavenumber Q ≡ qRg0. If we treat H0 and H1 at fixed χa*N as functions of Q alone, the peak wavenumber Q* ≡ q*Rg0 is determined by the condition 0 = H′0 (Q *) + N̅ −1/2H′1(Q *)

(28)

where prime denotes differentiation with respect to Q. Taylor expanding both functions about Q*0 = q*0 Rg0 to linear order in δQ* ≡ Q* − Q0*, using the fact that H0′ (Q0*) = 0, and retaining only terms of order N̅ −1/2 yields a dominant correction δQ * ≃ −

1 H′0 (Q 0*) N̅ 1/2 H″0 (Q 0*)

(29)

plus corrections of order N−1 and smaller. Here, H0″(Q) denotes the second derivative of H0(Q) with respect to Q. The sign of δQ* is determined by the sign of the derivative H′1(Q) of the deviation from the RPA. Figure 13 compares simulation results for the normalized deviation N̅ 1/2(q* − q*0 )/q*0 for all models except H32 and

Figure 11. Estimates of χe(α) obtained for model S (upper panel) and model H (lower panel) by fitting scattering results to the ROL theory. Solid curves show functions obtained by fitting results for the inverse peak intensity cNS−1(q*)/2, which give the fits shown in Figure 10. Solid circles on the curves indicate values of α for which χe(α)N = 15 for the indicated values of N.

Figure 13. Normalized shift N̅ 1/2(q* − q0*)/q0* in the scattering peak wavenumber vs χ*a N, scaled by N̅ 1/2, for several chain lengths in model H and S. The solid curve is the asymptotic prediction of ROL theory for infinitely long chains.

H16. The fact that the data from widely disparate chain lengths almost collapse in this representation shows that simulation results for q* − q*0 do indeed scale as N̅ −1/2 for sufficiently large enough N̅ . The solid line in this figure shows the ROL theory for this quantity in the limit N̅ → ∞, which is given eq 29. We show only the N̅ → ∞ asymptote in this figure for clarity, because the ROL predictions for this quantity do show a weak dependence on N̅ and would crowd the figure if included. The simulation results do not agree as well with this N̅ → ∞ asymptote as with predictions of the full ROL theory, which are shown in Figure 12, particularly at large values of χ*a N. A close look at the data also reveals, however, that data from the longest chains lie closest to the predicted asymptote, consistent

Figure 12. Normalized shift (q* − q0*)/q0* of the peak wavenumber q* vs χ*a N, for several chain lengths in models S (upper panel) and H (lower panel). The solid curves are predictions of ROL theory.

862


Macromolecules

Article

with the existence of slow but systematic convergence to the predicted asymptote.

8. INTRA- AND INTERMOLECULAR CONTRIBUTIONS TO S−1(q) The incompressible OZ equation provides a natural way of separating corrections to the RPA prediction for S−1(q*) that arise directly from deviations from random walk statistics, which can be quantified by measuring F(q), from more complicated collective effects, which enter via the function χa(q). As already noted, both F(q) and χa(q) can be extracted directly from simulation data. By measuring F(q) and χa(q), we gain additional insight into the physical origin of corrections to RPA predictions for S(q) and can test ROL predictions more rigorously than is possible by measuring S(q) alone. Simulation data for the quantity H(q) = cNS−1(q) may be expressed as a sum of three contributions H(q) = H0(q) + δF(q) − 2Nδχa (q)

in which H0(q) ≡ in which

cNS0−1(q)

(30)

− 2χeN is the RPA predictions and

δF(q) ≡ F(q) − F0(q) δχa (q) ≡ χa (q) − χe

(31)

Figure 14 shows the results of applying this decomposition to simulations of models S16 and H64, which have matched values of N̅ ≃ 240. Simulations for systems with α = χeN = 0 are shown in the upper panel and simulations of systems with nearly matched values of χeN ≃ 9.7 are shown in the lower panel. The quantities cNS−1(q), F(q), and χa(q) have been obtained directly from measurements of S(q) and Ωij(q), while cNS0−1(q), F0(q), and δχa(q) have been calculated using the estimates of b and χe(α) obtained in sections 4 and 6, respectively. Throughout this analysis, we use the nonlinear expressions for χe(α) obtained by fitting peak intensity data. In the lower panel, in which χeN ≃ 9.7, we show the total cNS−1(q) (triangles), the RPA contribution cNS0−1(q) (dashed line), and the corrections δF(q) (squares) and −2δχa(q)N (circles). In the upper panel, in which χeN = 0, we do not show −2δχa(q)N because it is known21 that χa(q) and δχa(q) must vanish in this limit. (We have confirmed that simulation results for χa(q) are equal to zero within statistical errors.) In each plot, we have also marked both the RPA peak wavenumber q*0 , which is the minimum of cNS0−1(q), and the actual peak wavenumber q*, which is the minimum of cNS−1(q). Plotting data from two different models (S16 and H64) with matched values N̅ and χeN allows us to confirm the universality of the results. We observe very good agreement between the two models at wave numbers near q* for all three functions δF(q), Nδχ(q), and cNS−1(q). Discrepancies between the models become noticeable at higher q, particularly in δF(q). 8.1. Physical Origin of Shift in q*. The experimental observation of a decrease in q* near the ODT was originally interpreted by Almdal, Rosedale, and Bates44 as being a result of change in the conformations of individual chains from Gaussian to “stretched”. Erukhimovich and Dobrynin10 noted, however, that the relationship between chain dimensions and q* is actually quite indirect. The physical origin of deviations from RPA predictions for q* can be clarified, however, by examining how δF(q) and −2Nδχa(q) affect the location of the minimum in cNS−1(q). Contributions to H(q) that increase

Figure 14. Individual contributions to the inverse structure factor cNS−1(q) from intra- and intermolecular correlations vs qRg0, at χeN = 0 (upper panel) and χeN ≃ 10 (lower panel). Circles represent data for the correction δF(qRg) to Leibler’s intramolecular correlation function F(qRg); squares represent data for the correction δχ(qRg) to the apparent χ-parameter χa. The intermolecular corrections have been obtained as described in the main text, using a value for χe(α)N from Figure 11. The dashed curve is the RPA prediction for cNS−1(qRg), and its minimum q0* is indicated by an arrow. The peak location q* of the actual structure factor is also indicated by an arrow. Data are shown for chains of matched value N̅ ≃ 240 in models S (N = 16, filled symbols) and H (N = 64, open symbols). The intermolecular part δχa(q) is identically zero at χeN = 0 (not shown).

with increasing wavenumber q near q*0 tend to shift q* to lower values, while contributions that decrease with increasing q tend to increase q*, as indicated in eq 29. The relevant features are visible in the two panels of Figure 14, which show typical examples of the behavior at small and large values of χeN, respectively. In both cases, δF(q) is negative and decreases with increasing q for q ∼ q*. The effect of δF(q) alone, in the absence of δχa(q), would thus be to decrease the minimum value of S−1(q*), or increase S(q*), and to increase q*. This is the behavior that is observed for χeN = 0, where both simulations and the ROL theory yield enhanced peak intensity S(q*) and a peak wavenumber q* > q*0 . For χeN > 0, the contribution −2Nδχa(a) is instead always a positive function that increases with increasing q for q ∼ q*. The effect of this contribution, by itself, would thus always be to decrease S(q*) and decrease q*. The observed behavior of q* and S(q*) is the result of a competition between these two contributions, which yields q* > q0* and enhanced values of S(q*) at very low values of χeN, where δF(q) dominates, but q* < q0* and suppressed scattering (relative to the RPA) close to the ODT, where −2Nδχa(q) dominates. The monotonic decrease of q* with increasing χeN is a direct result of the increasing importance of the contribution −2Nδχa(q) with increasing χeN. 863


Macromolecules

Article

The physical significance of this mathematical decomposition can be made clearer by considering how intra- and intermolecular interactions enter the RPA and OZ equations. The RPA expression for S−1(q) is a sum S−1(q) = F0(q)/cN− 2χe. The contribution F0/cN is a measure of the free energy cost per volume of imposing a composition fluctuation on a gas of noninteracting Gaussian copolymers and depends only on single-chain correlation functions. The contribution −2χe/c arises from corresponding change in a Flory−Huggins-like excess free energy of mixing. This term is taken to be independent of q and N because this SCF excess free energy is assumed to be a spatially local, N-independent functional of composition.18 The incompressible OZ equation has the same form but allows both for introduction of non-Gaussian chain statistics, which enter through the function F(q), and of nonlocal and N-dependent effective interactions, which enter through the apparent interaction parameter χa(q). To understand the expected physical consequences of a change in chain statistics alone, it is useful to consider the predictions of a simple modification of the RPA in which we allow for nonGaussian chain conformations by using the observed ensemble of chain conformations to calculate F(q), but still retain the local self-consistent field approximation for intermolecular interactions, and thus continue to approximate χa(q) by a constant χe. The simulation results for δF(q) show that this approximation would predict a peak wavenumber q* that is always greater than q*0 and that increases with increasing χeN because of the increase in the absolute magnitude of δF(q). In fact, q* decreases with increasing χeN and is less than q0 near the ODT. The fact that this simple modification of the RPA would actually give the wrong sign for the dependence of q* upon χeN makes it clear, we think, that attribution of the decrease in q* to changes in single chain statistics is misleading at best. Previous simulation studies have compared the fractional shift in q* to corresponding changes in the polymer radius of gyration Rg, the root-mean-squared distance RAB between the A and B blocks, and the radii of gyration of individual blocks. Fried and Binder45 were the first to show in simulations that Rg and RAB increased with increasing χeN, in qualitative agreement with the idea that the shift in q* was the result of chain stretching. These authors also pointed out, however, that the fractional shift in q* is larger than the change in either Rg or RAB, indicating that the naive association of the shift in q* with changes in overall chain dimensions is not quantitatively consistent. This comparison of q* to Rg is, moreover, arguably irrelevant. The polymer radius of gyration is related to the behavior of Ωij(q) at very low q, qRg ≪ 1. The free energy cost of a composition modulation at a nonzero wavenumber q* of order the inverse coil size is instead controlled by the value of F(q) at q*, where q*Rg ∼ 2. There is thus no inconsistency in our conclusion that Rg and RAB increase with increasing χeN but that δF(q) varies with q near q* in a manner that (in isolation) would actually tend to increase q*. 8.2. Comparison to ROL Predictions. Figure 15 shows a comparison of simulation results and ROL predictions for δF(q) and −2Nδχ(q) for models S16 and H64 at χeN ≃ 9.7. For both models, q*Rg0 ≃ 1.8. Agreement with the ROL theory is good, but not perfect, in the range of q in which the two simulation models give similar results. Figures 16 and 17 show how the values of δF(q*0 ) and δχa(q0*)N at the RPA peak wavenumber q0* depend on χa*N. In both plots, these deviations are multiplied by N̅ 1/2 in order to

Figure 15. Intra- and intermolecular corrections to the RPA for cNS−1(qRg) vs qRg0, for models S16 and H64 at χeN ≃ 10 at N̅ ≃ 240. The intramolecular correction δF(q) is indicated by circles and the intermolecular corrections −2Nδχ(qRg) by squares. Solid curves are ROL theory predictions for both quantities.

Figure 16. Normalized intramolecular part cNN̅ 1/2δF(q0) of the corrections to the RPA inverse peak scattering intensity S−1(q0) vs χa*N. Shown is data for chains of length N = 16, 32 in model S (filled symbols) and chains of length N = 64, 128 in model H (open symbols), representing values of N̅ ≃ 240 (circles) and N̅ ≃ 480 (squares). The solid curve is the prediction of ROL theory.

approximately collapse results from different chain lengths. The ROL prediction is shown by a solid line in each of these plots. If the self-consistent ROL theory were exact, this representation would exactly collapse the data from all chain lengths onto each other and onto this theoretical prediction. We have chosen to plot values of both quantities at q*0 rather than at the true peak position q* in order to simplify comparison with theory but have confirmed that almost identical plots are obtained if we plot simulation results evaluated at the true peak position q*. Figure 16 shows N̅ 1/2δF(q*0 ) vs χ*a N for all four chain lengths of model S and for models H64 and H128. The collapse of data from different chain lengths and models is quite good, confirming that the corrections are universal and that they do indeed scale as N̅ −1/2. The agreement between simulation results and ROL predictions is excellent for χeN ≃ 0 but deteriorates significantly near the ODT. Here, as in Figure 13, data from the longest chains lie closest to the ROL prediction, consistent with slow convergence toward to the ROL prediction. Because both axes in this graph are measurable quantities, it provides a test of the ROL theory that does not rely on the choice of values for any adjustable parameters. Figure 17 shows corresponding plots of N̅ 1/22Nδχa(q*0 ) vs χa*N. The upper panel shows results for model S and the lower panel for model H, showing all four chain lengths in each case. 864


Macromolecules

Article

RAB2 between the centers of mass of the two blocks. For symmetric diblock copolymers, for which Rg,AA = Rg,BB, these quantities are related by (see eqs 36−38 in ref 21) R g 2 = R g,AA 2 + RAB 2/4

(32)

It is known that Rg and RAB increase with increasing χeN but that Rg,AA shrinks. Here, we show results for Rg2 and Rg,AA2 in order to test the universality of the results from different models, the scaling with N̅ , and the level of agreement with ROL predictions. Results for the chain length dependence of Rg2 for homopolymers (or copolymers with α = 0) have already been shown in Figure 4, where we used this data to determine b. We found that the ROL accurately described the dominant 6(N̅ −1/2) corrections to random walk behavior but that we needed a nonuniversal 1/N correction to accurately fit the data for Rg2(N) vs N for short chains, particularly for model H. This 1/N correction is model dependent because it arises, in part, from the effects of discrete chain structure. For example, the difference between using N (the number of beads) and N − 1 (the number of bonds) in expressions for the mean-squared end-to-end length or squared radii of gyration would show up as a nonuniversal correction of 6(1/N ), with a prefactor that depends on details of discrete chain structure. We expect the most important effects of discrete chain structure to be the same in systems with zero or small but nonzero values of α. To focus attention on more universal phenomena, we have thus chosen to plot the difference

Figure 17. Normalized intermolecular contribution N̅ 1/2δχN(q0) to the corrections to the RPA inverse peak scattering intensity S−1(q0) vs χ*a N, for four different chain lengths N = 16, 32, 64, 128 in model S (upper panel) and model H (lower panel). To determine δχ(q0)N from χa(q0)N, the function χ(α) has been used, which is plotted in Figure 11.

ΔR g 2(N , α) = R g 2(N , α) − R g 2(N , 0)

(33)

between squared radii of gyration at nonzero and zero values of α. Figure 18 shows results for the fractional change ΔRg2/Rg02 vs χa*N for two pairs of simulations with matched values of N̅ = 240 (upper panel, S16 and H64) and N̅ = 480 (lower panel, S32 and H128). The collapse of the results from the two models in each pair confirms the universality of the results. Next, we compare the data for ΔRg2(N;α) to the ROL prediction. Because results from models S and H appear to be equivalent, we focus for this purpose on data from model S. Figure 19 shows results for the normalized deviation N̅ 1/2ΔRg2/ Rg02 vs χa*N for all four chain lengths, N = 16, 32, 64, 128, as well as the ROL prediction (solid line). If the ROL theory were exact, the data would collapse in this representation. Differences between theory and simulation data are larger for this quantity than for most others. The results for the shortest chains show qualitative differences from ROL predictions: The ROL theory predicts an essentially linear decrease with increasing χa*N for χa*N < 7 and a rapid increase near the ODT. Simulation results from the two shortest chain lengths (N = 16 and 32) show the predicted increase near the ODT but are essentially independent of χaN at lower values of χ*a N. The two longest chain lengths exhibit a linear decrease with χaN with a slope that increases with increasing N but that remains significantly smaller than predicted. Results for ΔRg2 thus do not closely approach ROL predictions over this range of chain lengths, but exhibit behavior that appears consistent with slow convergence to ROL predictions with increasing N̅ . The ROL theory is believed to yield the first correction to random walk model within an expansion in powers of N̅ −1/2. In the representation used in Figure 19, we thus expect the data to approach the ROL

Here, δχe(q*) has been calculated using the estimate for χe(α) obtained in subsection 6.2 by fitting peak intensity data. The collapse of the results for different chain length is excellent. Agreement between ROL predictions and simulation results is (again) good but not perfect. Because δχa(q*0 ) is calculated by subtracting an estimate for χe(α) from the measured χa(q*0 ), the comparison shown in this plot is sensitive to our estimate of χe(α). The agreement between ROL predictions and simulation results for H(q*) in Figure 10 is noticeably better than the agreement shown in Figures 16 and 17 for δF and −2Nδχa. This is possible because of a cancellation of errors: The simulation data for −2Nδχa(q*0 ) lie systematically below the ROL prediction in Figure 17 by an amount that is similar to the amount that data for δF(q0*) lie above the ROL prediction in Figure 16, yielding a smaller error for the sum of these two terms. (Note the difference in scales for different plots.) This is a natural consequence of the fact that the free parameters in our model for χe(α) were chosen so as to optimize agreement between ROL predictions and simulation data for H(q*). Because there are substantial errors in ROL predictions for δF(q), which are unaffected by any changes in our estimate χe(α), our procedure for fitting χe(α) introduces compensating errors in −2Nδχa(α) in order optimize the fit for H(q*).

9. CHAIN AND BLOCK RADII OF GYRATION Previous simulation studies of measured changes in the overall polymer radius of gyration Rg, the radii of gyration Rg,AA and Rg,BB of the A and B blocks, and the root-mean-square distance 865


Macromolecules

Article

convergence is clearly not valid for all of the chain lengths shown in Figure 19, since it implies that the differences between subsequent values of N should decrease by factors of 1/√2 when N is doubled, which is clearly not the case for the shorter chains. It does appear, however, that the data for the two longest chains lengths, S64 and S128, are consistent with this pattern of convergence. We have tested this as follows: We fitted a line through the data for S64 (shown in the figure) and took the difference between this interpolation and the ROL prediction to estimate the function K2. By multiplying this difference by 1/√2 and adding it back to the ROL prediction, we obtained a prediction for the expected location of data for S128. The resulting prediction, shown as a dashed line, is in reasonable agreement with our data. This confirms that data for these longest two chains is consistent with the hypothesized pattern of convergence to the ROL theory in the limit N̅ → ∞. 9.1. Rg of Individual Blocks. Figures 20 and 21 show results for the change in the radius of gyration Rg,AA of the A (or

Figure 18. Deviation ΔRg2 of the squared radius of gyration at finite interaction parameter α from the value at α = 0, normalized by the ideal radius of gyration, as a function of the χ*a N. Shown are data for model S (N = 16, 32) and model H (N = 64, 128) at two different matched values of N̅ ≃ 240 (upper panel) and N̅ ≃ 480 (lower panel).

Figure 20. Deviation ΔRg,AA2 of the squared block radius of gyration for finite α from its value at α = 0, as a function of χa*N, normalized by the ideal block radius of gyration Rg0,AA, for model H and S. Shown are data for two different matched values N̅ ≃ 240 (upper panel) and N̅ ≃ 480 (lower panel).

ΔRg2/Rg02

Figure 19. Normalized deviation N̅ of the squared radius of gyration from its value at α = 0, as a function of χ*a N, for chain lengths N = 16, 32, 64, 128 of model S. The top solid curve is a piecewise linear interpolation of the data for N = 64. The dashed curve represents the extrapolation of that curve to N = 128. The extrapolation is calculated from the difference between the top solid curve and the bottom solid curve representing ROL theory, assuming the asymptotic N̅ −1/2 approach discussed in the main text. 1/2

B) blocks of symmetric diblock copolymer, in a form analogous to that given in Figures 18 and 19 for the overall radius of gyration. As found in previous work, we find that the radius of gyration shrinks with increasing χeN. Figure 20 shows a test of the universality of this quantity, in which we show results for ΔRg,AA2/R0,AA2 from simulations of models S and H with matched values of N̅ ≃ 240 (upper panel) and N̅ = 480 (lower). In contrast to the nearly perfect universality found in a similar comparison for the overall radii of gyration, we find a noticeable difference between results from models S and H with matched N̅ . The difference is most obvious for the shortest chains (N̅ ≃ 240, upper panel), for

prediction as N̅ → ∞, with differences between the one-loop theory and data for finite N that decrease as N̅ −1/2 for sufficiently long chains. That is, we expect that N̅ 1/2ΔR g 2(α) R g0 2

≃ K1(χa* N ) + N̅ −1/2K 2(χa* N ) + ...

(34)

where K1(χa*N) is the ROL prediction and K2 is a correction arising from two-loop corrections.17,20 This pattern of 866


Macromolecules

Article

this is a real physical effect, rather than a finite size artifact, since it is an obvious, very reproducible feature that appears in systems of very short chains (e.g., S16) at somewhat lower values of χa*N for which we have seen no evidence of finite size effects. The stretching of individual blocks, as well as the distance between blocks, is a generic feature of strongly segregated ordered phases in which comparatively sharp interfaces separate regions of nearly pure A or B. In this limit, individual blocks must stretch as the unit cell size increases in order to allow polymers that are tethered to interfaces at one end to uniformly fill domains of nearly pure A or B, as required by the incompressibility constraint. This is very different from the behavior found in the disordered phase farther from the ODT, which is predicted by the ROL theory, in which sufficiently long individual blocks respond very much like polymers in a poor solvent and shrink in order to increase their self-concentration and decrease contact with monomers of the other species. We thus conjecture that this upturn reflects a crossover to a structure of more strongly segregated disordered domains in which the relationship between the block size and domain spacing begins to more closely mimic that expected in any strongly segregated structure. Whatever the origin of this phenomenon, it is missed entirely by the ROL theory. Erukhimovich and Dobrynin10 and Barrat and Fredrickson7 have previously presented predictions for changes in the dimensions of polymers and individual blocks within a block copolymer melt. The calculation of Barrat and Fredrickson predicts an expansion of the radius of gyration of the chain as a whole and shrinkage of individual blocks, in qualitative agreement with those obtained from the ROL theory. This calculation is, however, based on an expression for the screened interaction that vanishes in the limit χ = 0 (their eq 7 for Φ), which we believe to be incorrect. The expression used in the ROL theory (eq 77 of ref 18) is the same as that obtained by Vilgis and Borsali,38 which reduces in the limit χ = 0 to the appropriate nonzero screened interaction for an incompressible homopolymer melt. It appears to us that Barrat and Fredrickson may have made an error between their eqs 6 and 7 when taking the incompressible limit of a product of matrices that become singular or divergent in this limit. The calculation of Erukhimovich and Dobrynin for both symmetric and asymmetric diblock copolymers is based on a starting point closely analogous to that of the ROL theory for single-chain properties, including the use of the expression of Vilgis and Borsali for the screened interaction (eqs 6 and 7 of ref 10). The Erukhimovich−Dobrynin calculation, like other calculations of that era, relies, however, on several mathematical approximations that limit its intended region of validity to the immediate vicinity of the ODT. The most obvious of these is the use of a Lorentzian approximation for the integrand in all Fourier integrals that become sharply peaked around q* near the ODT. Quantitative comparison of our predictions to those of ref 10 is complicated by the fact that these authors calculated corrections to the mean-squared end-to-end length of polymer and individual blocks, whereas we have calculated corresponding radii of gyration, and by the fact that these authors did not explain some aspects of the fully self-consistent version of their theory. We note, however, that the results plotted for the more straightforward perturbative version of their theory, in which the screened interaction is calculated using the RPA χ parameter (Figure 1, plot b, dashed line 1), show an increase in the size of individual blocks upon increasing χN near the

Figure 21. Normalized deviation N̅ 1/2ΔRg,AA2/Rg0,AA2 of the squared block radius of gyration Rg,AA at finite α from its value at α = 0 vs χ*a N for four different chain lengths N = 16, 32, 64, and 128 in model S. The top solid and the dashed curve represent a piecewise linear interpolation for N = 64 and an extrapolation for N = 128, respectively, as discussed in the main text. The bottom solid curve is the asymptotic prediction of ROL theory.

which the discrepancy is about 1% at χa*N ≃ 8. This difference is much less for the longer chains with N̅ ≃ 480. We suspect that this difference may be the result of a slight α-dependence of the statistical segment length b, which (if present) would be different for different models, but which we have ignored thus far. Because the difference is quite small, it appears to us that a discrepancy of this magnitude would have been difficult to detect in our earlier tests of universality of data for q* and Rg. As an example, we note that in the upper panel of Figure 18 values of Rg2 for model S (open circles) actually are slightly below those for model H (filled circles) and that a 1% increase in the value of Rg2 for model S and χ*a N ≃ 8 would not substantially change the apparent level of agreement in that figure. Despite this evidence of a slight nonuniversality, for the sake of simplicity, we have chosen not to refine the analysis given in this paper so as to allow for an α-dependent statistical segment length for each model. Figure 21 compares simulation results for the normalized deviation N̅ 1/2ΔRg,AA2/Rg0,AA for all four chain lengths of model S to ROL predictions for this quantity. The relationship between the simulations results and the ROL predictions is very similar to that obtained with the overall radius of gyration: There are substantial quantitative discrepancies over the range of chain lengths studied here, but the behavior is consistent with convergence to ROL predictions in the limit N̅ → ∞. To show this, we have applied the same analysis of the relationship between data from S64 and S128 as that applied to the overall radius of gyration. The solid line through the data from model S64 is an interpolation, and the dashed line is a prediction for where we expect results for S128 to lie based on the results for S64. The relationship of the data from these two longest chain lengths again seems to be consistent with our hypothesis of slow convergence toward the ROL prediction. Simulation results for the block radii of gyration also, however, show one qualitative feature that is not described correctly by the ROL theory: All of the results for ΔRg,AA2 show an upturn at large values of χa*N, very near the ODT, while the ROL theory predicts a monotonic decrease. The position of the resulting minimum in ΔRg,AA2 vs χ*a N moves to closer values of χ*a N with increasing N̅ . Much, but not all, of the data that show this feature are in the region 10.5 − χa*N ≲ 0.35 where finite size effects are a potential problem. Despite this, we think that 867


Macromolecules

Article

predicted by the ROL theory. We thus conclude the observed shift in q* is not a direct result of changes in single-chain statistics or overall coil size. Predictions for a variety of single-chain properties, including radii of gyration of chains and blocks and F(q0*), show behavior that appears to be consistent with convergence toward ROL predictions as N̅ → ∞, but with slower convergence than seen for other quantities. As reported in earlier simulations, the radii of gyration of individual blocks shrink with increasing χeN over most of the range of parameters that we studied. Interestingly, however, we also see evidence of an increase of the block radii of gyration with increasing χeN very close to the ODT. Because this behavior occurs in a region where finite size effects are a potential problem, further study using larger systems would be useful, but we suspect that this qualitative observation is robust. We conjecture that this represents a crossover to a strongly segregated disordered regime in which the relationship between the block radii and the domain size becomes similar to that found in ordered phases. This is the only phenomena observed in our simulations that was not reasonably well described by the ROL theory, and suggests a possible breakdown of the theory very near the transition in systems of relatively short chains. The ability of the ROL theory to accurately predict correlations in the disordered phase over a wide range of parameters and chain lengths is potentially useful in several contexts. Most importantly, it is potentially useful as a tool for analyzing experimental data, which we hope will allow a more consistent description of small-angle scattering from block copolymer melts and polymer mixtures. Further work along these lines should pursue careful comparisons to experimental data. The theory is also potentially useful as a tool for assigning interaction parameters in simulations of ordered phases of block copolymers, by fitting measurements of scattering in the disordered phase, to allow more meaningful quantitative comparisons of results obtained with different models. Finally, the ability of the theory to predict the behavior of rather short chains, for which the RPA is inadequate, makes it potentially useful as a tool for interpreting the chain-length dependence of results obtained from more atomistic simulations of mixtures of short chains.

RPA spinodal, in contrast to the monotonic shrinkage predicted by the analogous ROL calculation. We do not see any obvious reason for this unexpected qualitative difference but have not examined this earlier calculation in detail.

10. CONCLUSIONS We have compared the results of extensive simulations of disordered symmetric diblock copolymers to predictions of the renormalized one-loop theory. Comparison of results from two rather different models, which were designed to yield pairs of systems with matched values of N̅ , allowed us to test the universality of all results. The use of two different models and chains of up to N = 128 beads also allowed us to explore a range of values of N̅ = 60−1920 that overlaps with the experimentally relevant range. The ROL theory is found to yield accurate predictions for all quantities that we have studied for sufficiently long chains and remarkably accurate predictions for some. All of our simulation results appear to be consistent with the hypothesis that the ROL correctly predicts the dominant 6(N̅ −1/2) corrections to the RPA in the limit N̅ → ∞, though the rate of convergence is different for different quantities. The ROL theory gives particularly accurate predictions for the structure factor S(q) of symmetric copolymers in the limit α = 0. In this limit, the theory is found to remain accurate even for very short chains (e.g., model H16), if we explicitly correct for contributions that arise from local chain structure. Because S(q) is entirely determined by intramolecular correlations in this limit, this success of the theory is a testimony to the accuracy of the theory of universal corrections to the random walk model of single-chain correlations in dense onecomponent liquids that was pioneered by Wittmer, Semenov, Beckrich, and collaborators.14−16 The theory also provides very accurate predictions for the peak wavenumber q* for all but the shortest chains. Because we use the measured “apparent” peak intensity χa*N as a correlating variable for this purpose, our analysis of the behavior of q* does not involve any fitting parameters. To test ROL predictions for the dependence of the peak intensity S(q*) upon the RPA interaction parameter χe, we needed to estimate the function χe(α). A linear approximation for χe(α), which was used for model H in an earlier report,21 was found to be adequate for the longest chains studied here, but inadequate for the range of α used in simulations of short chains. It has been shown here, however, that the ROL theory can be brought into excellent agreement with data from both models, even for short chains, by allowing for the possibility of a nonlinear function χe(α). Model H was shown to yield a strongly nonlinear function χe(α) that exhibits saturation at large values of α, for reasons that are well understood. The physical origins of corrections to RPA predictions for q* and S(q*) have been clarified by separating intra- and intermolecular contributions to cNS−1(q). We discussed a simple modification of the RPA in which changes in singlechain correlations are taken into account by using the measured single-chain correlation functions to calculate Leibler’s function F(q), rather than random-walk statistics, but in which we retain the RPA treatment of χ as a wavenumber-independent parameter. It was shown that the resulting approximation would actually predict an increase in q*. The observed decrease in q* with increasing χeN is instead a result of the wavenumber dependence of the apparent interaction parameter χa(q)

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected] (D.C.M.). Present Addresses

J.G.: Department of Chemical Engineering, University of Michigan, 2300 Hayward St., Ann Arbor, MI 48109. J.Q.: Institute for Molecular Engineering, University of Chicago, 5747 S. Ellis Ave., Chicago, IL 60637. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This work was supported by National Science Foundation (NSF) Award DMR-0907338 and by a postdoctoral fellowship for J.G. from the Deutsche Forschungsgemeinschaft (DFG), GL 733/1-1. This research used resources of Minnesota Supercomputing Institute, and of the Keeneland Computing Facility at the Georgia Institute of Technology, which is supported by the NSF under Contract OCI-0910735. 868


Macromolecules

■

Article

REFERENCES

(1) Matsen, M. W. J. Phys.: Condens. Matter 2002, 14, R21−R47. (2) Leibler, L. Macromolecules 1980, 13, 1602−1617. (3) Fredrickson, G. H.; Helfand, E. J. Chem. Phys. 1987, 87, 697. (4) Brazovskiĭ, S. Sov. Phys. JETP 1975, 41, 85. (5) de la Cruz, M. O. Phys. Rev. Lett. 1991, 67, 85−88. (6) Mayes, A. M.; de la Cruz, M. O. J. Chem. Phys. 1991, 95, 4670. (7) Barrat, J.-l.; Fredrickson, G. H. J. Chem. Phys. 1991, 95, 1281− 1289. (8) Dobrynin, A. V.; Erukhimovich, I. Y. J. Phys. II 1991, 1, 1387− 1404. (9) Dobrynin, A. V.; Erukhimovich, I. Y. Polym. Sci. 1991, 33, 1012− 1027. (10) Erukhimovich, I. Y.; Dobrynin, A. V. Macromolecules 1992, 25, 4411. (11) Holyst, R.; Vilgis, A. J. Chem. Phys. 1993, 99, 4835−4844. (12) Kudlay, A.; Stepanow, S. J. Chem. Phys. 2003, 118, 4272−4276. (13) Wang, Z.-G. J. Chem. Phys. 2002, 117, 481. (14) Wittmer, J. P.; Meyer, H.; Baschnagel, J.; Johner, A.; Obukhov, S.; Mattioni, L.; Müller, M.; Semenov, A. N. Phys. Rev. Lett. 2004, 93, 147801. (15) Wittmer, J. P.; Beckrich, P.; Meyer, H.; Cavallo, A.; Johner, A.; Baschnager, J. Phys. Rev. E 2007, 76, 011803. (16) Beckrich, P.; Johner, A.; Semenov, A. N.; Obukhov, S. P.; Benot, H.; Wittmer, J. P. Macromolecules 2007, 40, 3805−3814. (17) Morse, D. Ann. Phys. 2006, 321, 2318−2389. (18) Grzywacz, P.; Qin, J.; Morse, D. Phys. Rev. E 2007, 76, 061802. (19) Qin, J.; Morse, D. C. J. Chem. Phys. 2009, 130, 224902. (20) Morse, D. C.; Qin, J. J. Chem. Phys. 2011, 134, 084902. (21) Qin, J.; Grzywacz, P.; Morse, D. C. J. Chem. Phys. 2011, 135, 84902. (22) Kremer, K.; Grest, G. S. J. Chem. Phys. 1990, 92, 5057. (23) Grest, G. S.; Lacasse, M.-D.; Kremer, K.; Gupta, A. M. J. Chem. Phys. 1996, 105, 10583. (24) Groot, R. D.; Warren, P. B. J. Chem. Phys. 1997, 107, 4423. (25) Qin, J.; Morse, D. C. Phys. Rev. Lett. 2012, 108, 238301. (26) Glaser, J.; Qin, J.; Medapuram, P.; Mueller, M.; Morse, D. Soft Matter 2012, 8, 11310−11317. (27) Erukhimovich, I. Y. Polym. Sci. U. S. S. R. 1979, 21, 470−476. (28) Erukhimovich, I. Y. Polym. Sci. U. S. S. R. 1982, 24, 2223−2232. (29) Benoit, H.; Benmouna, M. Polymer 1984, 25, 1059−1067. (30) Schweizer, K.; Curro, J. Adv. Polym. Sci. 1994, 116, 321−377. (31) Schweizer, K. S.; Curro, J. G. Adv. Chem. Phys. 1997, 98, 1−142. (32) Ganesan, V.; Fredrickson, G. Europhys. Lett. 2001, 55, 814. (33) Fredrickson, G.; Ganesan, V.; Drolet, F. Macromolecules 2002, 35, 16−39. (34) Lennon, E.; Katsov, K.; Fredrickson, G. Phys. Rev. Lett. 2008, 101, 138302. (35) Edwards, S. F. Proc. Phys. Soc. 1966, 88, 265−280. (36) Doi, M.; Edwards, S. F. The Theory of Polymer Dynamics; Oxford University Press: Oxford, 1988. (37) Brereton, M. G.; Vilgis, T. A. J. Phys. (Paris) 1989, 50, 245−253. (38) Vilgis, T. A.; Borsali, R. Macromolecules 1990, 23, 3172−3178. (39) Anderson, J.; Lorenz, C.; Travesset, A. J. Comput. Phys. 2008, 227, 5342−5359. (40) Morse, D. C.; Chung, J. K. J. Chem. Phys. 2009, 130, 224901. (41) Mueller, M. Macromolecules 1995, 28, 6556−6564. (42) Bates, F. S.; Rosedale, J. H.; Fredrickson, G. H. Phys. Rev. Lett. 1988, 61, 2229−2232. (43) Bates, F. S.; Rosedale, J. H.; Fredrickson, G. H. J. Chem. Phys. 1990, 92, 6255−6270. (44) Almdal, K.; Rosedale, J. H.; Bates, F. S.; Wignall, G. D.; Fredrickson, G. H. Phys. Rev. Lett. 1990, 65, 1112−1115. (45) Fried, H.; Binder, K. J. Chem. Phys. 1991, 94, 8349−8366. (46) Fried, H.; Binder, K. Europhys. Lett. 1991, 16, 237−242. (47) Vassiliev, O. N.; Matsen, M. W. J. Chem. Phys. 2003, 118, 7700− 7713.

869


Collective and Single-Chain Correlations in Disordered Melts of

Recommend Documents