Introduction

Fuel cell deployable proton exchange membranes (PEMs) have been studied extensively over the past few decades due to their promise in technologies for clean and efficient power generation. [1][2][3][4][5][6][7][8][9][10][11][12] In recent years, nano-conned environments have been exploited in the study of cost-effective and reliable polymer architectures for electrochemical devices. [13][14][15][16][17][18][19][20][21] Understanding the behavior of water and ions in these conned environments is essential to gain insight into ion diffusion mechanisms within these devices.

The three structural facets that must be considered in the design of new PEMs under nano-connement are (a) the polymer backbone and associated mesoscale morphology, (b) an anionic group, and (c) a tether that connects the anionic group to the polymer. The sulfonate anionic functional end group (SO 3

À

) is one of the most widely used groups in PEM fuel cell devices. [22][23][24][25][26][27][28][29][30][31] The protonation state of SO 3

À and its dependence on water content and temperature in the system have been studied extensively using techniques such as NMR, [32][33][34] vibrational spectroscopy, 35,36 Raman spectroscopy 37 and X-ray scattering. 38,39 Specically, properties such as PEM morphology and the local structure of the sulfonate-terminated side chains have received increased attention both experimentally [40][41][42][43][44] and theoretically. 40,42,[45][46][47][48] Several studies suggest that the sulfonate end group appears in a deprotonated state (i.e., SO 3 À ) in hydrated systems, or in a protonated state (i.e., SO 3 H) at low-hydration conditions. In some cases, it was argued that the hydronium diffusion mechanism is mainly vehicular under low hydration conditions, whereas there is a signicant contribution from structural diffusion at high hydration conditions. 28,42,[49][50][51] Despite the abundance of studies in the eld of PEMs, ongoing discussion continues about the role of the protonation state of SO 3 À in the underlying atomistic picture governing the hydronium ion diffusion process. At a more fundamental level, elementary steps governing proton transport phenomena in hydrogen-bonded media continue to be of considerable interest to the physical chemistry community. Understanding and controlling the proton transport mechanism on the atomistic level is essential in order to develop PEMs with high hydronium conductivity. The computational methods most frequently used to study the SO 3 À end groups in PEMs include mesoscopic simulations, 40,42,45,46,48,[51][52][53][54][55] classical molecular dynamics (MD) simulations, 4,49,50,[56][57][58][59] and density functional theory (DFT) calculations. 5,22,52,60 On a more fundamental level, although some ab initio molecular dynamics (AIMD) simulations 28,[61][62][63][64] have been performed over the past decade, the high computational burden of such an approach has limited its use in this area. Nevertheless, the use of AIMD, in which the interatomic forces are computed "on the y" from DFT based electronic structure calculations as the simulation proceeds, is necessary when studying aqueous hydronium and hydroxide diffusion. As was shown in ref. 65, AIMD predicted a fourth weak hydrogen bond to the hydronium ion in the proton transfer (PT) process. This result was later used to parameterize new multistate empirical valence bond models. 66 Such empirical models cannot be easily transferred to different chemical environments such as those investigated here. The trade-offs in the use of AIMD over a reactive force eld are the usual limitations in accessible length and time scales. These limitations require careful selection of the systems to be studied. Recently, we used fully atomistic AIMD simulations to study hydroxide diffusion in model anion exchange membranes (AEMs) using nano-conned environments consisting of graphane bilayers (GBs) as mimics of the actual polymer architectures. [67][68][69][70] We nd that factors controlling hydroxide diffusion in these systems, including local coordination patterns and pre-solvation mechanisms, were shown to differ from those in bulk solution [71][72][73][74][75][76][77][78][79][80][81] in a way that is strongly inuenced by the shape and size of the conning volume, the hydration level, and the cation spacing.

In this study, we apply a similar protocol to explore hydronium diffusion in two architecturally distinct PEMs, employing nano-conned volumes inspired from our previous studies. [67][68][69][70] The choice of GBs to mimic a particular polymer architecture with a layered morphology in the study of PEMs was inspired by a recent study of Trigg, et al., 7 in which it was shown that wellordered and hydrated membranes with highly crystalline morphology have the potential to achieve high proton conductivity. We nd that the protonation state of SO 3

À changes during the course of the simulation, as under specic conditions, the following reaction occurs: SO À 3 + H 3 O + 4 SO 3 H + H 2 O. This study aims to uncover both the conditions required for this reaction to occur and the role of this reaction in the hydronium diffusion mechanism. In addition, we present a comparison of our ndings to previously studied hydroxide diffusion mechanisms in analogous AEM environments. We believe that identifying the differences in the diffusion mechanisms of ions in PEMs and AEMs will help reveal key principles for new stable membrane materials with high proton conductivity.

Description of systems

In this study, we explored two different GB systems representing two model PEM environments, in which the specic nano-conned structure is designed to mimic the layered arrangement recently reported in ref. 7. Each system contains two identical graphane layers parallel to the xy-plane, two SO 3 À anions attached to a graphane sheet using a (CH 2 ) 2 linker, two hydronium ions (whose oxygen cores are denoted O * 1 and O * 2 ), and a variable number of water molecules. The two anions are attached by the linkers to xed points on a graphane sheet but are otherwise free to move in the aqueous solution. The two attachment points dene the polymer electrolyte anion spacing in the x and y directions (see Fig. 1). Based on ref. 7 and 82, the tunable parameters for the two systems are: (i) the hydration level, l, chosen to be 3 or 4, (ii) the distance between the two carbon sheets, Dz, xed at 7.3 Å for all systems (see ref. 67 and 68 for rationale), and (iii) the polymer electrolyte cation spacing in the x and y directions, as measured between two sulfur atoms (Dx and Dy), in which Dx and Dy are xed at 10 Å and 6.6 Å, respectively, for the two systems. For clarity, we refer to the two systems as l3 and l4, in which the numbers represent the respective hydration levels.

Computational method

Once the desired starting structures were generated, AIMD simulations 83 were performed using the CPMD code. 84,85 Each system was equilibrated at room temperature using a massive Nosé-Hoover chain thermostat, 86 followed by 15-20 ps of canonical (NVT) dynamics, also using a massive Nosé-Hoover chain thermostat, and nally $80 ps of microcanonical (NVE) dynamics. In order to account for dispersion forces, we employed the Dispersion-Corrected Atomic Core Pseudopotentials (DCACP) scheme 87,88 within the Kohn-Sham formulation of Density Functional Theory using the B-LYP exchangecorrelation functional. 89,90 The performance of B-LYP + DCACP has previously been shown to give satisfactory results for wateracene interactions, 91 for liquid water, 92 and for hydronium diffusion in bulk water. 65,[71][72][73] A detailed description of the computational method can be found in our previous work. [67][68][69][70] 4. Results

Solvation structures

The primary aim of this work is gain insight into how changing the water content affects hydronium diffusion in the model systems employed. To this end, we begin by exploring the solvation structures of the water molecules, the hydronium ions, and the sulfonic acid end groups.

4.1.1 Water structure. Fig. 2 shows the spatial populations of oxygen atoms in the xy-plane generated from the trajectory. This allows us to glean the preferred locations of water molecules and provides a clear picture for the water density prole in the xy-plane. The results show that the water distribution in the cell of system l3 is not uniform. However, for system l4, which has one extra water molecule per cation than does l3, the water distribution was found to be uniform. Inspection of the congurations from the AIMD trajectories support these ndings (see Fig. 2). Specically, for system l3, we nd that void areas are formed in parts of the simulation cell. At this low level of hydration, all waters in the system can be regarded as interfacial, that is, in contact with some part of the "membrane", and inhomogeneously distributed throughout the system. Furthermore, we nd that the non-uniform/uniform water distribution for systems l3 and l4, respectively, persists throughout the simulation.

Unlike in bulk solution, in which the water oxygen has, on average, a fourfold-tetrahedral coordination pattern, the low hydration values in the two systems result in a rst solvation shell of approximately three for the water oxygens, in which two are water oxygens and one is an anion oxygen, as seen in the OO radial distribution function (RDF) presented in the ESI. † The differences in the water distribution between the two systems are pronounced in the second solvation shell of the water oxygens, for example, the integrated coordination numbers (CNs) of the second solvation shell are approximately 4 and 5 for systems l3 and l4, respectively.

These differences in water distribution at low hydration were previously seen in our recent work on low hydrated AEMs. 68 Furthermore, it was shown that the hydroxide ion diffusion is vehicular for non-uniform water distributions, and structural for uniform water distributions. Similar to the AEMs, we nd that for PEMs, the water distribution affects the hydronium ion diffusion. However, as will be discussed in the next sections, the hydronium diffusion mechanism in PEMs is fundamentally different from the hydroxide diffusion mechanism found for AEMs.

4.1.2 H 3 O + solvation structure. We turn next to an exploration of the hydronium ion solvation structure. For this purpose, we plot, in Fig. 3, the O*O RDF and CNs (O* represents the hydronium oxygens, and O represents SO 3 À and water oxygens). As shown, the rst solvation shell peak is located at 2.6 A for both systems, and the CN values of the rst solvation shell of the two systems is 3.1, as was previously found for bulk solution. 71 Furthermore, we nd that the oxygens taking part in the hydronium solvation complex are both water and SO 3 À oxygens (see inset of Fig. 3 for examples). To support these results, we calculated the population probabilities for the H 3 O + solvation complexes; these populations indicate that the most likely complex is 3A + 0D with 90% and 86% for systems l3 and l4, respectively (see hydrogen bond (HB) criteria in ESI †). While the rst solvation shell of the hydronium ion is approximately identical in the two systems, we nd that the difference between the hydronium solvation structures is pronounced in the second solvation shell, with CN values of 7.5 and 8.5 for systems l3 and l4, respectively (including the SO 3 À oxygens). Excluding the rst solvation shell oxygens and the SO 3

À oxygens, the numbers of water oxygens in the second solvation shell of the hydronium ions are 1.5 and 2.5 for systems l3 and l4, respectively. This suggests that the hydronium ions in system l3 are missing a complete second solvation shell as a result of the non-uniform water distribution that develops in this system (to support these results, the O*O w RDF and CNs are presented in the ESI †). 4.1.3 SO 3 À solvation structure. To explore the solvation structure of SO 3 À oxygens, we plot the OsO RDF and CNs for the two systems in Fig. 4 (Os represents the SO 3 À oxygens, and O represents all other oxygens in the system). As shown, the rst and second peaks, which represent the rst and second solvation shells, are located at approximately the same value of r for the two systems. A comparison of the CN values (Fig. 4b) shows elevated numbers for system l4 for both the rst and second solvation shells with values of 1.8 and 4.1, respectively, compared to values of 1.5 and 3.3 for l3. This result, which is a direct outcome of the higher hydration values in system l4, shows that the second solvation shell of the SO 3 À oxygens in system l4 contains one extra water oxygen compared to system l3.

H 3 O + diffusion mechanism

In order to shed light on the transport process of hydronium ions in this conned environment, we calculate water and hydronium diffusion constants along each of the axes separately (see Table 1). These components can be interpreted as the diagonal elements of the diffusion tensor, an important quantity in the calculation of ionic conductivities. [67][68][69][70] In the ESI, † we present a simple time trace of the coordinates of the hydronium oxygens along the trajectory.

A comparison of the diffusion constants of the two systems shows that increasing l from 3 to 4 results in a decrease in the average hydronium ion diffusion constant (0.070 Å2 ps À1 and 0.007 Å2 ps À1 for systems l3 and l4, respectively). Specically, system l4 is seen to have the lower hydronium diffusion constants along all axes. This agrees with the evolution of the coordinates of the two hydronium ions (ESI Fig. S2 †), where both ions are non-diffusive. In system l3, hydronium ion diffusion occurs along the x-and y-axes, with diffusion tensor components of 0.106 Å2 ps À1 and 0.095 Å2 ps À1 along each of these directions, respectively. According to Fig. S2 in the ESI, † both ions in system l3 become diffusion "activated" (at $50 ps) aer an initial quiescent period. The water molecules were found to be non-diffusive for both systems, suggesting that water diffusion is not required for hydronium ion diffusion (the non-diffusivity of the water molecules accords with the water prole presented in Fig. 2). Further insight into the conditions that enable the diffusion of the two hydronium ions in system l3 and suppress their diffusion in system l4 can be gleaned by turning to the RDF and CNs for SO* presented in Fig. 5a. As shown in the gure, the positions of the rst two peaks, which represent the rst and second solvation shells, are identical for both systems. The rst peak, located at $1.7 Å, corresponds to SO 3 H, in which an H 3 O + has transferred a proton to an SO 3 À , while the second peak is located at $3.7 Å. The SO* CN values (specied in Fig. 5b) for the rst peak are 0.7 and 0.02, respectively, and for the second peak, they are 0.8 and 0.4 for systems l3 and l4, respectively (note that, because of the limited system size and small number of S and O* species, the SO* CN does not go to 1 at large r, as it would in a larger system with a greater number of these atom types). This suggests that due to the reaction of hydronium with SO 3

À anion to create SO 3 H, the latter neutral species exists for a greater proportion of time in the l3 system than it does in l4.

To verify this, we calculate the percentage of time that the hydronium ions spent as H 3 O + $SO 3 À and H 2 O$SO 3 H (Fig. 5c).

As expected, we nd that for system l3, the hydronium ion appears as H 2 O$SO 3 H for 48.78% of the simulation time, while for system l4, it appears as H 2 O$SO 3 H for only 3.6% of the simulation time. Our results support those of a recent classical molecular dynamics study reported by Sengupta et al., 55 who observed that the hydronium ion diffusion decreases with increasing degree of deprotonation. It is important to note that the use of classical molecular dynamics requires xing the protonation state a priori, whereas it can vary naturally in the present AIMD simulations, which means that a direct comparison with the results of ref. 55 is not possible. It, nevertheless, appears that the protonation state of the sulfonate group plays an important role in the hydronium diffusion process at low hydration. Table 1 Diffusion constants obtained from the slope of the mean square displacement as presented in ESI (all in units of Å2 ps À1 )

System l3 0.070 0.106 0.095 0.010 0.039 0.083 0.036 0.001 System l4 0.007 0.003 0.011 0.006 0.009 0.006 0.004 0.017 Beyond the protonation state, however, we conclude, more specically, that the reaction SO 3

plays an important role in the hydronium ion diffusion mechanism in system l3. In order to shed additional light on the conditions that enable this reaction, we shi our focus to O next , where O next is the closest water or hydronium oxygen to the SO 3 À oxygens (see, also, ref. 65). In Fig. 6a we present the O next O RDF and CNs for the two systems. The rst peak is located at 2.7 Å and 2.6 Å for systems l3 and l4, respectively (see further discussion in ESI †). The CN values for the rst and second peaks are 1.3 and 5.1 for system l3 and 1.48 and 6.0 for system l4. The CN values found for system l3 suggest that before PT occurs between the protonated anion (i.e., SO 3 H) and a nascent water molecule, the latter acquires a rst solvation shell of one water oxygen and an incomplete second solvation shell.

In order to garner additional support for this claim, we investigate system l3, where the reaction is mainly observed, in greater depth. We dene a displacement coordinate, d ¼ |R OaH À R OwH |, where R OaH and R OwH are the distances between a shared proton of SO 3 H and the nearest water oxygen (i.e., O next ). Values of d > 0.5 are considered to be inactive complexes with respect to PT, while values of d < 0.1 are considered to be "active" and are associated with PT events. 67,74,75,93,94 In Fig. 6b, we present the O next O RDF and CNs for system l3 for d < 0.1 and d > 0.5, where, in this context, O next represents the rst neighbor oxygen to the SO 3 H oxygen, and O represents water and hydronium oxygens. We nd that for d > 0.5, the peak is located at 2.8 Å with a CN value of 1.54. However, for d < 0.1, which is associated with a PT event, the peak is located at 2.7 Å with a CN value of 1.14. This suggests that in order for the reaction to take place, O next is required to have a CN value of $1. Combining the results presented in Section 4, we conclude that the higher reactivity of the hydronium ion seen in system l3 can be explained in terms of the non-uniform water distribution in the system, which results in: (i) a higher probability of obtaining a CN value of $1 for O next , (ii) an incomplete second solvation shell for the hydronium ions, and (iii) fewer water molecules in the vicinity of the oxygen atoms of the anions (i.e., SO 3 À ).

Based on the results above, combined with inspection of congurations from the trajectory, we propose, in Fig. 7, an idealized diffusion mechanism for hydronium ions in PEMs under idealized hydration conditions (i.e., system l3). First, three water molecules in the simulation box solvate the hydronium ion, which is located in the center of the cell (Fig. 7a). Next, a PT occurs from the hydronium ion to a nearby water molecule (Fig. 7b). A hydrogen bond (HB) is formed between the nascent hydronium ion oxygen and the anion, SO 3 À , while the hydronium has only one water oxygen in its rst solvation shell (Fig. 7c). Finally, a PT occurs between the hydronium ion and the anion (i.e., SO 3 À + H 3 O + ), resulting in SO 3 H + H 2 O (Fig. 7d).

This procedure cycles back to the initial condition and restarts. The next PT occurs once SO 3 H donates its hydrogen to a nearby water molecule with a rst solvation shell consisting of only one water oxygen (Fig. 7f), which results, again, in SO 3 À + H 3 O + (see details in Fig. 7e-h).

Discussions and conclusions

The recent work of Trigg, et al. 7 reported a combined experimental and theoretical study demonstrating that a hydrated layered membrane can promote proton conductivity on par with current benchmark materials such as Naon. An innovative and precise polymer design controls polymer folding to achieve a well-ordered layered membrane. X-ray scattering, transmission electron microscopy and all-atom classical molecular dynamics simulations were used to reveal and characterize this and related layered morphologies, and electrochemical impedance spectroscopy probed the conductivity of the sulfonated polyethylene, p21SA, membrane. In addition, classical MD simulations were performed to explore the dynamics of water and hydronium ions in this particular material. For this purpose, chain folded molecules of atactic p21SA arranged in layers were studied under hydration values of 3 # l # 5.5. The MD simulations suggested that the hydronium ions are close to both the sulfonate groups and to water molecules. Furthermore, the simulations indicated that for a system with l ¼ 3, hydronium ions are coordinated by zero or one water molecules, while for l ¼ 5.5, the hydronium coordination environments have either two SO 3 À groups and two waters or one SO 3 À group and three waters (see Fig. 3 in ref. 7).

Comparison of the structural characterization obtained from the classical MD results reported in ref. 7 with our AIMD graphane bilayer models conrm the picture of the hydronium ion coordination structure in these well-ordered layered membranes. Specically, both studies show that the hydroniums are located in the vicinity of both the anions and the waters. Secondly, and more importantly, both studies nd that the hydronium ions in systems with l ¼ 3 are more likely to be coordinated by only one water molecule, while for higher hydration values the hydronium ion is more likely to be found with a full threefold coordination shell.

While classical MD simulations enable investigation of relatively large systems and are consider to be a useful tool for structural characterization, their use of empirical force elds with xed charges means that they cannot capture chemical bond breaking and forming events, nor do they include manybody polarization, both of which are critical for describing the hydronium structural diffusion process. 95 This limitation motivated our use of AIMD, which allow us to elucidate the atomistic/mechanistic details of the structural diffusion process that drives hydronium ion transport in PEM materials.

Specically, the AIMD simulation results presented here suggest, perhaps somewhat counterintuitively, that the low coordination state of the hydronium ions in system l3 is actually critical for achieving high hydronium conductivity, as the low water content ultimately promotes the reaction: SO

, which is properly captured in the AIMD simulations and which we nd to be an essential part in the hydronium diffusion mechanism in the l3 system. We believe that the discovery of such ion diffusion mechanisms has broad implications on future characterization of new stable polymer electrolyte materials with high ion conductivity.

In order to demonstrate the importance of revealing key principles in the hydronium diffusion mechanism, we nd it useful to compare our ndings to previously studied hydroxide diffusion mechanisms in analogous AEM environments. While both AEMs and PEMs have been studied extensively over the last decade, it is well known that hydroxide ion conductivity and cation stability remain key hurdles to realizing the full potential of fuel cell based AEMs. In our previous study on hydroxide ion diffusion in model AEMs under low hydration conditions (2 < l < 5), 67,68 we found that the water distribution is non-uniform. Comparing to the present study, we nd that while the non-uniform water distribution is a common feature of both AEMs and PEMs at low hydration, the inuence of the non-uniform distribution on the hydroxide and hydronium diffusion mechanisms is different. For AEMs, we found that hydroxide ion diffusion is mostly vehicular. This type of diffusion occurs when the water distribution is non-uniform but gives rise to both rst and second solvation shells for the hydroxide ions. 68 However, in this study, we nd that for PEMs, hydronium ion diffusion is structural rather than vehicular, with the participation of the anions, according to the reaction SO 3 À + H 3 O + 4 SO 3 H + H 2 O (as discussed in detail in the previous Section). Furthermore, we nd that the differences between AEMs and PEMs lie in the essence of the membrane materials. The region between each pair of cations in AEMs was found to create a bottleneck for hydroxide diffusion such that only specic solvation complexes are mobile, leading to a suppression of hydroxide ion diffusion. 68,70 In contrast to AEMs, the model studied here indicates that the anions in PEMs, rather than creating such a bottleneck for hydronium diffusion, become active participants in the hydronium diffusion mechanism via the reaction SO 3 À + H 3 O + 4 SO 3 H + H 2 O previously mentioned (see Fig. 7), suggesting that under the right hydration conditions, the presence of the anions in the PEM model would promote, rather than suppress, hydronium diffusion. We believe that elucidating the differences between the diffusion mechanisms of the hydroxide and hydronium ions in AEMs and PEMs is the rst step towards the discovery and determination of key design principles of new, stable cation or anion conductive membrane materials with high ion conductivity for use in emerging fuel cell technologies and other electrochemical device applications.

In conclusion, in this study, we aimed to uncover atomistic details of the hydronium ion diffusion mechanism in PEMs under conned environments with a layered morphology 7 in order to elucidate the inuence of the hydration value on the hydronium ion transport process and to compare our ndings to previously studied hydroxide diffusion mechanisms in analogous AEM environments. For this purpose, we simulated two different idealized PEM environments under two hydration conditions (l ¼ 3 and 4). We found that the water distribution is uniform only for system l4. Reducing the number of water molecules per cation by one (system l3) results in a water distribution that is non-uniform, which is associated with a dearth of water molecules and results in void areas within the simulation cell. We nd that the non-uniform water distribution results in an incomplete second solvation shell for the hydronium ion, fewer water molecules in the vicinity of a sulfonate oxygens (i.e., SO 3 À ), and a higher probability of obtaining a CN value of $1 for the oxygen located next to SO 3 À . The existence of these conditions increases the probability that the hydronium ion will react with the anion according to the reaction SO 3

À + H 3 O + 4 SO 3 H + H 2 O, which was found to be an essential part of the hydronium ion diffusion mechanism in system l3. Furthermore, we nd that under optimal hydration conditions (l3) the anions in the model PEMs promote hydronium conductivity by playing an active role in the hydronium diffusion mechanism. The results presented in this study enable us to suggest idealized hydration conditions and diffusion mechanisms for achieving high hydronium ion conductivity in high-performance PEM fuel cell devices. We believe this work is the rst to provide atomistic insight and a preliminary fundamental understanding of the unique hydronium ion diffusion mechanism in idealized PEMs by using AIMD simulations.