Ordering Transition in Salt-Doped Diblock ... - ACS Publications

Apr 26, 2016 - the so-called “order−disorder transition temperature” of .... the nonbonded interactions, bead positions must be converted .... n...
0 downloads 0 Views 1MB Size
Download PDF
Article pubs.acs.org/Macromolecules

Ordering Transition in Salt-Doped Diblock Copolymers Jian Qin†,‡ and Juan J. de Pablo*,†,‡ †

Institute for Molecular Engineering, University of Chicago, Chicago, Illinois 60637, United States Argonne National Laboratory, Argonne, Illinois 70439, United States

‡

S Supporting Information *

ABSTRACT: Lithium salt-doped block copolymers offer promise for applications as solid electrolytes in lithium ion batteries. Control of the conductivity and mechanical properties of these materials for membrane applications relies critically on the ability to predict and manipulate their microphase separation temperature. Past attempts to predict the so-called “order−disorder transition temperature” of copolymer electrolytes have relied on approximate treatments of electrostatic interactions. In this work, we introduce a coarse-grained simulation model that treats Coulomb interactions explicitly, and we use it to investigate the ordering transition of charged block copolymers. The order−disorder transition temperature is determined from the ordering free energy, which we calculate with a high level of precision using a density-of-states approach. Our calculations allow us to discern a delicate competition between two physical effects: ion association, which raises the transition temperature, and solvent dilution, which lowers the transition temperature. In the intermediate salt concentration regime, our results predict that the order−disorder transition temperature increases with salt content, in agreement with available experimental data.

■

INTRODUCTION Conventional lithium ion batteries suffer from safety concerns related to the high reactivity of liquid electrolytes with lithium ions.1 Solid polymer electrolyte materials that conduct ions and exhibit sufficient mechanical strength offer an interesting alternative to liquids.2−4 Several material design strategies have been proposed to improve lithium ion conductivity while offering mechanical strength; among these, lithium salt-doped block copolymers offer considerable promise. 2,3 These materials take advantage of the fact that microphase separation may occur under favorable conditions (typically triggered by temperature), leading to formation of periodically modulated microdomains. Both the characteristic domain size and the type of morphology can be tailored by suitable material choices. In this regard, PEO−PS diblock copolymers doped with lithium salts (e.g., LiTFSI) have been studied extensively. The PEO block is able to dissolve lithium salts,1 and the PS block is used because it exhibits a relatively high glass transition temperature (Tg ≃ 100 °C). When microphase separation occurs, at room temperature, the ion-associating PEO-rich domains conduct ions and the PS-rich domains confer mechanical support. Conductivity measurements have demonstrated the materials’ potential for ion-conducting membranes, and significant efforts have been devoted to optimize device performance by varying the types of anions and the underlying phase morphologies.3 An important parameter for this class of materials is the order−disorder transition temperature (TODT) at which microphase separation occurs. Experimentally, it has been found that TODT is extremely sensitive to the amount of added © XXXX American Chemical Society

salt in the system. Addition of a few percent of salt raises TODT by tens of degrees.5,6 Such dramatic effects have also been examined theoretically.4,7,8 In particular, a combination of selfconsistent field theory with a Poisson−Boltzmann treatment for electrostatic interactions has been able to explain a number of the experimentally observed features seen in salt-doped diblock copolymers, including the variation of TODT and microdomain spacing with addition of salt.4 More recently, a different approach has relied on a combination of liquid-state theory and self-consistent field theory to examine the phase behavior of charged copolymers and polymer blends, leading to intriguing predictions on the effects of charge density and electrolyte concentration on phase behavior.9−11 Given the mean field nature of past studies, however, it is of interest to approach the problem from a different perspective and, importantly, to introduce fluctuations by relying on molecular simulations. Ion correlations are known to give rise to phenomena in polymeric materials, such as attractions between chains of the same charge,9,10 that cannot be explained in their absence. While several simulations of ion-containing block copolymers have been reported in the recent literature,12−14 these have not considered Coulombic interactions explicitly. Furthermore, past simulations have not had the precision required to clearly delineate competing contributions to the order−disorder Received: December 5, 2015 Revised: April 10, 2016

A

DOI: 10.1021/acs.macromol.5b02643 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules

∑kδk,iδ(r − ri) and normalize bead volume fractions by ϕi(r) ≡ ρi(r)/ρ, with k running over all particles. In the limit of a homogeneous liquid without composition fluctuations, the nonbonded interaction becomes ρVχijϕ̅ iϕ̅ j, where ϕ̅ i is the number fraction for beads of type i. The components χij vanish for i = j and are symmetric with respect to exchange of i and j. Since the system contains four species, six independent parameters are required to fully specify the incompatibility. We note, however, that the key physical ingredients are the incompatibility between species A and B and the preferential association of the PEO-like block (A) with the cations (C). So to arrive at a minimal model, we set the values of χij to 0 for AD, BC, BD, and CD pairs and assign a positive and negative value for χAB and χAC, respectively. The term χAB is allowed to vary in our simulations, and it represents the effects of block incompatibility. The value of χAC needs to be a large negative number to mimic the strong association between the PEO block and the cations. The choice of χAC is discussed in the Supporting Information. To evaluate the nonbonded interactions, we discretize the simulation box into a grid of N3g sites, and regularize the δ function by mapping the particles onto the nearby lattice points

transition. Our purpose in this work is to introduce a coarsegrained simulation model that treats electrostatic interactions explicitly and use it to investigate the effects of salt addition and fluctuations on TODT. We do so through a powerful new simulation methodology that enables direct sampling of the free energy by relying on a density-of-states based approach. With this algorithm, we can distinguish a “dilution” effect, created by addition of ions, and an ion-association effect, which arises from Coulombic interactions. Depending on the conditions of interest, these two effects can compete with each other, leading to increases or decreases of the ODT. The next two sections present the details of the simulation model and the methodology adopted here to calculate the free energy. Our results on morphological characterization and variation of the ODT are then discussed, followed by a brief summary of our main findings.

■

SIMULATION MODEL The model considered here is based on an extension to a theoretically informed coarse-grained model that has been used extensively for studies of neutral polymeric systems.15,16 We consider a system of diblock copolymers doped with fully dissolved salts and label the two blocks, the cations, and the anions by A, B, C, and D, respectively. Conceptually, we may view the species A, B, C, and D as the PEO block, the PS block, the lithium ion, and the anion, respectively. The size of these four types of beads is the same in our model, which results in a highly coarse-grained representation that is conceived to capture the essential features of the system but that is not designed to represent detailed atomistic behaviors. Polymers are modeled as neutral beads connected by springs, and cations and anions are modeled as beads of opposite charge. The number of polymer chains, the number of beads per chain, and the number of cations and anions are denoted by M, N, Nc, and Na, respectively. Throughout this work, we set Na = Nc to ensure charge neutrality. The model Hamiltonian includes the polymer stretching energy, bead nonbonded interactions, and Coulomb interactions between charged species: H = Hb + Hnb + Hc. The stretching term Hb represents the conformational entropy associated with polymer bead connectivity and can be written 2 2 as Hb/(kBT) = 3/(2b2)∑N−1 i=1 (Ri+1 − Ri) , where 3/b is the effective spring constant and b is analogous to the chain statistical segment length. For simplicity, we use a symmetric model and assume identical b values for A and B blocks. Throughout this work, we use the polymer end-to-end distance Re as the unit of length, so the statistical segmental length b is given by 1/(N − 1)1/2. The nonbonded interaction energy Hnb accounts for the excluded volume interaction and the bead−bead incompatibility. We treat it as pairwise additive, in the spirit of the Flory− Huggins model, and express it as an integration over the bead density fields: Hnb/(kBT) = ρκ/2∫ dr (∑iϕi(r) − 1)2 + ρ/2∫ dr ∑i,jχijϕi(r)ϕj(r). Here ρ ≡ (NM + Nc + Na)/V is the bead number density, and V is the system volume. The spatially varying bead volume fractions are denoted by ϕi(r), where i = A, B, C, or D, representing the four types of beads. The parameter κ plays the role of a bulk modulus, which suppresses the overall density fluctuations. The parameters χij quantify the strength of bead repulsion and are physically related to the Flory−Huggins parameters between species i and j. To evaluate the nonbonded interactions, bead positions must be converted into density fields. We define the bead densities by ρi(r) =

Figure 1. Schematic illustration of theoretically informed coarsegrained model for salt-doped diblock copolymers. The density fields as required by the Hamiltonian are constructed by mapping bead positions onto an underlying grid. The mass and charge of beads are mapped onto the eight nearest-neighbor grid sites (a, b, c, d, etc.) and the contributions of different sites to the free energy are determined by the relative distance of each bead to a nearby grid site.

linearly, as illustrated in Figure 1.15,16 The field integration is then converted to a lattice summation, of the form Hnb = χAB ρv ∑

nA, i nB, i

i

=

N ρv

∑ χAB N i

ρv ρv nA, i nB, i N N

+ +

⎞2 ⎛ nA, i nB, i κ ρv ∑ ⎜ + − 1⎟ 2 ρv ⎠ i ⎝ ρv N ρv

nB, i ⎛ nA, i ρv ⎞ 2 ⎟ + − 2 N N N⎠

∑ κN ⎜⎝ i

Note that we have used a pure AB diblock copolymer to explain the discretization detail and generalizing the model to saltdoped systems by including additional density fields ϕC and ϕD is straightforward. Here v ≡ V/Ng3 is the grid site volume, and nA,i and nB,i are numbers of beads of type A and B on the ith site. On average, nA,i and nB,i are proportional to N, and ρ is proportional to N by definition. So the dependence of Hnb on N is weak for fixed values of χABN and κN. In simulations, this model is fully specified by the values of the interaction strength χABN, the incompressibility κN, the number of molecules M, and the grid site length v1/3. The magnitude of Hnb is proportional to N/(ρv), or Ng3/M, the inverse of the number of chains per site. By using Re as the length unit, the term N/(ρv) can be written (N/ρRe3)(Re3/v) = B

DOI: 10.1021/acs.macromol.5b02643 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules 1/[N̅ 1/2(v/Re3)], where N̅ ≡ (ρRe3/N)2 measures the degree of chain overlap and is used in our simulations to control the number of chains in the simulation box. In practice, for a fixed value of N, if one needs to vary the value of N̅ systematically, such that the density of nonbonded interactions is a constant, the proper way is to (1) decrease the value of N̅ and (2) increase the site volume v/Re3 accordingly, to ensure that the number of beads per site remains constant. Essentially this corresponds to reducing the N̅ value at fixed N/(ρv) (which is set to 1.6875), for parameters N = 32, M = 128(V/Re3), and v = (Re/6)3, the typical value used in our previous simulations.15,16 The above order-of-magnitude analysis applies to pairwise nonbonded interactions, and can be generalized to the nth order interactions, such as the cubic interactions needed by describing the process of solvent annealing.17 A generic nth order term may be converted into a summation as follows: ⎛n ⎞ ε N ⎛R ⎞ dr ϕ(r) = εnρv ∑ ⎜ i ⎟ = (nn− 1)/2 ⎜ e ⎟ ⎝ v ⎠ ⎝ ρv ⎠ N̅ i n

εnρ

∫

n

n−1

⎛n ⎞ ∑ ⎝⎜ i ⎠⎟ N

grid sites in all periodic cell copies. The prime sign excludes terms with i = j. The above expression essentially provides the Coulomb energy of 3D periodic grid charges. This energy can now be evaluated through Green’s function appropriate for Poisson’s equation with 3D periodic sources.19 The result is given by Hc =

1 1 4π ϵ 2

Na + Nc

∑ ∑′ i,j=1

n

G(r) =

e2 1 4π ϵ 2

∫ dr dr′ ρ|(rr)−ρ(rr′|′)

8 × πLxLyLz

∑ k ,l ,m=0

i

γklm

cos(2kπrx /Lx) cos(2lπry /Ly) cos(2mπrz /Lz) (k /Lx)2 + (l /Ly)2 + (m /Lz)2

where the values of symmetry coefficients are such that γ000 = 0, γk00 = γ0k0 = γ00k = 1/4 for k > 0, γ0kl = γk0l = γkl0 = 1/2 for k,l > 0, and γklm = 1 for k,l,m > 0. Note that the maximum wavenumber is Ng/2, and G(r) is invariant with respect to shifting r by integer multiples of lattice vectors, i.e., G(r) = G(r + L). The expression provided by eq 3 only describes the Coulomb interaction among different sites. To correct for the selfinteraction of charges belonging to the same site, one may incorporate an additional single-site term by assuming that the charges are uniformly distributed within the grid cell and evaluate the self-energy for such a cloud. The result is given by Hc,self = (e2/8πϵ)(Cu/δL)∑icellQi2, where δL is the site dimension and Cu is a constant that depends on the site shape, which equals 1.882 312 645 for a cubic site.18,20 Alternatively, a Gaussian smearing function may be used to regularize the Coulomb energy.21,22 By calculating the Coulomb interaction on a lattice, a short wavelength cutoff v1/3 is introduced; the model is therefore only expected to capture the long-wavelength behavior. The computation cost is high even though the Green’s function G(r) can be tabulated prior to simulations. The cost of the Coulomb interaction calculation scales with Ng6, since eq 3 involves a double summation over Ng3 lattice sites. This is in sharp contrast to the computation cost of the nonbonded interactions, which scales with Ng3. However, since we are relying on Monte Carlo simulations, whenever the bead position is updated, one only needs to keep track of the local change in the bead density or charge density, and the evaluations of nonbonded and Coulomb interactions for one bead move scale with Ng0 and Ng3, respectively. To facilitate the sampling of configuration space, a variety of Monte Carlo moves are implemented, including bead random displacements, chain reptation move, and chain head-to-tail flip moves for the copolymers. Treating Coulombic interactions by using eq 1 assumes that the medium’s dielectric permittivity does not vary with position. This prevents us from examining ordered phases having well-defined domains, with spatially varying dielectric profiles, and from considering the effects of Born solvation, which matters only for dielectrically heterogeneous media. Given this limitation, we focus on the disordered phases, and for the ordered phase, we focus on the weak segregation regime.

(1)

where ϵ is the mean dielectric constant, qi is the charge carried by the ith particle, +e for cations and −e for anions (we focus on monovalent ions throughout), and rij ≡ ri − rj. The n summation runs over all periodic cells, and L = {Lx, Ly, Lz} is the vector representing the simulation box dimensions. The prime sign in the second summation excludes the cases i = j and n = 0. The expression in eq 1 includes contributions from the simulation cell and its all periodic images, whose magnitude is controlled by the Bjerrum length, lB ≡ e2/(4πϵkBT), the length scale at which the Coulomb interaction becomes comparable to kBT. Calculation of the Coulomb interaction energy is expensive because it involves an infinite summation. To circumvent this problem, and in analogy to our treatment of the nonbonded interactions, we convert the Coulomb energy calculation into a field integration.18 We first define a charge density field by ρ(r) = ∑iqiδ(ri − r), where qi is the particle charge and the summation runs over all particles in the system. The charge density ρ(r) is periodic a simulation cell with the periodic boundary condition. Then the Coulomb interaction eq 1 can be written as Hc =

(3)

i,j

Ng /2

n

qiqj |rij + n ·L|

cell

∑ ′Q iQ jG(ri − rj)

where the summation is with respect to lattice sites within one cell, and G(r) is the Green’s function defined by

where εn is the parameter controlling the interaction strength, ϕ(r) is the volume fraction per site of a generic species, and ni is the bead number per site. The total contribution is also controlled by N̅ −(n − 1)/2 . The apparent divergence with v is compensated by the vanishingly small ni in the same limit. We now discuss how Coulombic interactions are built into the model. In its exact form, the Coulombic interactions among the Na + Nc charged particles in our system can be written as Hc =

e2 8π ϵ

(2)

where the domain of integration is the entire 3D space. With the lattice used here to map out particle number density, one can also map particle charges onto the grid using the scheme of Figure 1, which leads to a 3D periodically charged grid, with Ng3 grid sites and with Qie charges per site. Then the Coulomb interaction is given by Hc = (e2/8πϵ)∑i,j′(QiQj/|ri − rj|), where ri and rj are lattice site vectors. The indices i and j run over all C

DOI: 10.1021/acs.macromol.5b02643 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules

Figure 2. Free energies as a function of the multimode order parameter at different values of χABN for a pure diblock copolymer. Parameters: N̅ = 256, N = 8. The minimum near ψ = 0.15 corresponds to the disordered phase, and the minimum near ψ = 0.25 corresponds to the lamellar phase. The extrapolated value of (χABN)ODT is 18.3.

however, the system is brought increasingly close to the order− disorder transition, and the free energy cost of certain composition modulation modes decreases progressively until, in the end, the free energy cost vanishes at the ODT and becomes negative beyond the ODT. This suggests that to calculate the free energy, we may choose an order parameter based on the composition fluctuation modes of A and B monomers, ψ(q) = (1/ NM)∑km(k)eiq·rk, where q is the wave vector corresponding to the mode to be sampled, index k runs over all beads in the system, and rk is the position vector of the kth bead.30 The mode selection function m(k) equals 1 and −1 for A- and Btype beads, respectively, and equals 0 for C- and D-type beads. This definition of ψ(q) selects the composition fluctuation mode associated with the segregation of A and B monomers on the length scale of 2π/|q|. A nearly vanishing ψ(q) value represents the disordered phase, and a finite ψ(q) value represents the lamellar phase. In this work, we always focus on the q values with wavelength close to the natural periodicity of the lamellar morphology q*. In practice, three subtle effects need to be considered to improve the order parameter identified above. First, this definition neglects the fact that the lamellar phase might be displaced along the perpendicular direction, which would contribute a phase factor to the exponential term in the definition. Physically, composition modes differing just by a phase factor should be treated on the same grounds. To take this into account, we only use the amplitude |ψ| as the order parameter. Second, a lamellar morphology with the same periodicity may be oriented along different directions, and that orientational degree of freedom should be included in the definition of free energy. So a true order parameter must contain an average over contributions from q vectors with the same wavelength, but with different orientations. Third, as suggested by the Brazovskii ̆ model,24−28 the free energy of a lamellar ordering transition must also include contributions from those wave vectors surrounding the most favorable q* value. These effects therefore suggest that the order parameter N may be defined as ψ = (∑i=1q |ψ(qi)|4/Nq)1/4, where the averaging is carried out with respect to many wave vectors q around |q| = q*. The exponents 4 and 1/4 are purely empirical and are designed to cleanly delineate the minima corresponding to the disordered and lamellar phases, respectively. Additional details have been discussed in the context of metadynamics simulation.30 To calculate the free energy by using the multimode order parameter identified above, we resort to the expandedensemble density-of-states (EXE-DOS) sampling approach,31 which is implemented by using the Wang−Landau algorithm.32−34 In this approach, a bias potential that depends on

Our model differs from that of refs 4 and 8 in several aspects. First, we neglect the spatial heterogeneity in dielectric response, which prevents us from investigating the effects of Born solvation. Second, we do not treat lithium ions as bounded to PEO-like blocks, as was done in refs 4 and 8; our model therefore captures the effects of simultaneous association between lithium ions and multiple ether groups. Third, our molecular simulations incorporate compositional fluctuations, which are absent from the mean-field treatment of refs 4 and 8. We emphasize that neglecting dielectric heterogeneity is an important approximation and note that including it in the free energy calculations described in the next section requires highly computationally demanding sampling calculations. We thus choose to focus on the weak segregation regime, where the strength of the dielectric heterogeneity is small. In forthcoming work, focused on ordered copolymer phases, the dielectric heterogeneity, Born solvation, and its effects will be described in detail. Further, to fully account for the effects of strong binding between lithium ions and ether groups, it would be necessary to pursue separate simulations that model lithium ions as bonded to PEO-like blocks (yet allowed to anneal, as in ref 4), whose results could then be directly compared with those reported in refs 4 and 8.

■

FREE ENERGY FROM EXE-DOS Our main concern in this work is to investigate the order− disorder transition in lamellae forming diblocks with sufficient precision to distinguish “dilution” contributions from ionic association effects on the free energy. These two effects can compete with one another, leading to opposing effects on the phase behavior of material. For concreteness, we focus on the value of χABN at the transition point, (χABN)ODT, and investigate systematically its variation with the addition of salt. In previous simulations, the ODT value has been estimated by means of a temperature ramp sweep, by monitoring the values of the heat capacity, or by visual inspection of morphologies.23 Since the disordered-to-lamellar transition is a weakly first-order transition,24−28 these approaches can suffer from the effects of hysteresis near the transition.23 To arrive at a more accurate estimate of the ODT, we adopt an approach based on the free energy, which is inspired by Leibler’s expansion to the free energy of mesophases in block copolymers.29 Specifically, we construct an order parameter based on the Fourier modes of the composition fields. In a disordered phase, these fields include contributions from Fourier modes of different wavelengths, which reveal themselves in real space as bead number fluctuations with different periodicity and amplitude. The free energy cost of these fluctuation modes is always positive, so the temporal composition modulation decays quickly. As χAB is increased, D

DOI: 10.1021/acs.macromol.5b02643 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules

like block to the cations χACN, and the Bjerrum length lB. We focus on symmetric diblock polymers of moderate length with N̅ = 256 and N = 8. The compressibility parameter is set to κN = 50. The box dimensions are 4 × 4 × 4 Re3, and the number of beads in the system is 8192. Depending on salt concentration, the typical number of polymers is 1000. Sampling the free energy is computationally demanding, and in this initial study we do not attempt to use higher resolution models (larger N values). Several representative density profiles for lamellar phases at different χABN values are shown in Figure 3. Panel a shows the distribution of all four species at χABN = 22. Both cations and anions are concentrated inside the PEO-like domains. The cations are attracted into these domains by associative interactions with PEO-like beads, and anions are dragged into these domains by Coulombic interactions. The distributions of cations and anions are nearly identical, resulting in a nearly neutral overall charge density profile. Panels b−d represent the variation of bead distributions with χABN. As χABN increases, the system becomes more strongly segregated, and ions become increasingly localized. In contrast to the reported trends for neutral and weakly selective solvents that are localized in the interfacial regions,36 our results suggest that charged, and more selective, ions tend to be concentrated in the middle of the domains. We have, however, verified that the neutral and nonselective solvents are slightly more concentrated in the interfacial regions. The structure factor defined by compositional fluctuation correlations of PEO and PS-type monomers are presented in Figure 4, for different χABN values at ϕsalt = 0, 0.031 25, and 0.0625. The structure factor exhibits a scattering peak at the characteristic wavenumber q*, which is inversely proportional to the polymer radius of gyration, similar to the behavior of neutral diblock copolymers.29 As temperature is decreased or χABN is increased, the peak intensity increases and the peak position shifts to smaller q values, which correspond to greater range of correlation for collective compositional fluctuations, consistent with results in neutral diblock copolymers.37 Adding salts into diblock copolymers shifts the peak positions to smaller q* values. This trend is consistent with experimental findings.6 Our main results on the variation of (χABN)ODT with ϕsalt and lB, obtained using the EXE-DOS free energy calculations, are presented in Figure 5. The curve labeled “dilution effect” shows the shift in (χABN)ODT for diblock copolymers diluted by a neutral solvent (χACN = 0). The solvent dilutes the repulsive interactions between the A and B blocks, thereby necessitating higher values of χABN to induce microphase separation. Our free energy methodology enables us to make a quantitative prediction of this effect in simulations. Conventional mean field theories based on the so-called dilution approximation38 predict that neutral solvent dilution shifts (χABN)ODT by a factor proportional to (1 − ϕsolvent)−1. Experiments, however, yield an exponent whose value is between −1.2 and −1.5,39 and the discrepancy has been attributed to the neglect of fluctuations. Our results, which include the effects of fluctuations, indicate that the effective dilution exponent is indeed −1.22, in agreement with experiment and independent simulations.17 For a selective solvent, the ODT also depends on the quality of the solvent. In our model, the affinity of the PEO-like block for the cations is represented by a negative pairwise interaction χAC. We do not treat the lithium ion as bound to the PEO block, thereby preserving the translational entropy of the

the value of the order parameter ψ and a histogram of visits to specific values of ψ are maintained. We implement the standard Metropolis Monte Carlo sampling with the system Hamiltonian and with the bias potential, and the histogram of ψ is accumulated during the course of simulations. Every time a particular value of order parameter is visited, the bias potential is updated to a fixed quantity Δ, which is typically initiated to Δ = 1 kBT. Once the histogram becomes sufficiently flat, we reduce the increment Δ by a factor of 2. In the end, when the histogram has become sufficiently flat and smooth, the inverse of the bias potential provides a direct measure to the free energy. For our systems, the free energy curves converge after approximately 18 iterations. The representative free energy curves for different values of χABN of a diblock copolymer are provided in Figure 2. For χABN values smaller than 17.8, the free energy curve exhibits only one minimum close to ψ = 0.15, corresponding to the disordered phase. For χABN values around 18, a well-defined second minimum appears near ψ = 0.25, corresponding to the lamellar phase. A distinct free energy barrier appears between the two minima, showing that the transition is indeed first order. From the figure, we can see that the (χABN)ODT is close to 18.6, at which point the two minima are nearly identical (the extrapolated actual value is 18.3). For higher values of χABN, the lamellar minimum exhibits a lower free energy value. The box size we used in all reported calculations is 4 × 4 × 4 Re3. We made this choice based on careful calibration studies on the effects of box dimension. The values of (χABN)ODT obtained from the same system (ϕsalt = 0.0625) using different box dimensions are tabulated in Table 1. The values of order Table 1. Order−Disorder Transition Temperature Obtained Using the Nonbiased Order Parametera N̅ 1/2

N

M

L/Re

16 16 16 16 16

8 8 8 8 8

303 303 303 1024 1024

2.7 2.7 2.7 4.0 4.0

max index (qRe) 3 4 5 4 5

(7.0) (9.4) (11.8) (6.3) (7.9)

(χABN)ODT

ψ range

18.4−18.6 18.5−18.6 18.4−18.5 18.2−18.4 18.2−18.4

0.1−0.2 0.1−0.2 0.1−0.2 0.1−0.2 0.1−0.2

a

The max index refers to the maximum integer wave index, which multiplied by 2π/L equals the maximum wavenumber. For instance, if the max index is 3, the max wave vector used for the definition of order parameter is (2π/L)[111]. The value of q* is between 3 and 4 Re−1.

parameters are also reported. It is apparent that a choice of 4 Re is sufficient for eliminating the finite size effect. The main limitation of our choice of order parameter for EXE-DOS free energy sampling is that it was conceived for systems in a uniphase thermodynamic state and therefore cannot be used to quantify the binodal curve. The latter has been found to exhibit interesting features in salt-doped diblock copolymers35 and could be investigated in future work by relying on more demanding, iterative calculations of the chemical potentials of all species in each phase.

■

VARIATION OF ODT We now examine systematically the variation of ODT with the addition of salt and discuss the effects of solvent dilution, ion selectivity, and the combined effects of salt concentration and long-range Coulomb interactions. The main controlling parameters are salt concentration ϕsalt, selectivity of the PEOE

DOI: 10.1021/acs.macromol.5b02643 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules

Figure 3. Bead density profiles for parameters N̅ = 256, N = 8, lB = 8.0, and ϕsalt = 0.0625. The distributions of cations and anions nearly overlap.

Figure 4. Structure factors of disordered morphologies at χABN = 2, 4, 6, 8, 10, 12, 14, and 16 for parameters N̅ = 256, N = 8, and lB = 8.0 for concentrations ϕsalt = 0, 1/32, and 1/16. The peak positions q* are determined by fitting the structure function to cubic splines and numerically finding the maximum position.

cations,4 and our formalism naturally allows the ions to be associated with multiple PEO strands because the contribution to density fields at each grid site may come from different chains. Our results (Figure 2b of Supporting Information) indicate that (χABN)ODT decreases with χACN. This behavior can be understood as follows: phase separation results in the formation of PEO-rich domains, and placing ions inside those domains reduces the energetic contribution to the free energy. In what follows, when examining the effects of Bjerrum length and salt concentration, we fix the value of χACN at −40, so that cations are preferentially bound to the PEO-rich domains, but remain mobile. In the Supporting Information, we also present results for χACN = −60 and for lB = 8. The variation of (χABN)ODT with salt concentration is analogous to that observed for χACN = −40, although the magnitude of the variations is more pronounced. We focus on modest values of χACN = −40 because the main purpose of this paper is to demonstrate the trends and because using a larger value of χAC perturbs the system’s morphology considerably (see Figure 1 of

Figure 5. Left: shift in (χABN)ODT as a function of salt concentration. The top curve corresponds to a neutral nonselective solvent. The data sets labeled by different lB values represent different Coulombic interaction strengths, all obtained for a solvent selectivity χACN = −40 (N̅ = 256). Right: representative lamellar morphologies from simulations at salt concentrations ϕsalt = 0.031 25 and ϕsalt = 0.1875 (lB = 8.0). Cyan domain: PS; pink domain: PEO; green and purple spheres: cations and anions. F

DOI: 10.1021/acs.macromol.5b02643 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules

exponent −1.22. For diblock copolymers diluted with salts, in the low dilution regime, (χABN)ODT decreases, consistent with experimental results. In the high dilution regime, (χABN)ODT increases due to the solvent dilution effects. When the Coulomb interaction is weak (e.g., when lB = 0.5), these two effects compete with each other. The calculations presented here provide a glimpse of the order−disorder transition in salt-doped polymers with a model explicitly conceived and designed for ionic systems. It treats the Coulomb interactions explicitly, thereby providing the ability to fully incorporate fluctuation effects, in contrast to previous theoretical approaches. An important limitation of our model (and most models to date) that must still be addressed is related to the absence of dielectric inhomogeneity. Wang et al. have addressed this issue in the context of a mean field approach;8 at this point, the two approaches yield qualitatively consistent results. However, these two results have different origins. In ref 8 the variation of the ODT was mainly attributed to the Born solvation energy, which favors the solvation of ions in the medium having a higher dielectric permittivity.4,7 In our work, the variation of the ODT mainly comes from the associative (χAC < 0) interaction between lithium ions and PEO-like segments that may belong to different chains. Within the framework of the current model, we cannot investigate the effects of Born solvation, whereas the earlier theory8 neglected the simultaneous complexation of lithium ions with multiple PEO segments. The consequences of these two effects are similar, but it is still unclear which factor is more important in realistic systems. Recently, the theory described in ref 8 has been extended to incorporate the effects of cross-linking lithium ions with multiple PEO segments in polymer blends,42 where it was found that the value of the effective χ increases, implying a decrease in (χN)ODT that is consistent with our observations. Further simulation work that incorporates explicitly the dielectric heterogeneity by solving the full Poisson’s equation with realistic dielectric profiles would help address this question. The results of additional calculations along those lines will be presented in a forthcoming publication that will focus on the morphological behavior of ordered phases. It has also been suggested in the literature that adding salts to block copolymers will blur the order−disorder transition boundary, giving rise to a two-phase coexistence region.8,35 This coexistence is dictated by the Gibbs phase rule. What is particularly interesting is that the coexistence window is found to be of the order of 11 °C, and the domain spacing of the ordered phase apparently increases with temperature in the coexistence window.35 We have not attempted to delineate that coexistence curve in our simulations, which have relied on a uniphase density-of-states approach. Different types of more demanding simulations that sample directly the chemical potentials at varying compositions would be needed for that. Such calculations would be pursued in future work. In summary, we have provided a study of the ordering transition in salt-doped symmetric diblock copolymers by relying on a highly coarse-grained model and a density-of-states based free energy sampling technique. Our results, alongside with those from the literature,8 suggest that the ordering transition in ionic polymeric systems is affected by multiple factors: association of Li-like cations with PEO-like blocks (χAC N), solvent dilution (ϕ salt), strength of Coulomb interactions (lB), and selective solvation due to Born solvation.8 Incorporating more microscopic details, such as ion clustering,43 nonlinear dielectric screening,44,45 and ion polar-

the Supporting Information for the morphology corresponding to χACN = −320). It is, however, worth pointing out that saturation of associative interactions between lithium ions and ether oxygens can not be described by a negative χ parameter. The shift in the ODT as a function of salt concentration for different values of the Bjerrum length is also shown in Figure 5. The dashed lines are a guide to the eye, connecting data points corresponding to the same set of parameters. The data set labeled “lB = 0” shows results for selective solvents (χACN = −40). The other three sets of data correspond to the results of different Coulomb interaction strengths, for lB values ranging from 0.5 to 16. The Bjerrum length is measured in terms of the polymer end-to-end distance Re, which may be related to the domain spacing d* close to the ODT by Re ≃ [(1.95√6)/2π] d* = 0.76d*.29 In experiments, the values of lB are of the order of d* or Re. For instance, the domain spacings for the two block copolymers used in ref 35 are close to 7.61 and 14.9 nm, respectively, whereas the Bjerrum length calculated at room temperature is 21.5 nm for pure PS with ϵr = 2.640 and 7.5 nm for pure PEO with ϵr = 7.4.41 We first examine the results in the low-salt regime (ϕsalt < 0.1), which corresponds to the upper limit of most experimentally accessible salt concentrations. For the two data sets with lB = 8 and 16, corresponding to strong Coulomb interactions, (χABN)ODT decreases with salt addition, in agreement with previous experimental and theoretical studies.5,6,8 The data set with lB = 0.5 sits between the results for lB = 8 and those for the neutral solvent, suggesting that when the Coulomb interaction is weak, the attractive interaction between the cations and anions is not strong enough to bring all anions into the PEO-like domain; in that regime, the solvent dilution effect competes with the ion-association effect, thereby causing (χABN)ODT to increase slightly. In the concentrated regime (ϕsalt > 0.1), it is clear that for all the systems considered here (χABN)ODT increases with salt concentration. This is directly related to the solvent dilution effect, which is dominant in the neutral solvent case. In this highly concentrated regime, the highly swollen PEO-like domains cannot retain all of the ions, and some of them leak to the PS-like domains, effectively diluting the interaction between PEO and PS blocks and causing (χABN)ODT to increase. The nature of ion distributions in the two domains can be appreciated in Figure 5 for ϕsalt = 0.031 25 and ϕsalt = 0.1875, respectively, from simulated configurations (lB = 8). In both cases, the cations condense into the PEO domains due to the strong ion association effect, and the anions are dragged into the PEO domains by the electrostatic attraction. In the high-salt case, for entropic reasons, some of the ions can also reside inside the PS domain. Note, however, that the limited solubility of lithium salts in PEO-like domains makes the highsalt regime less relevant for applications and is only included here for completeness.

■

SUMMARY In summary, we have proposed a coarse-grained model that incorporates Coulomb interactions and fluctuations explicitly. By relying on a density-of-states free energy sampling scheme built around a structure-factor based order parameter, the model has been applied to study the effects of salt addition and fluctuations on the order−disorder transition of diblock copolymers. We have found that when diblock copolymers are diluted with neutral solvents, (χABN)ODT increases according to a power law of the polymer concentration with G

DOI: 10.1021/acs.macromol.5b02643 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules izability,46 would require a simulation model with hard-core potentials and with reliable polarizable force fields. Future atomistic simulations along those lines will be needed to examine some of the assumptions implicit in our coarse-grained approach. We conclude by noting that a recent mean-field approach8 and our coarse-grained simulations, which incorporate fluctuations, are complementary. They both lack the molecular packing details that are needed for a quantitative description of specific experimental systems, and they should both be revisited as additional information from experiments and atomistic simulations becomes available.

■

(10) Sing, C. E.; Zwanikken, J. W.; Olvera de la Cruz, M. Ion correlation-induced phase separation in polyelectrolyte blends. ACS Macro Lett. 2013, 2, 1042. (11) Sing, C. E.; Zwanikken, J. W.; Olvera de la Cruz, M. Electrostatic control of block copolymer morphology. Nat. Mater. 2014, 13, 694. (12) Ganesan, V.; Pyamitsyn, V.; Bertoni, C.; Shah, M. Mechanisms underlying ion transport in lamellar block copolymer membranes. ACS Macro Lett. 2012, 1, 513−518. (13) Knychala, P.; Dziȩcielski, M.; Banaszak, M.; Balsara, N. P. Phase behavior of ionic block copolymers studied by a minimal lattice model with short-range interactions. Macromolecules 2013, 46, 5724−5730. (14) Knychala, P.; Banaszak, M.; Balsara, N. P. Effect of composition on the phase behavior of ion-containing block copolymers studied by a minimal lattice model. Macromolecules 2014, 47, 2529−2535. (15) Pike, D. Q.; Detcheverry, F. A.; Müller, M.; de Pablo, J. J. Theoretically informed coarse grain simulations of polymeric systems. J. Chem. Phys. 2009, 131, 084903. (16) Müller, M. Concurrent coupling between a particle simulation and a continuum description. Eur. Phys. J.: Spec. Top. 2009, 177, 149− 164. (17) Hur, S.; Khaira, G. S.; Ramírez-Hernández, A.; Müller, M.; Nealey, P. F.; de Pablo, J. J. Simulation of defect reduction in block copolymer thin films by solvent annealing. ACS Macro Lett. 2014, 4, 11−15. (18) Jha, P. K.; Zwanikken, J. W.; Detcheverry, F. A.; de Pablo, J. J.; Olvera de la Cruz, M. Study of volume phase transitions in polymeric nanogels by theoretically informed coarse-grained simulations. Soft Matter 2011, 7, 5965. (19) Marshall, S. L. A periodic Green function for calculation of coulombic lattcie potentials. J. Phys.: Condens. Matter 2000, 12, 4575− 4601. (20) Essén, H.; Nordmark, A. Some results on the electrostatic energy of ionic crystals. Can. J. Chem. 1996, 74, 885−891. (21) Villet, M. C.; Fredrickson, G. H. Efficient field-theoretic simulation of polymer solutions. J. Chem. Phys. 2014, 141, 224115. (22) Riggleman, R. A.; Kumar, R.; Fredrickson, G. H. Investigation of the interfacial tension of complex coacervates using field-theoretic simulations. J. Chem. Phys. 2012, 136, 024903. (23) Beardsley, T. M.; Matsen, M. W. Monte Carlo phase diagram for diblock copolymer melts. Eur. Phys. J. E: Soft Matter Biol. Phys. 2010, 32, 255−264. (24) Brazovskiĭ, S. A. Phase transition of an isotropic system to a nonuniform state. Sov. Phys. JETP 1985, 41, 85−89. (25) Fredrickson, G. H.; Helfand, E. Fluctuation effects in the theory of microphase separation in block copolymers. J. Chem. Phys. 1987, 87, 697−705. (26) Barrat, J.-L.; Fredrickson, G. H. Collective and single chain correlations near the block copolymer order-disorder transition. J. Chem. Phys. 1991, 95, 1281−1289. (27) de la Cruz, M. O. Transitions to periodic structures in block copolymer melts. Phys. Rev. Lett. 1991, 67, 85−88. (28) Mayes, A. M.; de la Cruz, M. O. Concentration fluctuation effects on disorder-order transition in block copolymer melts. J. Chem. Phys. 1991, 95, 4670−4677. (29) Leibler, L. Theory of microphase separation in block copolymers. Macromolecules 1980, 13, 1602−1617. (30) Glaser, J.; Medapuram, P.; Beardsley, T. M.; Matsen, M. W.; Morse, D. C. Universality of block copolymer melts. Phys. Rev. Lett. 2014, 113, 068302. (31) Singh, S.; Chopra, M.; de Pablo, J. J. Density of States-based molecular simulations. Annu. Rev. Chem. Biomol. Eng. 2012, 3, 369− 394. (32) Wang, F.; Landau, D. P. Efficient, multiple-range random walk algorithm to calculate the density of states. Phys. Rev. Lett. 2001, 86, 2050−2053. (33) Yan, Q.; de Pablo, J. J. Fast calculation of the density of states of a fluid by Monte Carlo simulations. Phys. Rev. Lett. 2003, 90, 035701.

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.macromol.5b02643.

■

Choice of model parameters (PDF)

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected] (J.J.d.P.). Present Address

J.Q.: Department of Chemical Engineering, Stanford University, Stanford, CA 93405. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS We acknowledge helpful discussions with Jens Glaser and David Morse for details regarding the choice of order parameter. This work was supported by the Department of Energy, Basic Energy Sciences, Materials Sciences and Engineering Division.

■

REFERENCES

(1) Meyer, W. H. Polymer electrolytes for lithium-ion batteries. Adv. Mater. 1998, 10, 439−448. (2) Hallinan, D. T.; Balsara, N. P. Polymer electrolytes. Annu. Rev. Mater. Res. 2013, 43, 503−525. (3) Young, W.-S.; Kuan, W.-F.; Epps, T. H. Block copolymer electrolytes for rechargeable lithium batteries. J. Polym. Sci., Part B: Polym. Phys. 2014, 52, 1−16. (4) Nakamura, I.; Wang, Z.-G. Salt-doped block copolymers: ion distribution, domain spacing and effective χ parameter. Soft Matter 2012, 8, 9356−9367. (5) Soo, P. P.; Huang, B.; Jang, Y.-I.; Chiang, Y.-M.; Sadoway, D. R.; Mayes, A. M. Rubbery block copolymer electrolytes for solid-state rechargeable Lithium batteries. J. Electrochem. Soc. 1999, 146, 32−37. (6) Gunkel, I.; Thurn-Albrecht, T. Thermodynamic and structural changes in ion-containing symmetricdiblock copolymers: a small-angle X-Ray scattering study. Macromolecules 2012, 45, 283−291. (7) Nakamura, I.; Balsara, N. P.; Wang, Z.-G. Thermodynamics of ion-containing polymer blends and block copolymers. Phys. Rev. Lett. 2011, 107, 198301. (8) Nakamura, I.; Balsara, N. P.; Wang, Z.-G. First-order disordereto-lamellar phase transition in lithium salt doped block copolymers. ACS Macro Lett. 2013, 2, 478−481. (9) Sing, C. E.; Zwanikken, J. W.; Olvera de la Cruz, M. Interfacial behavior in polyelectrolyte blends: hybrid liquid-state integral equation and self-consistent field theory study. Phys. Rev. Lett. 2013, 111, 168303. H

DOI: 10.1021/acs.macromol.5b02643 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules (34) Kim, E. B.; de Pablo, J. J. Potential of mean force between a spherical particle suspended in a nematic liquid crystal and a substrate: sphere size effects. Phys. Rev. E 2004, 69, 061703. (35) Thelen, J. L.; Teran, A. A.; Wang, X.; Garetz, B. A.; Nakamura, I.; Wang, Z.-G.; Balsara, N. P. Phase behavior of a block copolymer/ salt mixture through the order-to-disorder transition. Macromolecules 2014, 47, 2666−2673. (36) Huang, C.; Olvera de la Cruz, M. Adsorption of a minority component in polymer blend interfaces. Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top. 1996, 53, 812−819. (37) Qin, J.; Morse, D. C. Fluctuations in symmetric diblock copolymers: testing theories old and new. Phys. Rev. Lett. 2012, 108, 238301. (38) Fredrickson, G. H.; Leibler, L. Theory of block copolymer solutions: nonselective good solvents. Macromolecules 1989, 22, 1238− 1250. (39) Hamley, I. W. Block Copolymers in Solution: Fundamentals and Applications; John Wiley & Sons: 2005. (40) Yano, O.; Wada, Y. Dynamic mechanical and dielectric relaxations of polystyrene below the glass temperature. J. Polym. Sci. Part A-2 1971, 9, 669−686. (41) Porter, C. H.; Boyd, R. H. A dielectric study of the effects of melting on molecular relaxation in poly(theylene oxide) and polyoxymethylene. Macromolecules 1971, 4, 589−594. (42) Ren, C.-L.; Nakamura, I.; Wang, Z.-G. Effects of ion-induced cross-linking on the phase behavior in salt-doped polymer blends. Macromolecules 2016, 49, 425−431. (43) Hall, L. M.; Seitz, M. E.; Winey, K. I.; Opper, K. L.; Wagener, K. B.; Stevens, M. J.; Frischknecht, A. L. Ionic aggregate structure in ionomer melts: effects of molecular architecture on aggregates and the ionomer peak. J. Am. Chem. Soc. 2011, 134, 574−587. (44) Nakamura, I.; Shi, A.-C.; Wang, Z.-G. Ion solvation in liquid mixtures: effects of solvent reorganization. Phys. Rev. Lett. 2012, 109, 257802. (45) Nakamura, I. Ion solvation in polymer blends and block copolymer melts: effects of chain length and connectivity on the reorganization of dipoles. J. Phys. Chem. B 2014, 118, 5787−5796. (46) Freed, K. F. Perturbative many-body expansion for electrostatic energy and field for systems of polarizable charged spherical ions in a dielectric medium. J. Chem. Phys. 2014, 141, 034115.

I

DOI: 10.1021/acs.macromol.5b02643 Macromolecules XXXX, XXX, XXX−XXX