Hydration-Controlled Coupling of Water Polarization and Proton Transport in Nafion Under External Electric Fields

1. Introduction

Proton exchange membrane fuel cells (PEMFCs) are among the most promising energy conversion devices for sustainable power generation, relying on the efficiency and durability of proton-conducting membranes such as Nafion. Nafion, a perfluorosulfonic acid (PFSA) polymer, combines a hydrophobic polytetrafluoroethylene (PTFE) backbone with perfluoroether side chains terminated by sulfonic acid groups, forming a complex nanophase-separated morphology that governs water uptake and ion transport [1,2,3,4]. Understanding the molecular-scale structure and dynamics of Nafion is therefore essential for the design of improved membranes with enhanced proton conductivity and mechanical stability, an area that remains under intense development [5].

Atomistic and coarse-grained molecular dynamics (MD) simulations have played a crucial role in elucidating the relationship between hydration, morphology, and transport phenomena in Nafion. Early simulation studies revealed the existence of continuous water domains within the hydrophilic regions and provided insight into the formation of interconnected channels facilitating proton transport [1,6]. Subsequent works introduced more refined models of hydrated Nafion, demonstrating how sulfonic acid aggregation and water clustering drive nanophase segregation [2,3]. Jinnouchi et al. [7] further established that the mobility of hydronium ions is strongly coupled to the dynamic reorganization of the polymer–water interface, while later simulations clarified oxygen permeation and gas solubility in the membrane under operational conditions [8]. Atomistic simulation remains an active tool in this area: recent studies have resolved proton transport across the interfaces of the nanophase-separated domains in short-side-chain ionomers such as Aquivion [9] and characterized transport in Nafion-based composite membranes [10].

More recent studies have focused on the morphological response of Nafion thin films at electrode interfaces. Kang et al. [11] and Cheng et al. [12] employed high-resolution simulations and spectroscopic analyses to demonstrate that interfacial confinement, water content, and substrate hydrophilicity significantly alter the orientation of side chains and the connectivity of water domains, thus influencing proton transport near catalytic surfaces. These results reveal the multiscale complexity of Nafion morphology and point to the need for accurate atomic-level models that capture both electrostatic interactions and dielectric polarization effects. More recently, molecular dynamics has also been used to probe how external electric fields, together with temperature and hydration, affect proton transport in ordered and amorphous proton exchange membranes [13]; the field response of these materials is attracting growing attention.

In this work, we extend previous atomistic modeling efforts by performing density functional theory (DFT) calculations to obtain Hirshfeld-derived atomic charges for Nafion monomer units, which were subsequently used in classical molecular dynamics simulations with GROMACS. Rather than treating the dielectric, structural, and transport properties of hydrated Nafion as separate observables, we focus on how they are coupled under an external electric field. We show that the equilibrium dielectric response of the aqueous phase, the build-up and saturation of the field-induced water polarization, the reorganization of the hydronium coordination shell, and the onset of directed hydronium transport all trace back to one variable: the degree of electrostatic screening set by the water content. More water means stronger screening of the attraction between the mobile ions and the sulfonate network, and the characteristic field of each of these responses moves upward accordingly. By scanning hydration levels ($\lambda$) and external fields, we thus map the coupling between water polarization and proton transport that governs the field response of perfluorinated ionomer membranes.

2. Molecular Models and Method

The Nafion monomer employed in the present simulations consists of a perfluorinated PTFE-like backbone composed of CF_x units, bearing four perfluoroether side chains terminated by sulfonic acid groups ($-{\mathrm{SO}}_3\mathrm{H}$), as shown in Figure 1. The detailed chemical structure is presented in Figure 2. In this molecular model, four sulfonated side chains are covalently attached to the fluorocarbon backbone at regular intervals. The resulting repeat unit has a molecular weight of 3215 g mol⁻¹ and an equivalent weight of EW = 804 g eq⁻¹. This value is lower than that of standard Nafion 117 (EW ≈ 1100 g eq⁻¹). Importantly, our model retains the standard Nafion long side chain (Figure 2); the reduced EW arises from a higher density of these side chains along the backbone, not from side-chain shortening. It thus represents the high-sulfonation-density end of the long-side-chain Nafion family, and although this EW coincides with the range of short-side-chain ionomers such as Aquivion/Hyflon (EW ≈ 700–870 g eq⁻¹) [14], the two differ in side-chain chemistry. The higher density of sulfonic groups increases the concentration of exchange sites per unit mass, which is expected to promote water clustering and the percolation of ionic domains at a given hydration level $\lambda$. We do not vary the equivalent weight in this work; accordingly, the absolute field scales reported below may shift for membranes of higher EW, while the coupling mechanism we describe is expected to be general.

The bulk Nafion system was modeled by placing 36 Nafion monomers together with 10 hydronium (H₃O⁺) and 10 hydroxide (OH⁻) ions, together with a number of water molecules (H₂O) set by the hydration level $\lambda$, defined as the number of water molecules per sulfonic acid group. With four sulfonic groups per monomer and 36 monomers, the system contains 144 sulfonic acid groups in total. In the present topology these groups are modeled in their protonated −SO₃H form, so that each monomer is electrically neutral, and the corresponding number of water molecules is $144\lambda$, so that higher $\lambda$ corresponds to a larger water content.

The ionic composition itself deserves a remark. The mobile carriers consist of a small, fixed population of 10 H₃O⁺ and 10 OH⁻ ions, so that the system is charge-neutral by construction. This differs from the fully dissociated limit expected for a perfluorosulfonic superacid at these hydration levels, in which each sulfonic group would donate its proton to the aqueous phase, yielding of order 144 H₃O⁺ and 144 SO₃⁻ per box. What our setup probes, therefore, is the field response of a dilute carrier population within a predominantly protonated matrix, not the transport of a fully dissociated membrane. The diffusion coefficients and drift thresholds reported below should be read in that light: they characterize the field-driven mobility of the explicit carriers and the reorganization of the H₃O⁺–SO₃(H) coordination, and provide a lower bound on the proton transport of the fully hydrated, fully dissociated membrane.

All systems were contained in an orthorhombic simulation box of fixed dimensions $4.80\times 4.50\times 6.80$ nm³, common to every hydration level. Periodic boundary conditions were applied in all three spatial directions to mimic bulk behavior, as illustrated in Figure 3.

The bonded and Lennard–Jones parameters of the Nafion monomer were taken from the GROMOS 54A7 united-atom force field [15], the initial topology being generated with the Automated Topology Builder (ATB) [16,17]. Non-bonded interactions follow the GROMOS geometric combination rules. The atomic charges provided by the ATB were then replaced by charges derived from a Hirshfeld population analysis, following the strategy described in Refs. [18,19,20]. The electronic structure calculations were performed using Density Functional Theory (DFT) with the ORCA 5.0.4 package, employing the B3LYP functional and the def2-SVP basis set. Geometry optimizations were followed by single-point energy calculations, and the Hirshfeld method was applied to extract the partial atomic charges. These charges were uniformly scaled by a factor of 1.2 to reproduce the reference molecular dipole moment [18], and the scaled values were assigned to the Nafion topology.

The molecular dynamics (MD) simulations were conducted with the GROMACS 2024.5 software. Systems consisting of hydrated Nafion were constructed using different hydration levels, characterized by the hydration parameter $\lambda$ (number of water molecules per sulfonic acid group). Values of $\lambda$= 3, 9 and 15 were studied to represent different levels of membrane hydration. Water was represented by the rigid four-site TIP4P/Ice model [21], while the hydronium (H₃O⁺) and hydroxide (OH⁻) ions were represented by rigid, non-polarizable models carrying formal charges of $+1$ and $-1$, respectively, with Lennard–Jones interactions assigned only to the oxygen atom. Following the approach of Elliott et al. [6], who built their hydronium fragment directly from the water molecule to preserve internal consistency with the solvent, the ionic oxygen Lennard–Jones parameters were taken from the SPC/E water model [22], and the ion geometry was inherited from a standard water model. This entails a slight mismatch between the van der Waals description of the ions and that of the TIP4P/Ice solvent; because the ions form a dilute species and the dielectric analysis is restricted to the water subsystem, we expect this mismatch to have a limited effect on the reported trends, but state it explicitly as a limitation. We note that TIP4P/Ice was parametrized to reproduce the phase diagram of ice rather than the dielectric or dynamical properties of ambient liquid water; owing to its four-site charge distribution it yields, at 300 K, a lower static dielectric constant and a slower self-diffusion than real water. The dielectric and transport quantities reported below should therefore be read against the intrinsic bulk behavior of this model, not against experimental water directly. After energy minimization and a 2.5 ns equilibration ($5\times {10}^6$ steps), production runs were carried out in the NVT ensemble at 300 K, with the temperature maintained by a Nosé–Hoover thermostat (${\tau }_T=0.05$ ps). The equations of motion were integrated with a leap-frog algorithm and a time step of 0.5 fs, and each production trajectory spanned 10 ns ($2\times {10}^7$ steps). A Verlet cutoff scheme was used with a common cutoff of 1.0 nm for the Lennard–Jones and real-space electrostatic interactions; long-range electrostatics were treated with the Particle Mesh Ewald (PME) method and a long-range dispersion correction was applied to both energy and pressure. All bonds involving hydrogen were constrained with the LINCS algorithm, and the water and ion geometries were kept rigid. The main simulation parameters are summarized in Table 1.

An external electric field was applied along the z direction. In the presence of the field, the total potential energy of the system is expressed as

$$U_{\mathrm{total}}=\sum U_{ab}-E\sum _i{q}_i{z}_i,$$

(1)

where E is the magnitude of the external electric field, $q_i$ is the partial charge of particle i, and $z_i$ is its Cartesian coordinate along the z direction. The term $U_{ab}$ represents the intermolecular interaction energy including electrostatic and Lennard–Jones contributions. This external-field scheme is the same one we have employed in previous molecular dynamics studies of electric-field effects on clathrate hydrates [23,24].

The electric field strength was varied between 0 and 4.0 V nm⁻¹, with the specific values chosen for each hydration level listed in Table 2. The field magnitudes explored here (up to 4.0 V nm⁻¹) are several orders of magnitude larger than the average macroscopic field across an operating PEMFC membrane ($\sim {10}^{-4}$–${10}^{-2}$ V nm⁻¹). Such elevated fields are a standard device in non-equilibrium molecular dynamics of membranes and electrolytes [13,23,24]: they are required to drive a measurable ionic response within the accessible simulation time, and they probe the local electrostatic scale at which the sulfonate–ion coupling operates. The drift thresholds reported below should therefore be read as a probe of the intrinsic coupling between polarization and transport, not as literal operating conditions; the physically meaningful result is the systematic shift of that coupling with hydration, not the absolute field values. All simulations were carried out in the canonical (NVT) ensemble at constant temperature of 300 K and fixed volume. A single simulation cell of fixed dimensions was used for all hydration levels and field strengths, and no constant-pressure (NPT) equilibration was performed. This choice is deliberate. Under the strong static external fields applied here, the field continuously does work on the net-mobile ionic species and drives the system into a field-driven steady state; barostat-based pressure control, which presupposes mechanical equilibrium and an unbiased virial, is therefore not well defined under the field, and a fixed-volume ensemble is the controlled and reproducible setting for a systematic scan over $\lambda$ and E. Using an identical cell for every system further removes the volume as a confounding variable: differences in the computed observables reflect hydration and field, not concurrent changes in cell size, and the zero-field dielectric reference and the field response are evaluated in the same cell. A consequence of this protocol is that, since water is added to a cell of fixed size, the mass density increases with hydration, from $1.40$ to $1.58$ and $1.75$ g cm⁻³ for $\lambda =3$, 9 and 15. We therefore treat the systems as a controlled-volume series rather than as reproductions of the experimental swelling isotherm; because absolute transport coefficients depend on density, the cross-$\lambda$ comparison below should be read as an analysis of the coupling and of the trends at fixed volume, not as a prediction of absolute mobilities at each hydration’s equilibrium density. For each $\lambda$ value, the total dipole moment, its longitudinal component, and the Kirkwood factor were extracted from the trajectories. The static dielectric constant was obtained from the Kirkwood–Fröhlich dipole fluctuation formula implemented in GROMACS, $ϵ-1\propto (\langle {\mathbf{M}}^2\rangle -{\langle \mathbf{M}\rangle }^2)/\left(V{k}_{B}T\right)$, where $\mathbf{M}$ is the total dipole of the aqueous subsystem. We restrict the fluctuation formula to the water molecules because the mobile ionic species (H₃O⁺ and OH⁻) carry a net charge: their displacement contributes to the ionic conductivity rather than to the orientational polarization, and including them would conflate the two responses and make $\mathbf{M}$ ill-defined under periodic boundary conditions. This expression is a linear-response, equilibrium relation: it is physically meaningful only in the absence of an external field, where the spontaneous dipole fluctuations sampled by the simulation define the zero-field susceptibility. Under a finite field the dipoles align and saturate, the fluctuations are suppressed, and the formula no longer yields a dielectric constant. We therefore report $ϵ$ only at $E=0$; the field-dependent polarization response is characterized instead by $\langle M_z\rangle$ and the alignment fraction $|M_z|/|\mathbf{M}|$, with the finite-volume Kirkwood factor $G_k$ used only as an apparent indicator of dipolar correlation. These analyses enabled the evaluation of the structural, polarization, and transport response of Nafion under different hydration states and external field conditions.

Statistical uncertainties on the zero-field dielectric constant were estimated by block averaging: the 8 ns of production sampled for the dipole (the first 2 ns being discarded as equilibration) was divided into four 2 ns blocks, and the reported uncertainty is the standard deviation across blocks. Our full-trajectory $ϵ$ reproduces the built-in GROMACS estimate to within numerical precision at all three hydration levels, confirming the implementation. The block-to-block scatter is smallest at low hydration ($\lambda =3$: $ϵ=2.9\pm 0.3$) and largest at high hydration ($\lambda =15$: $ϵ=17.2\pm 2.9$). This trend follows the collective nature of the dielectric response: the autocorrelation time of the aqueous dipole grows from $\approx 0.19$ ns at $\lambda =3,9$ to $\approx 0.37$ ns at $\lambda =15$, so that each fixed-length block contains correspondingly fewer independent samples of the slowly relaxing total dipole $\mathbf{M}$. The high-$ϵ$, strongly correlated regime is therefore intrinsically the hardest to converge, and the same collective correlations that enhance $ϵ$ also govern its slower convergence.

3. Results and Discussion

The high sulfonate density has two direct consequences for ion transport. First, even in the absence of an applied field, the H₃O⁺ ions at low hydration are already embedded in a dense network of neighboring sulfonic groups, which templates the vehicular motion of the intact H₃O⁺ ions through the sulfonic-group array. This is reflected in the zero-field diffusion coefficients, which decrease monotonically with hydration: $D(E=0)\approx 0.023$, $0.006$, and $0.00008\times {10}^{-5}$ cm² s⁻¹ for $\lambda =3$, 9, and 15, respectively. The progressive dilution of the sulfonate network by water thus reduces the intrinsic, field-free proton mobility as hydration increases. We emphasize that this trend is a feature of the controlled-volume protocol (Section 2): because the cell size is held fixed, increasing $\lambda$ raises the mass density (from 1.40 to 1.75 g cm⁻³), so the added water crowds the fixed cell and reduces the free volume available for ion motion. The decreasing zero-field mobility therefore reflects this increasing confinement; it does not contradict the enhanced conductivity observed in constant-pressure, freely swelling membranes.

Second, the water content controls the minimum field required to drive H₃O⁺ out of the diffusive regime into field-driven drift. As additional water molecules are incorporated, they progressively screen the electrostatic interaction between the sulfonate groups and the mobile ions, so that a stronger external field is needed to overcome the structural constraints of the hydrated matrix. This screening is quantified by the systematic shift of the drift threshold: $E\approx 0.7$, $1.4$, and $3.7$ V nm⁻¹ for $\lambda =3$, 9, and 15, respectively (see Section 3.2). Below each threshold, however, the field does enhance diffusive transport. The relative enhancement (the factor by which D increases from zero field to the diffusive maximum) grows with hydration, reflecting the larger structural reorganization that the field induces in more water-rich environments.

These observations are consistent with the classical cluster–channel model proposed by Hsu and Gierke [25] for Nafion membranes, in which sulfonic acid groups aggregate into hydrophilic ionic clusters interconnected by narrow channels. In the present model the high sulfonate density (EW = 804 g eq⁻¹) promotes well-defined ionic domains even at moderate hydration levels. The applied field enhances the longitudinal connectivity of these domains by partially aligning the hydrated ionic regions, but the field strength required to achieve this alignment is set by the degree of electrostatic screening, which increases monotonically with $\lambda$. The quantitative thresholds reported here thus provide a direct measure of how membrane hydration controls the field-driven onset of directed hydronium transport.

3.1. Zero-Field Dielectric Response and Field-Induced Dipolar Alignment

At zero field, where the fluctuation formula defines a genuine equilibrium dielectric constant, the aqueous-subsystem $ϵ$ grows monotonically with hydration, from $ϵ\approx 2.9\pm 0.3$ at $\lambda =3$ to $ϵ\approx 6.0\pm 0.3$ at $\lambda =9$ and $ϵ\approx 17.2\pm 2.9$ at $\lambda =15$ (Figure 4). This growth is driven by collective dipolar correlations: the static Kirkwood factor $G_k$ exceeds unity and increases with hydration, whereas the water molecular dipole is fixed by the (rigid) water model. The dielectric enhancement with hydration therefore reflects the increasing mutual alignment of the water dipoles as the aqueous domains percolate, not a change in the single-molecule dipole.

These values should be read relative to the bulk dielectric constant of the water model itself. Owing to its four-site charge distribution, TIP4P/Ice yields a static dielectric constant below that of experimental water at 300 K, so the growth from $ϵ\approx 2.9$ at $\lambda =3$ to $\approx 17.2$ at $\lambda =15$ reflects the progressive build-up of the aqueous dielectric response toward the model’s own bulk value as the hydrophilic domains percolate; it is not an approach to the experimental bulk permittivity.

Under a finite field the dipole-fluctuation formula no longer defines a dielectric constant: the field aligns the dipoles and suppresses the very fluctuations on which the linear-response relation rests, so the quantity it returns loses its physical meaning as $ϵ$. We therefore characterize the field response directly through the polarization of the aqueous subsystem. The longitudinal water dipole $\langle M_z\rangle$ rises steeply at low fields and then saturates as the dipoles become fully aligned (Figure 5); for $\lambda =15$, for instance, $\langle M_z\rangle$ grows from essentially zero at $E=0$ to $\approx 3000$ D already at $E=1$ V nm⁻¹, and then approaches its saturation value of $\approx 4500$ D only gradually, with increments of a few percent per field step beyond $E\approx 3$ V nm⁻¹. Both the saturated dipole and the field needed to approach it grow with hydration: the $\lambda =3$ and $\lambda =9$ curves approach saturation by $E\approx 1.5$ and 2 V nm⁻¹, respectively, whereas for $\lambda =15$ the approach extends to $E\approx 3$–4 V nm⁻¹. The response is strongly anisotropic, with the longitudinal component dominating over the transverse ones for all $\lambda$ and E, and the two observables separate the orientation of the collective dipole from its magnitude. The alignment fraction $|M_z|/|\mathbf{M}|$ jumps from its near-isotropic zero-field value to $\gtrsim 0.99$ already at $E=0.5$ V nm⁻¹ for $\lambda =9$ and 15, and above $\approx 0.98$ by $E\approx 0.7$ V nm⁻¹ for $\lambda =3$ (Figure 6): the direction of the collective water dipole locks onto the field essentially immediately, at fields well below those at which its magnitude saturates. The continued growth of $\langle M_z\rangle$ at higher fields therefore reflects the progressive alignment of the individual molecular dipoles within an already field-oriented collective moment, and it is this magnitude saturation, not the orientational locking, that carries the hydration-dependent field scale, consistent with field-induced ordering amplified by the percolation of the aqueous domains.

This field-driven alignment does not occur in isolation: it is coupled to the structural and transport reorganization characterized in the following sections, all sharing a common field scale, which is again set by hydration through electrostatic screening: the more water is present, the stronger the field required to align the ionic network and to drive the carriers out of the diffusive regime. Concretely, for $\lambda =3$ the diffusion coefficient D reaches its diffusive-regime maximum at $E\approx 0.6$ V nm⁻¹, the drift onset occurs at $E\approx 0.7$ V nm⁻¹, and the S–H₃O⁺ first-peak height $g_1$ peaks at $E\approx 1.0$ V nm⁻¹, all within the same narrow window; for $\lambda =9$ the diffusive maximum of D ($E\approx 1.0$ V nm⁻¹) and the field-induced recovery maximum of $g_1$ ($E\approx 1.5$ V nm⁻¹) likewise cluster around the drift onset at $E\approx 1.4$–$1.6$ V nm⁻¹, and for $\lambda =15$ the same sequence is displaced to $E\approx 3.7$ V nm⁻¹. That the maxima of cooperative diffusion and of first-shell H₃O⁺–SO₃(H) occupancy track each other as $\lambda$ varies points to a common origin for the polarization, structural, and transport responses: a screening-controlled build-up of correlated ionic ordering.

3.2. Hydronium Transport Coupled to the Polarization Response

The polarization response identified above is directly coupled to ion transport. To make this coupling explicit, we extracted the diffusion coefficients D of hydronium ions from the linear regime of the mean squared displacement (MSD) curves, fitting the time interval between 2000 and 6000 ps. This window excludes the initial ballistic regime and minimizes the nonlinear effects that emerge at long times under strong external fields. Because a uniform field exerts a net force on the charged carriers, the long-time motion ceases to be purely diffusive once the field is strong enough; we therefore restrict the quantitative analysis to the diffusive regime and treat the large-field data separately as a drift threshold (see below). The resulting $D\left(E\right)$ curves for $\lambda =3$, 9 and 15 are shown in Figure 7, and the corresponding values are listed in Table 2.

Within the diffusive regime, $D\left(E\right)$ displays a nonmonotonic, hydration-dependent behavior that parallels the build-up of the water polarization described above. For $\lambda =3$, D first decreases weakly as the field rises from zero ($0.023\to 0.018\times {10}^{-5}{\mathrm{cm}}^2{\mathrm{s}}^{-1}$ between $E=0$ and $0.3\mathrm{V}{\mathrm{nm}}^{-1}$), then rises sharply to a local maximum near $E\approx 0.6\mathrm{V}{\mathrm{nm}}^{-1}$ ($1.74\times {10}^{-5}{\mathrm{cm}}^2{\mathrm{s}}^{-1}$). For $\lambda =9$ the same feature appears but shifted to higher field: D peaks at $E\approx 1.0\mathrm{V}{\mathrm{nm}}^{-1}$ ($0.66\times {10}^{-5}{\mathrm{cm}}^2{\mathrm{s}}^{-1}$) and is then suppressed at $E=1.2$–$1.3\mathrm{V}{\mathrm{nm}}^{-1}$ ($0.33$ and $0.27\times {10}^{-5}{\mathrm{cm}}^2{\mathrm{s}}^{-1}$). This local-maximum-then-suppression structure reflects the same cooperative ordering that underlies the growth of the water polarization, and progressively weakens with hydration: for $\lambda =15$ it survives only as a shallow plateau near $E=1.5$–$2.0\mathrm{V}{\mathrm{nm}}^{-1}$ ($0.023$ and $0.022\times {10}^{-5}{\mathrm{cm}}^2{\mathrm{s}}^{-1}$), where the variation falls within the scatter of the fits and should be regarded as indicative rather than conclusive.

The most robust trend is the systematic shift of the field scale with hydration. The onset of field-driven drift, identified by the abrupt rise of the apparent D by several orders of magnitude, moves to progressively higher fields as $\lambda$ increases: $E\approx 0.7\mathrm{V}{\mathrm{nm}}^{-1}$ for $\lambda =3$, $E\approx 1.4$–$1.6\mathrm{V}{\mathrm{nm}}^{-1}$ for $\lambda =9$, and $E\approx 3.7$–$4.0\mathrm{V}{\mathrm{nm}}^{-1}$ for $\lambda =15$. This monotonic displacement is a direct manifestation of the electrostatic screening discussed above: the additional water progressively screens the interaction between the sulfonate groups and the mobile ions, so that a stronger field is required to overcome the structural constraints of the hydrated matrix and drive the carriers into the drift regime. The same screening that sets the field scale of the water polarization thus also displaces the transport threshold: the two responses move together because they share their origin.

Beyond these onset fields, the apparent diffusion coefficients increase by several orders of magnitude and no longer describe equilibrium diffusion: the motion is dominated by field-driven drift rather than thermal diffusion. These large-field values are reported in Table 2 for completeness, marked with a dagger, but are excluded from the quantitative analysis above because they do not represent equilibrium transport properties.

Finally, a limitation intrinsic to the present approach should be stated explicitly. Because the simulations use a classical, non-reactive force field, proton transport is captured only through the vehicular motion of intact H₃O⁺ ions; the structural (Grotthuss) mechanism, in which the excess proton is relayed through the hydrogen-bond network, is not represented. As already emphasized by Elliott et al. [6], neglecting charge transfer in classical MD means that the computed diffusion coefficients should be regarded as a lower bound on the true proton mobility, particularly at high hydration where structural diffusion contributes an increasing fraction of the experimental conductivity. A reactive treatment, such as the multistate empirical valence bond (MS-EVB) method [26] or ReaxFF [27], would be required to recover the Grotthuss contribution.

3.3. Structural Signature in the First Coordination Shell

To connect the macroscopic transport behavior described above with microscopic structural changes, we computed radial distribution functions (RDFs) $g\left(r\right)$ for three atom pairs: sulfur–sulfur (S–S), sulfur–hydronium (S–H₃O⁺), and sulfur–water (S–W), for $\lambda =3$, 9, and 15 at the same electric fields studied in the MSD analysis. The resulting distributions are shown in Figure 8.

The S–S peak at $r\approx 0.41$ nm, which corresponds to the nearest-neighbor sulfonate distance, and the S–W peak at $r\approx 0.38$ nm remain at the same position across all fields and hydration levels, and their heights show no systematic field dependence below the drift onset: the S–S height fluctuates by ∼15–20% without a clear trend, and the S–W height is essentially constant, decreasing modestly (≈10–20%) only at the largest, drift-regime fields (Figure 8 and Figure 9). Neither the polymer backbone nor the primary hydration shell of the sulfonate groups therefore undergoes a significant field-induced structural reorganization below the drift onset.

In sharp contrast, the S–H₃O⁺ distribution exhibits a clear and systematic response. Figure 9 shows the height $g_1$ of the first peak of each distribution as a function of E for the three hydration levels. For all $\lambda$, $g_1$ follows a nonmonotonic pattern that directly mirrors the ion-transport behavior: it rises toward a maximum located near the drift onset and falls once the drift regime is fully developed. For $\lambda =3$, $g_1$ rises from 0.75 at $E=0$ to 1.08 at $E=1.0$ V nm⁻¹ and drops to 0.63 at $E=1.5$ V nm⁻¹ (post-threshold). For $\lambda =9$, $g_1$ dips from 1.32 at $E=0$ to 0.81 at $E=0.5$ V nm⁻¹, recovers to $\approx 1.1$–$1.2$ around the drift onset (peaking at 1.21 at $E=1.5$ V nm⁻¹), and then falls to 0.73 at $E=2.0$ V nm⁻¹, beyond the onset at $E\approx 1.4$–$1.6$ V nm⁻¹ where D has already jumped by two orders of magnitude. For $\lambda =15$, $g_1$ reaches its maximum value of 1.31 at $E=3.0$ V nm⁻¹ and drops to 1.01 at $E=4.0$ V nm⁻¹, coinciding with the abrupt rise of D from $0.79$ to $250\times {10}^{-5}$ cm² s⁻¹.

This behavior identifies a two-stage structural sequence. Below the threshold, the field progressively draws H₃O⁺ ions into closer coordination with the SO₃(H) groups, increasing the occupancy of the first coordination shell and building a field-aligned ionic network that enhances cooperative susceptibility. Once the threshold is crossed, the ions are no longer trapped at a single sulfonate site but transition to a mobile regime in which they migrate continuously between sites; the resulting reduction in residence time at any fixed site is reflected in the drop of $g_1$. The drift threshold thus marks the crossover from localized coordination to delocalized transport, and the structural signature in $g\left(r\right)$ confirms that this crossover is the common microscopic origin that couples the transport threshold in $D\left(E\right)$ to the saturation of the field-induced dipolar alignment. The systematic shift of this threshold to higher fields with increasing $\lambda$ is again a direct consequence of the electrostatic screening discussed in the preceding sections.

To complement the peak-height analysis with an integral measure, we computed running coordination numbers by integrating each RDF to a fixed cutoff $r_c$ taken, for every hydration level, from the position of the zero-field first minimum, so that the field dependence reflects changes in shell occupancy rather than in the integration boundary ($r_c\approx 1.0$ nm for $\lambda =9$ and $0.79$ nm for $\lambda =15$). Uncertainties were obtained by block averaging over four 2 ns blocks, as for the dielectric analysis. Because the S–H₃O⁺ RDF is referenced to the hydronium ions, this coordination number counts the sulfonic groups surrounding each hydronium. At $\lambda =3$ the S–H₃O⁺ correlation remains weak ($g_1$ stays at $\lesssim 1.1$ and the distribution lacks a clear first minimum), so no well-defined first coordination shell exists at low hydration and the field only begins to build incipient coordination; the coordination analysis is therefore meaningful only at $\lambda =9$ and 15, where the shell is well defined. There, the first-shell occupancy confirms the two-stage picture of $g_1$: it is high below the drift threshold and collapses once the drift regime is entered. For $\lambda =15$ the number of coordinating sulfonic groups rises from $1.16\pm 0.03$ at $E=0$ to $1.65\pm 0.13$ near the drift-onset field ($E=3.5$ V nm⁻¹) and falls to $1.10\pm 0.16$ at $E=4.0$ V nm⁻¹; for $\lambda =9$ it remains at $\approx 2.8$–$3.2$ up to the drift onset ($3.24\pm 0.15$ at $E=0$; $2.80\pm 0.22$ and $3.03\pm 0.46$ at $E=1.4$ and $1.5$ V nm⁻¹) and drops to $1.89\pm 0.37$ at $E=2.0$ V nm⁻¹, beyond the onset, where D has already jumped by two orders of magnitude. The nonmonotonic $g_1$ signature therefore reflects a genuine change in first-shell occupancy, not a normalization artifact, and the drift onset coincides with the loss of hydronium–sulfonate coordination. The S–S and S–W coordination numbers show no counterpart of the S–H₃O⁺ collapse, varying by $\lesssim 16\%$ over the full field range: the S–S number decreases monotonically with hydration ($2.72\pm 0.09$, $1.87\pm 0.07$, $1.77\pm 0.04$ for $\lambda =3$, 9, 15 at $E=0$), reflecting dilution of the sulfonate network, whereas the number of water molecules in the first sulfur shell is largest at high hydration ($\approx 7.0$, $6.6$, and $11.3$ for $\lambda =3$, 9, 15; the $\lambda =3$ value is inflated by its wider first shell, $r_c\approx 0.87$ nm versus $0.66$ nm at $\lambda =9$), reflecting the growing hydration of the sulfur environment, consistent with Figure 9.

4. Conclusions

At zero field, the static dielectric constant of hydrated Nafion is set by collective dipolar correlations: the aqueous-subsystem $ϵ$ grows with hydration (from $\approx 2.9$ at $\lambda =3$ to $\approx 17.2$ at $\lambda =15$) while the single-molecule dipole remains essentially constant, so the enhancement is cooperative in origin and tracks the percolation of the aqueous domains. Under an external field the dipole-fluctuation formula no longer defines a dielectric constant; the meaningful response is the polarization itself, which we follow through the longitudinal dipole $\langle M_z\rangle$ and the alignment fraction $|M_z|/|\mathbf{M}|$. The orientation of the collective water dipole locks onto the field already at low fields for every hydration level; the field then progressively saturates the magnitude of the water polarization, reorganizes the first H₃O⁺–SO₃(H) coordination shell, and drives a crossover from diffusive to drift-dominated proton transport. The central result is the coupling among these responses. The saturation field of the water polarization, the first-shell occupancy, and the transport threshold all move to higher fields as hydration increases, and they move together. We attribute this concerted shift to electrostatic screening: the added water weakens the sulfonate–ion interaction, and in doing so it sets at once the equilibrium dielectric response and the field required for directed hydronium transport. Because all systems share a fixed simulation cell, the mass density increases with hydration, and the absolute diffusion coefficients (in particular their decrease with $\lambda$) reflect this controlled-volume protocol rather than the behavior of a freely swelling membrane. Our conclusions concern the coupling between the polarization, structural, and transport responses and the concerted shift of their characteristic field with hydration, which are robust to this choice; they are not intended as predictions of absolute mobilities or of the hydration dependence of conductivity at each equilibrium density. Viewed this way, the dielectric, structural, and transport signatures form a single picture of the field response of perfluorosulfonated membranes, one in which hydration acts as a single knob that tunes both the polarization and the onset of proton conduction under bias.