Surface-grafted polyelectrolyte brushes (PEBs) are responsive to external stimuli, such as electric fields and ionic strength. 1 This responsiveness makes PEBs attractive for applications ranging from drug delivery to separations technologies, 2 but essential to PEB utilization is understanding how critical parameters like grafting density (σ) impact PEB structure and the dynamics of the PEB and counterions.
From a fundamental standpoint, counterion mobility in polyelectrolyte systems is vital because the interplay between the components (polymer, solvent, coions, and counterions) impacts transport in non-trivial ways. For instance, in Nafion (Dupont), the industry-leading polymer electrolyte membrane material, polymer phase separation between hydrophobic backbones and hydrophilic sidechains that provides ion pathways enables the high charge mobility required for applications such as fuel cell membranes and batteries. 3 Surface grafting of chain ends in PEBs adds complexity beyond phase separation by modifying polymer microstructure through topological constraints, especially at high σ. For these reasons, understanding counterion transport in PEBs is relevant for fundamental and application purposes and motivates this work.
Experimentally, counterion transport in PEBs has been investigated. Resnik 4,5 et al. used fluorescence correlation spectroscopy (FCS) to probe the diffusion of counterions using a cationic fluorescent probe (rhodamine-6G) in a strong polyanionic brush composed of poly(styrene sulfonate) (PSS). They demonstrated that counterion diffusivity within the PEB depended on electrostatic interactions between the cationic probe and anions on the PSS brush and steric hindrance in the grafted PSS brush. They noted counterion hopping between chains as the primary transport mechanism within the PEB. However, the relative importance of these competing factors (electrostatic interactions and steric hindrance) on counterion diffusivity as σ increases was not explored. In another study, Zhang 6 et al. again used FCS to study the effect of excess salt concentration on the diffusion of ionic fluorescent probes in two PEB systems: (1) a cationic rhodamine-6G derivative atop an anion PSS PEB and (2) an anionic rhodamine-6G derivative atop a cationic poly[(2-(methacryloyloxy)ethyl) trimethylammonium chloride] (PMETAC) PEB. Here, they found, in both systems, that as the salt concentration increased, osmotic pressure pushed counterions toward the substrate and deeper into the PEB at salt concentrations too low to introduce any detectable structural changes (quantified via thickness measurements) in the PEB.
They posited that the increased polymer concentration in the PEB closer to the substrate increased steric hindrance and reduced counterion diffusivity. Still, this study also did not directly investigate the effects of increasing σ on counterion diffusivity in PEBs, and neither of these studies looked at the effects of σ on polymer mobility in the PEB.
While simulation studies of PEBs have typically focused on the equilibrium PEB structure [7][8][9][10][11] and solvent transport 12 using generic coarse-grained (CG) and all atomistic (AA) molecular dynamics or self-consistent field theory, recent work has begun to investigate counterion dynamics. In particular, two recent molecular dynamics simulations were used to investigate counterion mobility as a function of σ. Yaun 13 et al. studied the effect of σ on counterion mobility in an electric field parallel to the substrate in a generic CG model of short PEBs (N = 10 monomers) in an implicit solvent over a range of σ. Yuan found non-monotonic behavior in counterion mobility with increasing σ where counterion mobility reached a maximum at intermediate σ. They attributed the trend in counterion mobility with σ to a balance between electrostatic interactions and steric hindrance. Namely, as σ increased, ions on the PEBs were pushed closer together lowering the energy barrier for counterions to hop from one ionic binding site to another on a different chain resulting in an increase in counterion mobility. However, as σ increased further, steric hindrance (between tightly packed polymer chains) dominates, resulting in a decrease in counterion mobility. Using AA molecular dynamics, Sachar 14 et al. investigated the effect of σ in fully charged poly (acrylic acid) (PAA) brushes (N = 36) on solvation structure and water/sodium counterion diffusivity. Their study concentrated on the effects of σ on water structure/dynamics. They found that the brushes induced "water-in-salt"-like scenarios that decreased the dielectric constant of the solvent and counterion diffusivity as σ increased. However, they did not investigate the relationship between σ and counterion diffusivity in the PEBs.
The experimental and simulation studies described above present some interesting findings but also leave questions to be answered. For instance, what is the effect of σ on counterion diffusivity in longer chain length (higher N) PEBs? Does the interplay between electrostatics and steric hindrance still lead to enhanced counterion diffusivity at intermediate σ in longer chain PEBs? And lastly, how does polymer mobility in longer chain PEBs affect counterion diffusivity?
Particle-based simulation methods are a powerful tool in answering these questions. However, while AA molecular dynamics studies provide the most detail, previous AA studies of PEBs were limited in the monomer number (N < 100) and the time scales (t < 400 ns) they could investigate.
On the other hand, generic CG models, capable of investigating longer PEB chains and time scales, lack the chemical specificity relevant to experiments. A better solution may lie between these two limits in a CG model that maintains chemical specificity relevant to the PEB, solvent, and counterions while allowing for investigations of dynamics at longer time and length scales.
To investigate these phenomena, we fine-tune a CG model of a strong polyelectrolyte that maintains the chemical specificity of poly[(2-(methacryloyloxy)ethyl) trimethylammonium chloride] (PMETAC) and explicitly models water and solvated chloride counterions using the MARTINI 15 forcefield. To better model electrostatics, we implement a refined version of polarizable MARTINI water 16 and polarizable MARTINI chloride counterions. 17 Using GPU accelerated GROMACS 18 simulations (version 2020.4), we map an AA simulation of a 100monomer PMETAC chain in solution, using previously published forcefield parameters, 7 to MARTINI beads using particle swarm optimization via Swarm-CG 19 to fine-tune the bonding parameters between mapped MARTINI PMETAC beads. Results of the MARTINI PMETAC model are verified by comparing the radius of gyration (Rg) and radial distribution functions (RDFs) from the simulations with AA reference simulations and experimental small angle X-ray scattering (SAXS) measurements 24 of the Rg of PMETAC chains of similar molecular weight in dilute solutions. We then build simulations of salt-free (only enough counterions to balance charge on the PEB) PMETAC brushes of N = 150 monomers at experimentally realizable σ of 0.05, 0.10, 0.20, and 0.40 chains/nm 2 , and investigate the brush structure, counterion diffusivity, ion dissociation dynamics, and polymer mobility. Results show that ion dissociation times are reduced at intermediate σ of 0.10 and 0.20 chains/nm 2 , but, after normalizing for excluded volume effects on counterion diffusivity within the brush, counterion diffusivity peaks at 0.10 chains/nm 2 , highlighting the competing effects of electrostatic interactions, steric hindrance, and polymer mobility on counterion diffusivity in PEBs.
Coarse-Graining PMETAC. The MARTINI 15 framework, previously used to coarse-grain similar methacrylate polymers 20 and polyelectrolytes, 21 was used to fine-tune a chemically specific CG model of PMETAC. Figure 1 shows how heavy atoms in the simulation were mapped to Lennard Jones (LJ)-type MARTINI CG beads based on their chemical constituents. representations of a PMETAC monomer (left, red), water (right-top, blue), chloride ions (rightmiddle, green), and silica substrate beads (right-bottom, grey). The MARTINI type is listed to the right of each CG bead. a All PMETAC beads are shown in red to match the color scheme in later figures, but each monomer consists of three distinct MARTINI beads to capture the chemical nature of the subgroups. b POL and PQa groups each contain two virtual charged beads with no LJ interactions to make them polarizable, DPOL and DPQa, respectively. They are not pictured here for clarity.
In each PMETAC monomer (red, figure 1, left), three backbone carbons were mapped to the SC1 MARTINI bead, a bead parameterized to replicate the apolar nature of hydrocarbons. The ester group was mapped to the Na MARTINI bead to represent a neutral chemical group with the ability to accept hydrogen bonds through its oxygen atoms. The quaternary amine group, which carries charge, was assigned to the Q0 MARTINI bead previously used to represent quaternary amine groups in choline 15 that carry positive charge but are incapable of donating or accepting hydrogen bonds. Water (blue, figure 1, right-top) was modeled using the refined polarizable POL MARTINI bead 16 to better capture screening effects in electrostatic interactions. Solvated chloride ions (green, figure 1, right-middle) were represented by a polarizable PQa MARTINI bead 17 representing chloride ions solvated by three water molecules. The polarizable water and chloride beads include two virtual charged beads, DPOL and DPQa (not pictured in figure 1), that allow the bead to polarize but have no LJ interaction parameters. To explicitly model a silica substrate in our simulations, the Nd MARTINI bead (gray, figure 1, right-bottom), a bead previously used to model silica 22 with neutral polarity and the ability to act as a donor in hydrogen bond interactions, was used. LJ bead interaction parameters, mass, and charge are summarized in Table S1.
Bonds between the CG beads in the PMETAC polymer were described by harmonic functions in eq. 1 and 2 𝑉 𝑏,𝑖 = 𝑘 𝑏,𝑖 (𝑙 𝑖 -𝑙 0,𝑖 ) 2
(1)
where kb,i is the bond length stiffness associated to bond 𝑖 with an equilibrium length li,0, and ka,i is the bond angle stiffness associated to bond i with an equilibrium angle of θi,0. LJ interactions between bonded beads and their two nearest bonded neighbors were excluded. The bonded parameters were optimized using a particle swarm optimization algorithm, Swarm-CG, 19 to best match the angle and bond distribution centers and breadths of all atomistic (AA) reference simulations of 100 monomer PMETAC chains in salt-free conditions using previously published AA parameters 7 using GROMACS (version 2020.2). Torsion angle potentials were not included in the MARTINI model of PMETAC because they prevented convergence of the optimization in Swarm-CG runs. LJ interactions between bonded beads and their two nearest bonded neighbors were neglected. The fine-tuned bonded parameters (Table 1) and details of the fine-tuning process, AA reference simulations, distribution overlap with the AA reference simulations, and CG-AA structural comparisons (radial distribution functions and radius of gyration) can be found in section S1. Salt-Free PMETAC Brush Construction and Simulation Protocols. We constructed brushes using MARTINI PMETAC chains with N=150 monomers (MW = 31.2 kg/mol) and studied experimentally realizable grafting densities of σ = 0.05, 0.10, 0.20, and 0.40 chains/nm 2 .
At these σ, all brushes simulated lie within the osmotic brush regime. 23 The substrate consisted of explicit substrate beads (MARTINI type Nd) assembled in an FCC lattice using the Atomic Simulation Environment 24 Python package with a lattice constant of 0.65 nm, slab thickness of 2 nm, and a length and width of 10 nm. Substrate beads were held in place by position restraints via a spherical harmonic potential with an elastic constant ksubs = 8000 kJ/mol-nm 2 following ref. 25 . Depending on the target σ, initiator beads of the same type as the MARTINI PMETAC backbone (SC1) were positioned at a height of 0.3 nm above the top substrate layer using an in-house, pseudorandom positioning algorithm that ensured the spacing between grafting points, including their periodic images, represented the average spacing expected at large length scales. Details about the pseudorandom positioning algorithm are detailed in section S2.
PMETAC monomers were stacked vertically from these initiator sites to create MARTINI PMETAC chains that extended from the initiator site in a fully stretched conformation. The z dimension of the simulation boxes was chosen to ensure that the box size was at least twice the height of the PEB after equilibration (~40 nm) and set to 105 nm. This height was chosen to minimize the effect of the explicit substrate and its periodic image on water structuring and to prevent the brush from interacting with its periodic image in the z direction. Using PACKMOL, 26 simulation boxes containing the brush and substrate were filled first by MARTINI chloride counterions confined to a region within the fully extended brush and subsequently by standard MARTINI water, fully defining the initial configuration of the simulation boxes for each σ.
After creating the box, all brush simulations contained ~300,000 beads and were carried out with the GPU-accelerated version of GROMACS 2020.4 at the Pittsburgh Supercomputing Center (PSC) on the Bridges-2 cluster. 27,28 We equilibrated the boxes as follows. Energy was minimized using steepest descent until the mean force reached a value lower than 10 kJ/mol-nm with periodic boundary conditions in the xyz directions. To expedite equilibration of the brush, boxes were annealed from 700 K (with velocities sampled from a Maxwell distribution) to 400 K in an NVT ensemble with a timestep of 20 fs over 125 ns using a v-rescale thermostat with a relaxation time of 1 ps. Following annealing, standard MARTINI waters were replaced with polarizable MARTINI water using Python scripts available at cgmartini.nl. Following the standard to polarizable water exchange, steepest descent energy minimization was again used to eliminate bead overlap until the mean force within the box reached a value lower than 10 kJ/mol-nm with periodic boundary conditions in the xyz directions. A short NVT equilibration of 50 ps initialized velocities after energy minimization to 400 K with a timestep of 5 fs. In an essential step, boxes were then compressed to remove vacuum bubbles that formed in the simulation using a semiisotropic NPzAT ensemble with a timestep of 10 fs for 10 ns to reduce the temperature from 400 K to 300 K using a v-rescale thermostat with a relaxation time of 1 ps, a Berendsen barostat with Analysis. The radius of gyration of the polymer chain, Rg, and the end-to-end distance, Re, were calculated at each timestep using GROMACS tools and with indexes for positions of the backbone carbon atoms. Snapshots of the simulation were captured using VMD. 29 Brush density profiles, mean square displacements (MSDs, eq. S3), and radial distribution functions (RDFs) were also calculated using GROMACS tools with relevant indexes. Diffusion coefficients (D) were determined by fitting eq. S4 to the MSD data in linear regions with the LMFIT 30 Python package.
Ion-pair correlations (𝐶(𝜏), eq. 3) were calculated using custom Python scripts and the MDAnalysis 31 Python package. The characteristic ion dissociation time (𝜏 𝐶 ) was calculated using eq. 3 using with parameters determined by fitting 𝐶(𝜏) to eq. 4 using the LMFIT 30 Python package.
(Eq. 3 and 4 are detailed in the Results Section). For both D and 𝜏 𝐶 , error bars were calculated as follows. The 5 µs production trajectories were binned into five independent 1 µs sub-trajectories.
MSDs and 𝐶(𝜏)s were calculated for each sub-trajectory, and those MSDs and 𝐶(𝜏)s were then fit to determine D normalized to by the free diffusivity (D0) and fit to equation 4 to determine 𝜏 𝐶 for each independent trajectory, respectively. The central point represents the mean of D and 𝜏 𝐶 from the five samples, and the error bars represent one standard deviation. Considering the simulations were periodic in the x and y directions, the pseudorandom positioning algorithm decreases the spacing between grafting points as σ increases. However, rather than uniformly spacing grafting sites on a grid, the random positioning of the grafting sites gave a closer approximation to the average spacing expected at experimental length scales. Thus, the pseudorandom positioning better represents the experimental conditions for growing polymer brushes from a surface.
Using this approach for positioning initiator sites, figure 2 (bottom) shows the brush conformations and the location of the water molecules and chloride ions after 5 µs of production time. As σ increased from 0.05 to 0.40 chains/nm -2 , chains stretched farther into the solution as the lateral spacing between the grafting sites decreased. At all σ, chloride ions (green) were localized within the brush to maintain local electroneutrality, with only a few ions transiently diffusing above the brush. We also note that, at the substrate, the first few monomers are adsorbed to the substrate, an effect that appears more pronounced at the lower σ. The density profiles shown in figure 3 quantify this observation and will be discussed later. σ after 5 µs of production simulation time. (top) Positions (xy) of the grafted backbone initiator beads (red) for the PEBs placed using the pseudorandom algorithm. All boxes are 10 nm x 10 nm. (bottom) Side-on images of the brushes in water as σ increases. Gray beads are position-restrained substrate beads, red beads make up the MARTINI PMETAC chains, and green beads are chloride ions. Water beads are shown transparently in teal to facilitate visualization. Images are 10 nm x 41.4 nm. Beads are not drawn to scale. The time-averaged density profiles of each simulation for PMETAC (left, red), chloride ions (middle, green), and water (right, blue) are shown in figure 3. For all cases, the profiles are perpendicular to the substrate in the z-direction (see figure 2). In agreement with the snapshots direction (perpendicular to the substrate) of PMETAC (left, red), chloride ions (middle, green), water (right, blue) for σ values from 0.05 to 0.40 nm -2 (increasing with darkness) of 150 monomer PMETAC MARTINI PEBs in salt-free conditions. The first 2 nm of the profile is zero because this is the substrate.
The initial peak in the PMETAC density profiles at the substrate shows the first few PMETAC monomers flattened against the substrate. The flattening of the first few monomers in polymer brush simulations has previously been seen in CG simulations of polymer brushes 25,32 and LJ simulations of liquids in contact with solid surfaces. 33 In these studies, the substrates and polymer had the same enthalpic affinity for the solvent, and the authors attributed the adsorption of chains to entropic effects that relieved the structuring of the solvent at the interface. In our study, the situation is different. The substrate beads and PMETAC backbone beads are less polar than water and are modeled with different interaction strengths with the solvent and substrate (εSC1-POL = 1.2 kJ/mol and εSC1-Nd = 2.7 kJ/mol). This difference in interaction strength causes the backbone beads to enthalpically interact more favorably with the substrate than the polarizable MARTINI water. Because the hydrophobic nature of the backbone reflects the chemical nature of the PEB system, flattening may reflect the preferential adsorption of polymer chains to the substrate in PEB systems that may be present in natural systems but are difficult to measure experimentally.
At all values of 𝜎 simulated, the salt-free PEBs lie well within their osmotic brush regime (σ = 0.0001 nm -2 to σ = 1.0 nm -2 ). 23 In the osmotic regime, counterions are strongly confined to the brush and PEBs, and PEB height is driven by the osmotic pressure of counterions within the brush. Here, PEB height scaling is predicted by Zhulina and Rubenstein 23 to be σ 0 . The PMETAC density profiles in our simulations also do not follow the parabolic profiles derived and measured for neutral polymer brushes at high σ. 34 Our profiles are box-like profiles that display exponential tails previously observed in other MD simulations of PEBs. 10,11 We also observe weak scaling of PEB height with σ. Following He 10 et al., we estimated the height of the brushes from the density profiles by defining the height of the PEB as the height at which the PMETAC density decreased to 10% of the average density within the brush. The brush heights (or distance above the 2 nm substrate) were calculated to be 30.3, 31.0, 35.2, and 37.3 nm for σ of 0.05, 0.10, 0.20, and 0.40 chains/nm -2 , respectively. The weak scaling of brush height with σ for PEBs in the osmotic regime, again seen in previous PEB simulations, 35 indicates other factors captured in these simulations such as excluded volume effects of the PEB on counterion entropy, osmotic pressure, and finite extensibility also contribute to PEB structure. Most importantly, figures 2 and 3 importantly show that counterions are confined to the PEB and support the conclusion that these PEBs are within the osmotic regime, a detail used to account for excluded volume effects on counterion diffusivity in the next section.
Counterion Diffusivity. In section S4, the ensemble average mean squared displacements (MSD) and diffusion coefficients (D) of chloride counterions (Cl -) in the PEBs were calculated at each σ. Figure S11 shows that prior to the onset of confinement from the finite PEB height, D in the lateral and axial directions was comparable, and we therefore focus on Dxy in this work. As σ increased, figure S11 shows the Dxy of the counterions decreased monotonically. However, one limitation of the MARTINI coarse-graining is that Cl -ions incorporate the ion and part of its solvation shell into a single bead, which retains its partial solvation shell throughout the simulation.
In developing polarizable ions for the MARITNI framework, excluded volume effects from the coarse-grained solvation shell were observed to decrease Cl -diffusivity in simulations of aqueous solutions with increasing salt concentration. 17 This reduction in diffusivity with increasing salt concentration is not an effect seen in experiments at the salt concentrations considered in this study. Thus, to normalize for the effects of excluded volume on diffusivity due to coarse-graining, we referenced the free diffusivity, D0, of Cl -in aqueous solutions of polarizable MARTINI Na + and Cl -ions at the same salt concentration as the Cl -ions confined to the PEB volume environment.
Because all Cl -ions in the simulation were confined to the brush, we estimated the salt concentration in a box of volume Vbrush = 10 nm x 10 nm x H(σ) where H(σ) is the height of the brush for each σ. The concentration of Cl -within the brush was then estimated as 𝐶 𝑏𝑟𝑢𝑠ℎ,𝐶𝑙 -(𝜎) = was acquired from ref. 17. Figure 4 shows Dxy/D0 values for each σ. σ. D0 corresponds to ion diffusion in an aqueous solution of equivalent ion concentration as that within the brush. The center of the points represents the mean of Dxy/D0 for independent 1 µs production trajectories and error bars represent one standard deviation. For all σ except 0.40 nm -2 , error bars are smaller than the plotted point. S12). First, at all σ, electrostatic interactions and steric hindrance in the brush suppressed the relative Cl -diffusivity to a fraction of 0.40 or less compared to the freely diffusing Cl -. This suppressed diffusivity is consistent with results from recent all atomistic simulations of counterion diffusivity in strongly charged polyelectrolyte brushes. 36 Second, there is a non-monotonic trend in diffusivity with increasing σ.
At low σ = 0.05 nm -2 , Dxy/D0 is slightly reduced with a value of 0.36. The maximum Dxy/D0 was observed at an intermediate σ = 0.10 nm -2 with a value of 0.40, then decreased to 0.31 at σ = 0.20 nm -2 and 0.19 at σ = 0.40 nm -2 . A non-monotonic dependence in counterion mobility with σ in a strong PEB was previously reported by Yuan 13 et al. in simulations containing 10 generic coarsegrained monomers. In that study, ion mobility peaked at a σ of 0.20 nm -2 after converting their LJ units to the real units used in our simulation. They attributed this peak in ion mobility to a competition between electrostatic interactions between counterions and ions on the PEB and steric hindrance from increased polymer density in the PEBs at higher σ. In this study, with longer polymer chains (150 monomers), the peak in mobility occurs at a lower σ of 0.10 nm -2 . We investigate the origin of this difference in the following sections.
Ion Pair Dissociation. To investigate the nonmonotonic trend in Dxy/D0 with σ, we characterized ion dissociations in the brush by calculating the structural relaxation time, 37,38 𝜏 𝐶 , of dissociations between the cationic amine groups on the MARTINI PMETAC chain and anionic chloride counterions. We calculated the intermittent ion-association autocorrelation function, C(τ), to track whether ion pairs associated at time t0 were still associated at a time t0 + τ using eq. 3
where p(t0) is unity when ions are associated at a distance less than the cutoff distance, rcut, p(t0 + τ) is unity if ion pairs are associated at a distance less than rcut some time t0 + τ and zero otherwise, and p is the number of associations at t0. We set rcut = 0.652 nm, corresponding to the position of the first minimum in the amine-chloride ion RDF plot in figure S13.
Following Bollinger 38 et al., C(τ) was then fit to a stretched exponential function (eq. 4), and the characteristic ion dissociation relaxation time, 𝜏 𝐶 , was calculated using eq. 5.
C(τ) for each 𝜎 is shown in figure S14 and 𝜏 𝐶 verses 𝜎 is shown in Figure 5. PMETAC monomer. (b) The characteristic ion dissociation structural relaxation time (τC) as a function of σ. The center of the points represents the mean of τC for independent 1 µs production trajectories, and error bars represent one standard deviation. For σ of 0.10 and 0.20 nm -2 , error bars are smaller than the plotted point.
Figure 5 shows that the ion-dissociation relaxation time depends on grafting density.
Initially, τC decreases by about 44% as σ increases and then increases by 124% at the maximum σ. its original ionic bonding pair. In other words, τC quantifies the ionic bond dissociation time. Thus, smaller values of τC indicate faster dissociations of the chloride counterions from their amine
partner, which, in turn, would enhance diffusivity. At the lowest σ = 0.05 nm -2 , chloride anions remained associated with their initial cationic partner for a τC of 1544 ps. At intermediate grafting densities of σ = 0.10 nm -2 and σ = 0.20 nm -2 , τC decreased to 851 and 891 ps, respectively. This behavior is consistent with faster dissociations compared to the low σ case. Finally, at the highest σ = 0.40 nm -2 , τC increased to its highest value of 1996 ps, representing the slowest dissociation time.
Because τC quantifies the intermittent ion-dissociation time, it captures effects from the electrostatic interactions and the effect of the local environment (steric hindrance) and its relaxation (i.e., the polymer and water relaxation) on ion pair dissociations. Thus, to interpret the trend in τC, we follow Yuan 13 et al., who described counterion mobility in PEBs as a tradeoff between electrostatic interactions and steric hindrance and consider the electrostatic energy landscape a counterion experiences in a PEB. Here, the polymer chains are approximated to be confined within a cylinder with a radius equal to the length of their side chains centered on a grafting site. For a counterion to hop from one chain to the next, Yuan 13 et al. postulated that counterion must overcome an electrostatic barrier 𝑉 𝑏𝑎𝑟𝑟𝑖𝑒𝑟 = 〈𝑉 𝑐𝑜𝑢𝑙 (𝑟 𝑚𝑎𝑥 ) -𝑉 𝑐𝑜𝑢𝑙 (𝑟 𝑚𝑖𝑛 ) 〉 to dissociate from their current bonding pair, where 〈𝑉 𝑐𝑜𝑢𝑙 (𝑟)〉 = 𝑞 𝐶𝑙 -𝑞 𝐴 + /4𝜋𝜀 0 𝜀 𝑟 𝑟 is the ensembleaveraged electrostatic energy of the counterion. rmax is the spacing between a charged group on a polymer side chain and the nearest charged group on a different grafted chain in the PEB. rmin is the minimum separation of the associated charged pairs in an energy well whose value is determined by the first peak position of the chloride-amine RDF (0.445 nm) in figure S13. We can approximate rmax(σ) by calculating the spacing between edges of cylinders as 𝑟 𝑚𝑎𝑥 (𝜎) =
-𝑙 𝑠𝑖𝑑𝑒 𝑐ℎ𝑎𝑖𝑛 ) where lside chain = 0.63 nm and probe the effect of σ on Vbarrier in figure 6. cylinder edges approximating the positions of side chains of fully stretched grafted MARTINI PMETAC chains used to calculate Vbarrier. (bottom) The electrostatic barrier (Vbarrier) to ion hopping between cylinder edges that approximates the average position of charge groups on fully stretched grafted chains as a function of σ. The dashed black line indicates where thermal energy is equivalent to the electrostatic barrier.
For chains confined to vertical cylinders at their grafting points, figure 6 (bottom) shows that Vbarrier decreases as σ increases and dips below thermal energy at grafting densities greater than 0.20 nm -2 . At the lowest σ = 0.05 nm -2 , the large spacing between chains increases Vbarrier, and the energy required for counterions to hop from one chain to another is greater than kbT in the system. This lowering of Vbarrier makes it less likely for ions to disassociate from their initial association site and leads to high τC (long dissociation times), limiting the ability of diffusing ions to hop from chain to chain. At intermediate grafting densities, Vbarrier approaches kbT. However, polymer brush, longer chains in the PEB can undergo greater lateral excursions and, therefore, interact with more neighboring chains. These increased lateral excursions suggest that the mobility of the polymer, and thus the ability of the neighboring chains to move closer together to lower the electrostatic barrier to ion dissociation transiently, may enhance the hopping of counterions from chain to chain. Furthermore, when ions are associated in ionic bonds to a given chain, they will move with the chain as it moves. In that way, the increased lateral excursions of the polymer chains could also directly contribute to the diffusivity of ions as the chains carry ions with them.
To investigate the mobility of the polymer chains in our simulations, figure 7 shows the lateral ensemble average mean square displacement (MSDxy) of polymer backbone beads (SC1) in the MARTINI PMETAC brush simulation at each σ calculated using eq. S3. backbone beads (SC1). σ increases with the darkness of the line, as indicated by the legend.
The mean square displacement of the backbone beads exhibits two main trends. First, as σ increases, the MSDxy decreases in magnitude at all lag times, indicating that neighboring chains sterically limited the lateral mobility of the backbone beads monotonically as σ increased. Second, when considering the slopes of the MSDxy curves, the slope of the curves at σ = 0.05 and 0.10 nm - 2 remain constant from 𝜏 equal to 10 -2 to 10 0 µs. In contrast, the slopes of the σ = 0.20 and 0.40 nm -2 MSDxy curves decrease at 𝜏~10 -1 µs, indicating slower dynamics at long 𝜏. This slope change also interestingly coincides with the decrease in diffusivity of the counterions at grafting densities greater than or equal to 0.20 nm -2 in our simulations. The enhanced reduction in polymer chain dynamics may indicate a change in the segmental motion mechanism from Zimm to Rouse as the concentration of polymer in the brush increases, but a complete and systematic comparison to the polymer scaling physics is needed in future studies.
When considering the effects of polymer chain dynamics with ion dissociation times, we see the complex effects of polymer chain dynamics on ion diffusivity. At the lowest grafting density, mobile chains carry ions with them, but the spacing between chains limits interchain ion hopping. At intermediate grafting densities, polymer dynamics begin to slow, but the proximity of the chains to one another decreases, lowering the electrostatic activation barrier and facilitating hopping of ions between chains. These competing effects lead to increased ion diffusivity when relatively mobile chains are spaced closely enough to lower the electrostatic barrier to interchain ion hopping and then decreased ion diffusivity when the polymer chain mobility is significantly hindered. Lastly, at high grafting densities, steric hindrance and slowed polymer chain dynamics dominate, and ion diffusivity is slow despite the low activation barrier. As grafting density increases, there is a complex interplay between reducing the electrostatic barrier to enhance hopping, mobile polymer chains transporting associated ions, and increased steric hindrance. The maximum value of counterion diffusivity at intermediate σ (figure 4) likely results from a compromise between ion hopping and coupling with the polymer chain dynamics. Counterions that dissociate from their initial ionic bonding site do not diffuse but instead rattle in their local environment. Specifically, the counterions are confined by the concentration of polymer in the brush and undergo reverse hopping to and from their initial ionic association site. We note that while this study investigates the effect of σ on counterion diffusivity, it would be interesting to investigate these mechanisms in future studies by probing the effects of increasing molecular weight and polydispersity on counterion diffusivity in polyelectrolyte brushes.
In this work, we fine-tuned a coarse grained (CG) model of a strong polyelectrolyte (poly[(2-(methacryloyloxy)ethyl) trimethylammonium chloride]) [PMETAC]) that retained chemical information about the components of the system and matched the local structure and radius of gyration in all-atom (AA) simulations and experiments using the polarizable MARTINI forcefield. We then built simulations of salt-free PMETAC brushes with 150 monomers at experimentally realizable grafting densities, σ, of 0.05 chains/nm 2 , 0.10 chains/nm 2 , 0.20 chains/nm 2 , and 0.40 chains/nm 2 and investigated the brush and counterion structure, ion dissociation dynamics, and diffusivity. Results showed that ions dissociate faster at intermediate σ of 0.10 chains/nm 2 and 0.20 chains/nm 2 , but, after normalizing for excluded volume effects on counterion diffusivity within the brush, diffusivity reaches a maximum at 0.10 chains/nm 2 . The peak in diffusivity at intermediate σ highlights the competing effects between electrostatic interactions, polymer mobility, and steric hindrance on counterion diffusivity in polyelectrolyte brushes.
displacements, non-normalized ion diffusivity, radial distribution functions, and the intermittent ion-association autocorrelation function.
Portions of this work were part of a PhD thesis deposit with DOI:29395821.