The evolution of acidic and ionic aggregates in ionomers during microsecond simulations

Single-ion conducting polymers, such as ionomers which contain ionic groups covalently bonded to the polymer backbone, continue to be of interest as electrolytes since in these polymers the conductivity is solely due to the unbound ion.^1,2 In typical ionomers, the ionic groups self-assemble into aggregates that are nano-phase separated from the polymer backbone.³ Ion transport is strongly affected by the morphology of these ionic aggregates. This has perhaps best been demonstrated in hydrated ionomers such as Nafion and related proton-conducting membranes, where percolation of the hydrated ionic domains is needed to obtain good conductivity.⁴ 

In melt ionomers (with no solvent), experimental relationships between ion conductivity and ionic aggregate morphology are mostly unexplored. This is largely because it is difficult to determine the morphology from experiments alone. Microscopy images of ionomers, such as those from scanning transmission electron microscopy (STEM), are complicated by projection issues, although in a few cases the morphology can be determined.^5–8 X-ray scattering from the ionic aggregates generally shows a single broad peak at low wavevector, known as the ionomer peak.^9,10 In general, there is not sufficient information in this peak to determine details of the morphology.¹¹ One solution to this problem is to use validated molecular dynamics (MD) simulations to determine the morphology.

We have been following this strategy of using MD simulations and X-ray scattering together to understand the morphology in precise poly(ethylene-co-acrylic acid) ionomers.^8,11–16 These precise materials are made by acyclic diene metathesis (ADMET) polymerization.¹⁷ They are strictly linear polyethylenes, with carboxylic acid groups separated by a precise number of carbon atoms along linear polymer backbones. In a set of poly(ethylene-co-acrylic acid) copolymers and ionomers with a range of neutralization levels and cations, we found good agreement between the structure factors calculated from MD simulations and those measured by X-ray scattering.^11,15 The MD simulations revealed a range of ionic aggregate morphologies, from isolated, compact aggregates, to stringy aggregates, to systems with percolated aggregates that spanned the simulation box in all directions. These distinct morphologies give rise to very similar ionomer peaks, making it impossible to distinguish different morphologies by scattering alone.

We have also performed MD simulations on a coarse-grained model of these polymers. These simulations show that ion transport is faster in systems that form percolated ionic aggregates than in systems with isolated ionic aggregates.^18,19 A similar improved transport has been seen in other simulations of ionomers with percolated ionic aggregates.^20,21 Additionally, recent simulations and experiments on precise polyethylene-acid copolymers at lower temperatures showed a unique layered morphology in semi-crystalline systems,²² which have good conductivity when sulfonated and hydrated.²³ It thus seems that spatially continuous, extended aggregates are advantageous for fast ion conduction. Consequentially, understanding the detailed morphology of ionic aggregates in ionomers is essential for designing improved single-ion conducting polymers.

Our previous atomistic MD simulations^11,14,15 of these precise polyethylenes and ionomers were quite short, on the order of 30 ns duration, and we did not attempt to investigate the dynamics of the ionic aggregates. These dynamics can be very slow, as demonstrated by Lu et al. in a similar ionomer system.²⁴ Here, we extend our previous atomistic simulations of precise Li ionomers to more than 1 μs at two temperatures. We find that at the higher temperature, the ionic aggregates rearrange and slowly grow larger with time. Additionally, the composition of the ionic aggregates is somewhat different than what we had reported previously. Here we find that in partially neutralized systems, the acid groups decorate the outside of a core of lithium ions and oxygen atoms in the ionic aggregates and also form an independent population of acid-only aggregates. These longer simulations provide a more complete view of the ionic aggregate morphology in precise poly(ethylene-co-acrylic acid) ionomers.

In this paper, we focus on a set of poly(ethylene-co-acrylic acid) precise polymers neutralized with lithium ions. The system nomenclature is pxAA-y%Li, where p indicates a precise polymer, x is the carbon spacing length (x = 9, 15, or 21) between pendant acrylic acid (AA) groups, and y is the percent the system is neutralized with Li⁺ ions. (In figure legends, we will drop the “AA” and refer to the systems as px-y%Li for brevity.) Here we simulate systems neutralized either partially or fully with Li⁺; the six systems simulated are listed in Table I. The backbone chain lengths were kept roughly constant, while preserving the precise spacing of acid groups along the chains. The partial neutralization levels were chosen to match those of a set of ionomers for which we have X-ray scattering data and also quasi-elastic neutron scattering (QENS) data, which will be the subject of a future paper. Simulations were performed at both 423 K and 600 K. All simulation boxes were cubic, with box lengths L, as reported in Table I. Figure 1 shows the structure of the acid copolymer and one of the partially neutralized ionomers.

Our previous atomistic MD simulations of precise acid-functionalized polyethylenes and precise ionomers, neutralized to varying degrees with metal cations (Li, Na, Cs, and Zn), were performed at two temperatures (393 K and 423 K) using the all-atom optimized potentials for liquid simulations (OPLS-AA) fully atomistic force field.^25,26 Here, we follow our more recent simulations of the acid copolymers²⁷ and perform similar simulations with the recent LOPLS-AA force field, which was developed by Siu et al. in 2012 to improve properties of long alkanes.²⁸ We showed previously that this force field gives dynamics in excellent agreement with QENS data for polyethylene and for the precise acid copolymers.²⁷ We note that 600 K is a high temperature, above the physical degradation temperature of the ionomers; we employed this temperature to allow sufficient ion motion to analyze the time dependence of the ionic aggregate distributions. The LOPLS-AA force field has been shown to be accurate up to at least 477 K (e.g., for polyethylene dynamics as compared with QENS data²⁷) and its accuracy is not likely to be strongly temperature dependent for these precise ionomers.

All simulations used the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) software package,^29,30 with a real-space nonbonded cutoff of 12 Å and the particle-particle-particle-mesh (PPPM) solver for electrostatics. A Nosé-Hoover thermostat with a 100.0 fs damping parameter was used to maintain the temperature at either 423 K or 600 K. The integration time step was 1.0 fs in all cases. Initial states for each system were created using configurational-bias Monte Carlo simulations via the Enhanced Monte Carlo (EMC) package,^31,32 which builds the system near the desired density. Simulations were run in the NPT ensemble at 423 K and 1 atm for 2 ns, followed by at least 10 ns in the NVT ensemble, after which trajectories were analyzed. After approximately 100 ns at 423 K, the temperature was ramped up to 600 K for each system. Systems were further equilibrated at 600 K in the NPT ensemble for 5 ns to equilibrate the density. Further equilibration then took place in the NVT ensemble for at least 15 ns, after which trajectories were analyzed.

Structure factors were calculated from the Fourier transform of the radial distribution functions, g_ij(r), for each atomic type, using

The revised Lorch window function given by

was used to remedy the cutoff ripple artifact due to the Fourier transform.^33–35 This was important as not all g_ij(r) → 1 at large r in these simulations. In Eq. (1), i and j are the indices used to count atomic species (here, carbon, hydrogen, sulfur, oxygen, and lithium), ρ is the total number density of the system, c_i and c_j are the mole fractions of species i and j, and f_i is the atomic scattering function for species i. The atomic scattering functions f_i typically depend on q. This dependence has been empirically fitted for all atom types for X-ray scattering, using the following form:

We use the data of Waasmaier and Kirfel³⁶ for the coefficients (c_i, a_ik, b_ik) in Eq. (1). S(q) is calculated for wavevectors q ≥ 4π/L, where L is the simulation box length. This is because the information in g(r) extends only to r = L/2, limiting the longest mode to L/4 and hence the smallest accessible wavevector in S(q) to 4π/L.

As a check on the structure factors, they were also calculated directly using

with the total structure factor given by

Here, N_α is the number of atoms of type α and N is the total number of atoms in the simulation. The results are equivalent to those of Eq. (1) but noisier, especially at low q; see Fig. S1 in the supplementary material for two examples. Note that the lowest accessible wavevector in the direct calculation is 2π/L, but S(q) is fairly noisy at low q due to poor statistics at low q values.

Cluster analysis was performed as described previously.¹⁵ In brief, all Li ions, oxygen atoms, and acid hydrogen atoms that are within a given cutoff distance of each other are considered to belong to the same aggregate; furthermore, all oxygens and acid hydrogens that belong to the same carboxylate/carboxylic acid group are also considered to be in the same aggregate. The cutoffs are based on the trough between the first and second peaks in the relevant g_ij(r). Visualizations were performed with VMD.^37–39 

As in our previous work, our MD simulations produce a variety of ionic aggregate shapes depending on the spacer length and the degree of neutralization, at both temperatures. In particular, the p9AA-100%Li system, which has the highest ion content, forms a percolated ionic aggregate at both temperatures, whereas the p21AA-38%Li system, with the lowest ion content, forms compact, isolated clusters. Examples of these two morphologies are shown in Fig. 2 for T = 600 K; similar aggregate morphologies were observed at 423 K.

A previously unexplored question is, how equilibrated are these ionic aggregates? The dynamics in these systems are slow, which motivated us to do simulations at the higher temperature of 600 K, in addition to the more experimentally relevant temperature of 423 K. To quantify the time evolution of the aggregates, we calculated the average sizes of the aggregates as a function of time, sampled every 10 ns. The aggregate size changes over time with substantial changes occurring even 100 s of ns into the simulation. All three of the 100%Li systems at 600 K eventually evolve percolated ionic aggregates that span the simulation box in all three directions. The p9AA-100%Li system is percolated at all times, while the p15AA-100%Li and the p21AA-100%Li systems only become percolated after about 120 ns and 480 ns, respectively. Figure 3 shows the size of the largest cluster in each of these systems as a function of time. This largest cluster size has been normalized by the total number of Li ions in the system, so we plot this largest cluster size as a fraction of the largest possible size (which would contain all the Li ions). By the end of the simulation trajectories, these three systems appear to have reached a steady state cluster distribution, in which >90% of the Li ions are in a single percolated cluster, with 1 or 2 smaller isolated clusters in p15AA-100%Li and p21AA-100%Li, respectively. For the latter parts of the simulation, these systems have reached a steady state with regards to the aggregate morphology, namely, percolated, box-spanning aggregates.

The partially neutralized systems also show an evolution in the aggregate distribution with time. These systems are more complicated because there are two kinds of ionic aggregates. First, there are aggregates containing Li ions and COO⁻ groups; these aggregates typically also include some unneutralized acid groups (COOH). Second, there are aggregates consisting solely of unneutralized COOH groups that hydrogen bond to each other. We note that our previous study of the acid copolymers (with no lithium ions) showed that at 423 K, the acid aggregates rearrange on an approximately 2 ns time scale.²⁷ Aggregates involving lithium ions have much longer rearrangement time scales. To characterize the time dependence of the aggregate morphologies, we look at aggregates containing at least one Li ion. Figure 4(a) shows the fractional maximum cluster size at 600 K, again measured by the number of Li ions in the largest cluster divided by the total number of Li ions in the simulation box, as a function of time. Here the largest aggregate only contains a relatively small fraction of the total number of ions in the system. This largest aggregate grows in each ionomer with time, by the merger of smaller aggregates, as demonstrated by the decrease in the number of clusters with time, as shown in Fig. 4(b). The partially neutralized systems also appear to have reached a steady-state in the latter part of the simulation trajectories, particularly with regard to the largest Li ion-containing aggregates in the system, although the aggregates continue to slowly coarsen by occasional mergers of smaller aggregates into larger ones.

The specific ions in the aggregates change with time, due to a process of aggregates merging and breaking up, similar to what we have observed in coarse-grained simulations.¹⁸ Figure 5 shows only the Li ions and O atoms (in COO– groups) at t = 0 (measured from the start of the production run, after equilibration) and 1300 ns later, for p9AA-35%Li. Here, the aggregates are identified at t = 0 and each distinct aggregate is given its own color. We then retain the color of each ion or atom at t = 1300 ns; we find that there has been significant scrambling among the aggregates over the course of the simulation. Also, the trends shown in Fig. 4 are evident in Fig. 5, with an overall decrease in the number of Li-containing aggregates and an increase in aggregate size.

At the lower temperature of 423 K, the aggregates are less mobile and so the aggregate distributions show less change with time. The p9AA-100%Li system is also percolated at the lower temperature, with the largest percolated aggregate including nearly all of the ions in the system. However, the p15AA-100%Li system is initially only partially percolated and after about 200 ns develops a percolated aggregate, although other smaller aggregates remain. The p21AA–100%Li system only shows stringy aggregates that do not coarsen much during the simulation. These observations are quantified by the fractional maximum cluster size, as shown in Fig. 6. Recall that the initial states for the simulations were created by a short Monte Carlo simulation at 423 K, followed by equilibration for about 12 ns. The p9AA-100%Li system has so many ions that it immediately forms a percolated structure. The other two 100%Li systems, in contrast, are initially stuck in a nonequilibrium distribution of aggregates. The p15AA-100%Li system is able to coarsen somewhat, while the p21AA-100%Li system at 423 K cannot coarsen in the time frame of the simulations. We would expect the true equilibrium state of these three systems to consist of percolated aggregates; these are the steady states at the higher temperature of 600 K and surely must have lower free energy than the smaller aggregates that we observe at 423 K since the electrostatic interactions are very strong in the low dielectric medium of the polymer. We therefore conclude that the aggregate distributions for the 100% neutralized ionomers at 423 K are not in equilibrium, even after more than 1 μs of simulation time. For the partially neutralized systems, there is even less change in the aggregate distributions with time at 423 K (see Fig. S2 in the supplementary material).

The typical distance between aggregates, as manifested by the ionomer peak in the structure factor, is not particularly sensitive to the distribution of aggregate shapes. As examples of this, Fig. 7 shows the structure factor S(q) at 600 K for p9AA-35%Li and for p21AA-100%Li, calculated over two different time frames during the simulations. Scattering from these precise ionomers leads to two peaks. The higher wavevector peak at around 14 nm⁻¹ is known as the amorphous halo, and is the result of scattering between polymer backbones. The lower wavevector peak is the ionomer peak that comes from scattering between ionic aggregates.^11,15 In p9AA-35%Li, during the earlier time window, the system contains on average 31 aggregates, with mean size 8.5 Li ions, while during the later time window, the system contains on average 26 aggregates with mean size 10 Li ions. In p21AA-100%Li, during the earlier time window, the system contains 9 aggregates, ranging in size from 4 to 121 Li ions, while during the later time window, the system contains only 3 aggregates, the percolated aggregate with 290 Li ions and two smaller aggregates (with 25 and 5 Li ions). Nevertheless, the structure factor from these two distinct aggregate morphologies is nearly identical in both cases, as shown in Fig. 7. Here S(q) has been averaged over snapshots every 10 ns. As indicated by the location of the ionomer peak, the average distance between aggregates remains the same, regardless of the number and size of the aggregates (in p9AA-35%Li) or whether the aggregates are isolated or connected into one percolated aggregate (in p21AA-100%Li). Thus, the structure factor does not detect the degree of equilibration of the ionic aggregates, and the X-ray scattering data will be similarly limited.

This justifies comparing the structure factors calculated from the simulations at 423 K with our previously published X-ray scattering data.⁴⁰ We calculate S(q) by averaging over later times in the simulations, when the aggregate distributions are closer to the steady state; the specific averaging times for each system can be found in Table S1 in the supplementary material. Figure 8 shows S(q) calculated from the MD simulations for the partially neutralized systems at 423 K, compared with the measured X-ray scattering at 348 K.⁴⁰ The X-ray scattering curve intensities are in arbitrary units, so they have been scaled to the height of the amorphous halo in the MD simulations. Overall we see very good agreement between the simulations and the X-ray scattering, especially given that the temperatures are somewhat different. The amorphous halo is well-captured, as we have seen previously with these materials.¹¹ The ionomer peak is also remarkably well captured by the simulations, particularly for p9AA-35%Li and p21AA-38%Li. Overall the level of agreement between the MD simulation and X-ray scattering is similar to what we found for other pxAA copolymers and ionomers, neutralized at varying levels with a variety of cations.^11,15,27 Note that the location of the ionomer peak is well-captured by the simulations, but the accuracy of its height and width are limited by the small size of the simulation box and noisy statistics at low q; as shown in Fig S1, the ionomer peak intensity is not well-determined by the simulations.

The ionomer peak shifts to lower wavevector with increasing spacer length at similar neutralization levels, as seen in Fig. 9(a), as we have observed before.¹⁴ For these partially neutralized systems at 423 K, the ionomer peaks occur at q^* = 4.94, 3.94, and 3.32 nm⁻¹, for 9, 15, and 21 carbon spacers, leading to real space distances between isolated aggregates of d^* = 2π/q^* = 1.27, 1.59, and 1.89 nm, respectively. A larger spacer length leads to aggregates that are on average further apart. Increasing neutralization leads to a somewhat higher ionomer peak intensity for p9AA-y%Li systems and little change in ionomer peak position, as shown in Fig. 9(b) (d^* = 1.31 nm for p9AA-100%Li, compared to d^* = 1.27 nm for p9AA-35%Li, an increase of only 0.4 Å). For comparison, we include S(q) for the acid copolymer from our previous work.²⁷ In the acid copolymer, the ionomer peak is at somewhat higher wavevector, corresponding to a smaller average distance of d^* = 1.14 nm between acid aggregates. For the larger spacer lengths (15 and 21), increasing neutralization leads to only small changes in the ionomer peak (see Fig. S3). We see similar trends at the higher temperature of 600 K. A comparison between the scattering at 423 K and 600 K is shown for the p15AA systems in Fig. 10 [again calculating S(q) only at later times]. The amorphous halo shifts to somewhat lower wavevectors at higher temperature due to the slightly lower density at higher temperature. The ionomer peaks have somewhat larger intensity and a slight shift to lower q. The other two systems (p9 and p21) show similar trends (Fig. S4). At 600 K, where the aggregate morphologies appear to have reached steady-state, the location of the ionomer peak varies from q^* = 4.71 nm⁻¹ for p9AA-100%Li to q^* = 2.93 nm⁻¹ for p21AA-38%Li. Average distances between aggregates thus span from 1.33 to 2.14 nm.

Finally, we examine the nature of the aggregates in the partially neutralized systems at 600 K. As mentioned above, these systems contain two kinds of aggregates, mixed aggregates that contain lithium ions as well as acid groups, and aggregates consisting solely of acid groups. In our previous simulations carried out at 423 K over relatively short times (30 ns), we identified predominantly mixed aggregates with both lithium/oxygen repeating regions and hydrogen-bonded acid groups [see e.g., Fig. 8(b) in Ref. 15 and Fig. 1 in Ref. 11]. Many of the mixed aggregates consisted of clusters of lithium ions and oxygen atoms connected by hydrogen-bonded acid groups to another cluster of lithium and oxygens. A few aggregates of this type were also seen in the simulations reported here; an example is shown in Fig. 11(a). This aggregate is comprised of two lithium/oxygen aggregates connected by a single acid group (see arrow). As also mentioned above, the time scale for rearrangement of acid groups is much faster than for lithium-containing aggregates. In particular, our MD simulations of the acid copolymers at 423 K showed that the time scale for acid group rearrangement is on the order of 2 ns. This time scale is likely to be even shorter at 600 K. Since the simulations reported here were equilibrated and run for about 100 ns at 423 K and then equilibrated for an additional 20 ns at 600 K before analyzing properties at 600 K, the acid aggregate distributions are essentially at equilibrium at the start of the property analysis. At the higher temperature of 600 K, the lithium ions tightly pack with the charged oxygens in the central regions of the aggregates, while the acid groups are found on the perimeters of the aggregates (at the interface between the Li⁺/COO⁻ aggregates and the PE matrix). The result is long, stringy aggregates consisting of lithium/oxygen repeating motifs, decorated by acid groups at the ends and sometimes along the sides of the lithium/oxygen clusters, as shown in Fig. 11(b). In recent coarse-grained simulations of these ionomers that modeled the acid groups with “sticky” sites, the ionic aggregates had very similar morphologies, with a tightly-packed central core of ions and sticky groups decorating the perimeters of the clusters.⁴¹ 

Statistics for the composition of the Li⁺-containing aggregates are shown in Table II. Here we report the average number of Li ions and the average number of acid groups in the mixed Li-acid aggregates, over two different time windows for each system. The reported error is the standard deviation over that time window. In two cases, the standard deviation is zero because there was no change in the number of lithium ions in the aggregates over that time window. The last two columns in Table II report the average minimum and maximum sizes of the aggregates, measured by the number of Li ions. In all three partially neutralized ionomers, consistent with Fig. 4, the number of lithium ions in the aggregates increases with time as the aggregates merge. However, the number of acid groups associated with lithium aggregates remains constant within error. With increasing spacer lengths and hence decreasing ion concentration, both the mean and the maximum sizes of the aggregates decrease.

Overall, the morphologies in the partially neutralized systems tend to evolve to long, stringy aggregates consisting mostly of Li and O, with a few acid groups on the aggregate perimeters, coexisting with a population of purely acidic aggregates, as illustrated for p21AA-38%Li in Fig. 12. The purely acidic clusters tend to be considerably smaller than the aggregates containing lithium ions. The mean size of acidic aggregates is 2.5, 2.4, and 2.2 acid groups per aggregate for p9AA-35%Li, p15AA-38%Li, and p21AA-38%Li, respectively. These mean sizes do not change over the entire course of the simulations. The average maximum acid aggregate sizes are 5.6, 4.5, and 3.7 acid groups, for p9AA-35%Li, p15AA-38%Li, and p21AA-38%Li, respectively. On the other hand, the number of acid aggregates in the system does increase somewhat with time, as shown in Fig. S5. This is because the lithium-containing aggregates grow larger and their number decreases with time, while the number of acid groups associated with each lithium ion aggregate remains about the same (see Table II). The extra acid groups thus form additional acid-only aggregates.

Figure 13 shows the distribution of acid cluster sizes for the acid-only aggregates on both linear and semilog scales. The distributions are calculated over the entire simulation time; averaging over just later times gives very similar results. Here we have plotted the fraction of acid aggregates which consist of a given number of acid groups. The fraction of aggregates of a given size decreases monotonically with increasing size (numbers of acid groups in the aggregates). There are many single acid groups in all three systems, with a fairly substantial number of dimers, a smaller number of trimers, and decreasing numbers of larger acid aggregates. The number of aggregates with more than one acid group is larger for shorter spacer lengths, presumably because it is easier to form acid aggregates in those systems. The trends are similar to what is seen in the acid copolymers, with no lithium ions, except that in that case the fraction of single acid groups is lower (around 0.35 for p21AA), and there are more large acid aggregates, with up to 18 acid groups in an aggregate as opposed to about a maximum of 8 acid groups here, for the 9 carbon spacer.⁴² Thus, the presence of the lithium ions leads to ionic aggregates that incorporate some acid groups into the lithium ion aggregates, but the unneutralized acid groups then form their own distribution of acid aggregates that is similar to what is seen in the unneutralized acid copolymers.

We have performed microsecond-long simulations on a series of precise (polytethylene-co-acrylic acid) ionomers partially and fully neutralized with lithium, with three spacer lengths between acid groups on the ionomer and at two temperatures. At the lower temperature of 423 K, the ionic aggregate morphologies do not reach a steady-state distribution over the course of the simulations. At the higher temperature of 600 K, the aggregates are sufficiently mobile that they rearrange and appear to reach steady-state after 100 s of ns, depending on the system. For all three spacer length systems that are 100% neutralized with lithium, the ions form percolated aggregates that span the simulation box. In the partially neutralized systems, the morphologies include lithium ion aggregates that may also include some unneutralized acid groups, along with a coexisting population of acid aggregates. In the lithium ion aggregates, the acid groups decorate the perimeters of the central core of tightly-packed lithium ions and carboxylic groups. These distinct morphologies give rise to different ion and polymer dynamics, which will be the subject of a future manuscript.

See supplementary material for a comparison of the direct and Fourier methods for computing the structure factor, a table of averaging times, and additional cluster analysis and structure factor results.

We thank Jason Koski for running the MD simulations to long times and Benjamin Paren for helpful discussions. K.I.W. acknowledges support from NSF No. DMR 1506726. This work was performed, in part, at the Center for Integrated Nanotechnologies, an Office of Science User Facility operated for the U.S. Department of Energy (DOE) Office of Science. Sandia National Laboratories is a multi-mission laboratory managed and operated by National Technology and Engineering Solutions of Sandia, LLC., a wholly-owned subsidiary of Honeywell International, Inc., for the U.S. DOE’s National Nuclear Security Administration under Contract No. DE-NA-0003525.