Simulation Study of Ion Diffusion in Charged Nanopores with Anchored Terminal Groups

TL;DR

This study uses coarse-grained simulations to explore how ion diffusion in charged nanopores depends on pore size, sidechain length, and water structuring, revealing conditions for enhanced ion mobility.

Contribution

It demonstrates that ion diffusion is maximized with medium-length sidechains and specific pore parameters, providing new insights into nanopore design for improved ion transport.

Findings

Ion diffusion depends strongly on pore size and sidechain length.

Maximum ion diffusion occurs with medium-sized sidechains and specific pore geometries.

Radial charge separation suppresses ion diffusion in certain configurations.

Abstract

We present coarse-grained simulation results for enhanced ion diffusion in a charged nanopore grafted with ionomer sidechains. The pore surface is hydrophobic and its diameter is varied from 2.0 nm to 3.7 nm. The sidechains have from 2 to 16 monomers (united atom units) and contain sulfonate terminal groups. Our simulation results indicate a strong dependence of the ion diffusion along the pore axis on the pore parameters. In the case of short sidechains and large pores the ions mostly occupy the pore wall area, where their distribution is strongly disturbed by their host sulfonates. In the case of short sidechains and narrow pores, the mobility of ions is strongly affected by the structuring and polarization effects of the water molecules. In the case of long sidechains, and when the sidechain sulfonates reach the pore center, a radial charge separation occurs in the pore. Such charge separation suppresses the ion diffusion along the pore axis. An enhanced ion diffusion was found in the pores grafted with medium-size sidechains provided that the ions do not enter the central pore area, and the water is less structured around the ions and sulfonates. In this case, the 3D density of the ions has a hollow-cylinder type shape with a smooth and uninterrupted surface. We found that the maximal ion diffusion has a linear dependence on the number of sidechain monomers. It is suggested that the maximal ion diffusion along the pore axis is attained if the effective length of the sidechain extension into the pore center (measured as twice the gyration radius of the sidechain with the Flory exponent 1/4) is about 1/3 of the pore radius.

Tables

Table 1. Table 1: Simulation parameters for the case A , case B , and case C runs. All distances are given in nm units. The pore length is L 𝐿 L =16 nm and the number of grafting points on the pore surface is N s subscript 𝑁 𝑠 N_{s} =200. The parameters d ∗ superscript 𝑑 d^{*} and r s ​ s ∗ superscript subscript 𝑟 𝑠 𝑠 r_{ss}^{*} for the case A correspond to the reduced value of N s = 100 subscript 𝑁 𝑠 100 N_{s}=100 . L s subscript 𝐿 𝑠 L_{s} is the number of monomer units per sidechain. The physical length l s subscript 𝑙 𝑠 l_{s} of a stretched sidechain with zero dihedral angle is l s subscript 𝑙 𝑠 l_{s} = L s ​ b ​ sin ⁡ ( θ / 2 ) subscript 𝐿 𝑠 𝑏 𝜃 2 L_{s}b\sin(\theta/2) . Overlined quantities on the left column represent the fixed pore parameter for the corresponding runs.

case A runs	A1	A2	A3	A4	A5	A6	A7
${\overline{L}}_s$	4	4	4	4	4	4	4
$d$	2.03	2.33	2.70	2.89	3.12	3.37	3.55
$\lambda$	1	3	6	8	10	13	15
$r_{ss}$	0.65	0.71	0.77	0.80	0.84	0.87	0.90
$d^{\ast }$	1.54	1.75	2.01	2.17	2.30	2.49	2.61
$r_{ss}^{\ast }$	0.78	0.84	0.92	0.96	0.99	1.04	1.07
case B runs	B1	B2	B3	B4	B5	B6	B7
$L_s$	2	4	6	8	10	12	16
$\overline{d}$	2.89	2.89	2.89	2.89	2.89	2.89	2.89
$\lambda$	9.3	8	6.7	5.7	4.2	2.8	1
$r_{ss}$	0.80	0.80	0.80	0.80	0.80	0.80	0.80
case C runs	C1	C2	C3	C4	C5	C6	C7
$L_s$	4	8	10	12	14	16	18
$d$	2.89	3.15	3.30	3.40	3.50	3.60	3.70
$\overline{\lambda }$	8	8	8	8	8	8	8
$r_{ss}$	0.80	0.84	0.86	0.88	0.89	0.91	0.92

Equations

$$r_{ss}\approx \frac{1}{\sqrt{n_s}}$$

$$U(\vec{r})=\sum _i{U}_b^i+\sum _j{U}_{\theta }^j+\sum _m{U}_{\phi }^m+\sum _{k,l}{U}_{nb}\left(|{\vec{r}}_k-{\vec{r}}_l|\right)$$

$$6Dt=\frac{1}{N_s}\sum _{i=1}^{N_s}⟨|{\vec{r}}_i(t)-{\vec{r}}_i(0){|}^2⟩$$

$$2D_{x}t=\frac{1}{N_s}\sum _{i=1}^{N_s}⟨|x_i(t)-x_i(0){|}^2〉⟩$$

$$g_{jj}=\frac{1}{4\pi r^2\Delta r{N}_s}\sum _{i=1}^{N_s}\delta \left(\vec{r}-{\vec{r}}_i^{(j)}\right)\theta \left(|\vec{r}-{\vec{r}}_i^{(j)}|-\Delta r\right)$$

$$U_i(r)=-\frac{e^2}{{ϵ}_{s}r}{e}^{-\frac{r}{r_D}}$$

$$r_D=\sqrt{\frac{k_{B}T{ϵ}_0{ϵ}_s}{n_e{e}^2}}=\sqrt{\frac{k_{B}T{ϵ}_0{ϵ}_{s}d}{8n_s{e}^2}}\approx r_{ss}\sqrt{\frac{k_{B}T{ϵ}_0{ϵ}_{s}d}{8e^2}}\approx 0.008r_{ss}\sqrt{\frac{{ϵ}_{s}d}{\sigma }}$$

$$\text{mbox}CN={\rho }_H^0{\int }_{\sigma }^{R_m}{g}_{S{O}_3-H}(r)4\pi r^{2}dr$$

$${\rho }_j(r^{(j)})=\frac{1}{2\pi r\Delta rL{N}_s}\sum _{i=1}^{N_s}\delta \left(r-r_i^{(j)}\right)\theta \left(|r-r_i^{(j)}|-\Delta r\right)$$

$$F=\frac{3⟨{\cos }^2(\theta )⟩-1}{2},$$

$$D=\frac{2}{3}{D}_r+\frac{1}{3}{D}_x$$

$${\rho }_O(x)=\frac{1}{\pi d^2\Delta x}\sum _{i=1}^{N_s}\delta \left(x-x_i\right)\theta \left(|x-x_i|-\Delta x\right)$$

$$R_g=b{L}_s^{\nu }$$

$$\Gamma =\frac{2l_s}{d}=\frac{4R_g}{d}=\frac{4}{d}b{L}_s^{\nu }$$

Full text

1 Introduction

Over the past decade there has been a growing interest in ion transport phenomena in restricted geometries and porous networks. The development of many applications relies on a fundamental understanding of the diffusion properties of particles in cylindrical pores. For example, transport properties of the catalyst layer of polymer electrolyte fuel cells are largely regulated by the electrostatic interaction of the ion with the charges on the cylindrical pore in the Pt electrode [1, 2, 3]. The charge storage in conducting narrow nanopores relies on the voltage-controlled accumulation of ions in a narrow metallic nanopore [4, 5, 6]. Charged nanopores are used for the partitioning of ions and proteins [7, 8, 9], for electrolyte nanofiltration [10], and for ion-current rectification [11]. A charged-pore model also has been adopted as a convenient study model for water sorption [12] and structural and kinetic characteristics of ion conductivity in fuel cell ionomers [13, 14, 15], for the self-diffusion of ions in nanoporous media [16, 17], for the permeability of the ion channels to water and ions [18], and for the electro-osmotic flow of ionic solutions in charged channels [19]. Nanopores with anchored ionizable surface groups have been used for advanced membrane separations [20], for ion selectivity [21], and for the enhanced transport of proteins [22, 23]. They also serve as a model for the ion diffusion in cylindrical pores of the membrane-ionomer interface [24], and for the permeability of the ionomer pores to specific ions [25, 26].

Ongoing research on new polyelectrolyte membranes (PEM) for fuel cell applications also focuses on the ion diffusion properties of porous networks [27, 28, 29, 30]. Ionomer chains of the PEM are mostly composed of hydrophobic inert backbones grafted with pendant side chains that are terminated by anionic headgroups. When hydrated, the terminal group ions dissociate and the sulfonates self assemble into connected clusters, creating hydrophilic pathways inside the hydrophobic backbone matrix. For effective ion transport through the membrane, the PEM needs to be hydrated and swollen in order to form a connected network of hydrophilic pathways. The hydration water, however, makes the PEM vulnerable to icing, boiling and water evaporation at low and high temperatures and to loss of its elastic properties. In this sense, designing new ionomers with low water uptake and, at the same time, with acceptable ionic conductivity, is the goal of many ongoing activities in developing advanced fuel cell membranes.

One of the possible ways to proceed in this direction is the concept of a matrix-reinforced membrane, where a porous hydrophobic matrix film is impregnated with an ion-conducting ionomer. As a support material for the matrix, ceramic films of metal oxides [31, 32], silica xerogels [33], surface-modified silica films [34], polypropylene, polytetrafluoroethylene (PTFE), polycarbonate, polysulfone, and microglass fiber fleece (see Refs. [35, 36] and references therein) have been used. Matrix reinforcement improves the mechanical stability of the modified PEM against water flooding, and increases its resistance to the crossover of fuel molecules. However, the pore walls of the matrix restrict macromolecular motion of the ionomer, causing a lack of connectivity of the hydrophilic pathways. As a result, the ion diffusion

rate

in matrix reinforced PEMs is smaller than the diffusion in the corresponding bulk ionomers. A simple increase in the pore size will not solve this issue, because the increased connectivity of the ionomer will also increase membrane permeability to the feeding gases.

An alternative approach for the matrix reinforced ionomer might be a concept of a regular porous matrix with the pore walls grafted with pendant sidechains, as schematically illustrated in Figure 1. Here the left-side picture represents a hard-wall porous matrix with regularly ordered and parallel cylindrical pores, and the right-side picture represents a top-view snapshot of the ordered macrophase separation in a poled Nafion®-like ionomer (perfluorinated polymer from DuPont) from our previous work [37]. The ionic diffusion in such poled membranes takes place along the cylindrical sulfonate aggregates [37, 38, 39, 40]. The formation of similar parallel and inverted-micelle like cylinders in a Nafion ionomer has been reported in Ref. [41]. This finding, however, was later disputed in Refs. [42, 43] where an alternative slit geometry for the sulfonate aggregates was suggested.

In narrow pores the confinement effect of the walls suppresses the freezing point of water with respect to pure water [44, 45, 46, 47, 48, 49]. As a result, the water absorbed in the membrane tends to be in a liquid state at sub-zero temperatures. In ionomer membranes at low water contents \lambda$$\leq4.8 the confined water bonds tightly to the sulfonates and becomes strongly structured. This additionally prevents the water from freezing at low temperatures [50]. In such pores the proton transport at sub-zero temperatures probably occurs according to the Grotthuss mechanism.

In this work we consider a single pore with hydrophobic and stiff walls which represents the pores in the matrix reinforced membrane shown in Figure 1. The pore size is assumed to be in the range of 2-3.7 nm, and the pore walls are grafted with terminal groups to provide similar sulfonate aggregates as in the snapshot picture in Figure 1. The proposed single pore, in addition to being totally resistant to swelling, water flooding, and gas crossover problems, has two additional remarkable features.

First, no external electric field is needed to generate continuous lanes of sulfonates. A continuous cluster network is formed self-consistently through the sulfonate-sulfonate aggregation supplemented by the hydrophobic wall-hydrophilic sulfonate and wall-water repulsive interactions [51, 52]. Second, the water content can be regulated by the pore size, the grafting density and the number of sidechain monomers. These parameters set up a non-homogeneous distribution of the ions and water molecules in the pore [51, 53, 54].

Our objective is to analyze systematically the role of the sidechain monomer number $L_{s}$, the pore diameter $d$, and the water content $\lambda$ (the number of water molecules per terminal group) in getting enhanced ion diffusion along the pore axis

compared to the ion diffusion in a bulk Nafion-like ionomer at the same hydration levels.

We suggest that the flexibility of the sidechains and the dipolar nature of the sidechain sulfonates will greatly assist the ion transport along the pore axis, provided that the pore center is free from ion and sulfonate clusters. We report high ion diffusion coefficients for water contents \lambda$$\approx8–10 and propose a new scaling rule for the pore parameters which guarantees maximal ion diffusion along the pore axis. According to our findings, the sulfonate group protrusion length into the pore center, measured as twice the gyration radius of the sidechain with the Flory exponent 1/4, should be about 1/3 of the pore radius for securing maximal diffusion rates for the ions.

The rest of the paper is organized as follows. In section 2 we describe the main parameters of the accepted model for the pore geometry and anchored sidechains. The details of the simulation runs are given in section 3. In section 4 we consider three different sets of pore parameters for the detailed analysis of the ion diffusion in charged pores, and discuss the subsequent results. Our conclusions are presented in section 5.

2 Simulation model

A schematic picture describing our simulation model for the pore geometry is shown in Figure 2. The inner surface of the cylindrical pore of length $L$ and diameter $d$ is randomly grafted with $N_{s}$ anchoring centers. These centers serve as the attachment points to the sidechains modeled as a spring-bead polymer chain of $L_{s}$ monomers. This number, $L_{s}$, will be referred to hereafter as the sidechain protrusion length. A typical structure of a sidechain with $L_{s}$=8 monomers is shown at the bottom of Figure 2. The shaded sidechain particles are hydrophobic and electrostatically neutral with zero charge, while the hydrophilic terminal group monomers at the sidechain tip, inscribed with the letters S for the sulphur atom and O3 for the oxygen group atoms, are charged.

We assume for them partial charges $q_{S}$=+1.1$|e|$, and $q_{O_{3}}$=–2.1$|e|$, [39, 55] where $e$ is the electron charge. The assigned charges $q_{S}$ and $q_{O_{3}}$ depend on the fractional charge distribution along the sidechain. For the latter there are different force-field approaches [56, 57, 58, 59, 60, 61] which put the charge $q$=–$|e|$, necessary to compensate the hydronium charge $q_{H}$=+$|e|$, either entirely on the terminal group SO3, or smear it over the sidechain monomers. For example, the following terminal group charges have been considered in the literature: $q_{S}$=+0.932$|e|$ and $q_{O_{3}}$=–1.757$|e|$ [60], $q_{S}$=+1.7$|e|$ and $q_{O_{3}}$=–2.25$|e|$ [61], $q_{S}$=+2.08$|e|$ and $q_{O_{3}}$=–2.57$|e|$ [59], $q_{S}$=+1.284$|e|$ and $q_{O_{3}}$=–1.862$|e|$ [58], or $q_{S}$=+1.0817$|e|$ and $q_{O_{3}}$=–1.852$|e|$ [57]. All of these do not fully compensate the hydronium charge. The residual charge $\Delta q$=$q_{S}+q_{O_{3}}-|e|$ is thus distributed between the other sidechain monomers. In Ref. [56] with $q_{S}$=+1.19$|e|$ and $q_{O_{3}}$=–2.19$|e|$, $\Delta q$=0, and thus the hydronium charge is fully compensated by the charge of the terminal groups. In this sense, our choice for $q_{S}$ and $q_{O_{3}}$ is close to the terminal group charges adopted in Ref. [56].

In our model the total negative charge of the $N_{s}$ terminal groups SO${}_{3}^{-}$ is compensated by the positive charge of the $N_{s}$ hydronium ions modeled as spherical blobs of diameter $\sigma$. The pore also holds $\lambda N_{s}$ water molecules H2O which hydrate the terminal groups and liberate the ions. In bulk Nafion-like ionomers the parameter $\lambda$ is varied between 3 and 25.

The pore geometry parameters $L$ and $d$, and the number of grafting points $N_{s}$ define the grafting density of the sidechains $n_{s}$=$N_{s}/(\pi dL)$, which can be also referred to as the surface charge density of the pore. The average neighbor-to-neighbor distance between the anchoring points is

[TABLE]

For short protrusion lengths with $L_{s}$=2 monomers, when terminal groups are in the vicinity of the pore walls, the parameter $r_{ss}$ essentially has the meaning of the average distance $a_{1}$ between the neighboring sulfonate groups. For Nafion-like ionomers both these quantities usually take values from 0.6 to 1.2 nm. However, for long sidechains with L_{s}$$\geq3 monomers, when the sulfonate groups tend to protrude into the pore volume and form compact clusters, $a_{1}$ is defined from the position of the first maximum of the pair distribution function of sulfonates $g_{OO}(r)$ (here $O$ represents terminal group oxygens O3).

In general, the smaller $a_{1}$ is, the more robust is proton hopping among neighboring sulfonate groups via intermediate water molecules. Considering that the minimal value of $a_{1}$ may not be smaller than 0.6 nm because of the electrostatic repulsion between negative SO${}_{3}^{-}$ groups, at least one water molecule should be present between neighboring sulfonate groups. This should enable accepting a proton from one sulfonate and passing it to another sulfonate [62, 63, 64, 65, 66, 67].

3 Simulation Details

We employ a coarse-grained approach for the sidechain in the framework of the united-atom representation for the CF2 groups (from the blob number 1 to the blob number $L_{s}$-2 of the sidechain) [61, 68, 69], and for the sulphur atom S and the oxygen group O3 of sulfonates [69, 70]. All united-atom groups are modeled as Lennard-Jones (LJ) monomers with a diameter $\sigma$=0.35 nm and 6-12 LJ interactions among them. The sidechain constituents are additionally subjected to stretching, bending and dihedral forces. The force-field details are given in our previous work [55, 69, 71] and agree in most instances with the Nafion sidechain model of Paddison [72]. A brief description of the force field includes four components. First, the total potential energy of the sidechain polymers is

[TABLE]

where $(\vec{r}_{1},\vec{r}_{2},...,\vec{r}_{N})$ are the three-dimensional position vectors of the $N$ particles in the system, the index $i$ in the two-body bond-stretching potential $U_{b}^{i}$ runs over all bonds, the index $j$ in the three-body angle-bending potential $U_{\theta}^{j}$ runs over all bond angles, the index $m$ in the four-body dihedral component of the interaction energy runs over all torsional angles, and indices $k$, $l$ in the non-bonded (Lennard-Jones and Coulomb) potential run over all force-center pairs in the system.

With the sidechain equilibrium bond length $b$=0.44$\sigma$ and the equilibrium bond angle $\theta$=110o, the fully stretched sidechain with a dihedral angle $\alpha$=0 has a physical length $l_{s}$=$L_{s}b\sin(\theta/2)$. For the maximal number of sidechain monomers $L_{s}$=18 we get $l_{s}$=2.27 nm.

Second, the water is modeled as a TIP3P liquid [69, 73], which has explicit charges $q_{H}$ = +0.417$|e|$ on the hydrogen atom, and $q_{O}$ = -2$q_{H}$ on the oxygen atom. The distance between hydrogen and oxygen atoms is $r_{OH}$=0.0957 nm, and the angle between OH bonds is $\theta_{HOH}$=104.52o. The TIP3P model does not allow for inclusion of the non-classical Grotthuss transport of protons [42, 74, 75]. The ion transport in our simulation model represent a classical en-masse diffusion of ions which is affected by the dynamical rearrangement of flexible sidechains [49]. Third, standard Lorentz-Berthelot combining rules are used in the LJ interactions between monomers $i$ and $j$ such that the interaction potential between units $i$ and $j$ has a minimum at the separation distance $r_{ij}=2^{1/6}\sigma_{ij}$ where $\sigma_{ij}$ is defined as $\sigma_{ij}=\left(\sigma_{ii}+\sigma_{jj}\right)$/2, and the depth of the interaction at this minimum is $\epsilon_{ij}=\sqrt{\epsilon_{ii}\epsilon_{jj}}$. For hydrophobic monomers $\epsilon_{ii}$=0.2 kcal/mol, and for hydrophilic monomers $\epsilon_{jj}$=0.1 kcal/mol coefficients were used in the LJ interaction potentials.

The numerical part of our study consists of three series of simulation runs. In the first series of runs, a pore with grafted sidechains is created and the system is subsequently equilibrated during runs of 1 ns at constant volume and temperature. In these runs the charges of the ions, water atoms and terminal groups were set to zero, and the equilibration was done following the hydrophobic and hydrophilic interactions in the system. In the second series of runs, the charges were turned on, and the system was equilibrated for another 1 ns of simulation time under a constant pressure condition $P$=1 atm. The details of this procedure are explained in section 4. Finally, in the third series of runs, the necessary statistics on the sidechain conformations and the ion diffusion were gathered during the final runs of 50 ns duration.

During the production runs, the simulated system was coupled to a Langevin thermostat with a friction coefficient $\gamma=2$ ps*-1* and a Gaussian white-noise force of strength $6k_{B}T\gamma$. The equations of motion were integrated using the velocity Verlet algorithm with a time step of 0.5 fs. One dimensional periodic boundary conditions were imposed on the system along the pore long axis, which coincides with the $x$-axis of the system. The translational replicas of the fundamental cell pore correspond to an infinite cylindrical pore. The long-range electrostatic interactions between charged particles are handled using a one-dimensional Lekner-like summation [15, 76, 77, 78].

For simplicity, it is assumed that there is no dielectric discontinuity between the interior of the pore and the wall material. Otherwise, the interaction of the ions with their image charges will add more complexity to our simulations.

The wall-monomer interaction is treated using the particle-micropore interaction potential given by Tjatjopoulos et al. in Ref. [79], further developed in Ref. [80], and described in detail in our previous work [51]. The wall is considered smooth with no atomic structure on it [19]. The wall-ionomer interaction is hydrophobic for the hydrophobic segments of the sidechain. The wall-terminal group, wall-ion, and wall-water molecule interactions are hydrophilic.

For the dielectric constant of the pore interior we use a distance-dependent effective permittivity function $\epsilon(r)=1+\left(\epsilon_{B}-1\right)\left(1-r/\sigma\right)^{10}/\left(1+(r/\sigma)^{10}\right)$ similar to Refs. [51, 69]. This function gradually increases from $\epsilon=1$ at small monomer-monomer separations to the bulk ionomer value $\epsilon_{B}$=4

for $r/\sigma\gg$1.

Although the best approximation for $\epsilon(r)$ can be a matter of debate [14, 15, 81], it is clear that the permittivity must increase with distance to account for the polarization effects of the neutral and hydrophobic parts of the sidechain when the interaction between the ions is calculated. The suggested $\epsilon(r)$ does not overestimate the Coulomb interactions in the system when explicit water exists in the pore. All water confinement, polarization, and screening effects are explicitly taken into account in our set-up through the explicit ion-water and water-water Coulomb interactions.

4 Simulation Results

In total, there are five system parameters, namely the sidechain protrusion length $L_{s}$ monomers, the pore diameter $d$, the water content $\lambda$, the pore length $L$, and the number of grafting points $N_{s}$ on the pore surface, which fully define the pore structure. The fixing of the last two parameters, the pore length to $L$=46$\sigma$ ($L$=16 nm) and the number of grafting sites to $N_{s}$=200, allows us to reduce the number of system variables from five to three. Subsequently, we will analyze the following three separate cases. The case A with ($\overline{L_{s}},d,\lambda$), the case B with ($L_{s},\overline{d},\lambda$), and the case C with ($L_{s},d,\overline{\lambda}$), where the overlined variables will be kept constant while the two other parameters inside the parentheses will be treated as running variables. The set-up parameters for the cases A, B, and C are given in Table 1.

The pore parameters $d$ and $\lambda$ listed in Table 1 are coupled to each other by the constant pressure condition $P$=1.0 atm. For the one-component system in restricted geometries the pressure can be calculated using the virial theorem [52, 82]. However, this method is not efficient for multi-component systems like confined and hydrated ionomers. In these complex systems the water mostly forms hydration shells around sulfonate groups and hydronium ions where its density is about 10-20 percent higher than the bulk water [83]. As a result of this, the local pressure deviates strongly from the average pressure in the system [84], and thus the virial pressure method becomes less accurate.

We calculate the pressure $P$ from the cumulative force $F_{wall}$ of the hydronium-wall and water-wall interactions. For the cases A and C the pressure $P$=$F_{wall}/(\pi dL)$=1 atm was achieved by treating the pore size $d$ as a running parameter. For the case B all runs were started with the minimal water content $\lambda$=1. This parameter subsequently was increased through adding more water molecules until $P$=1 atm is reached in the pore.

It should be noted that the calculated water content $\lambda$ in the case B, and the pore size $d$ in the cases A and C under the condition $P$=1 atm might be different from experimental realizations because of the non-polarizable water model used in our simulations. The precise polarizable water models [85] might result in different water densities in the solvation shells of ions and terminal groups. This is still an open question that deserves close attention, though it is out of the scope of our current model.

The case A, where the dependence of the ion diffusion on the pore size at various water contents will be elaborated, closely mimics the pore swelling process in ionomer membranes. In the case B the principal role of the sidechain protrusion length $L_{s}$ will be revealed in pores with a fixed size $d$. In this case the change of $L_{s}$ is complemented by the change in the water content $\lambda$. Lastly, the case C is devoted to the role of the pore size $d$ in ion diffusion. For each of these cases we will determine optimal pore parameters for obtaining maximal ion diffusion coefficients $D$ and $D_{x}$, which are defined below.

A snapshot picture taken from the simulation cell for run B4 is shown in Figure 3. Here, along with the sidechains drawn in blue and the sulfonate groups SO${}_{3}^{-}$ shown as spheres, the pore also accommodates $N_{s}$ monovalent ions and $\lambda N_{s}$ water molecules, which are not shown in this figure. Here two different types of effect take place. The sidechain-specific effects are their strong electrostatic fields, which polarize the water and attract

hydronium H3O+ ions

to the sulfonates, and their affinity to clustering due to the dipolar nature of SO${}_{3}^{-}$ groups and mutual SO${}_{3}^{-}$–SO${}_{3}^{-}$ attraction mediated by the ions. The pore-specific effects are the pore wall hydrophobicity and the pore size which regulate the distribution of terminal groups, water molecules and ions.

For each simulation run listed in Table 1, the ion self-diffusion coefficients $D$ and $D_{x}$ were calculated using their mean square of displacements,

[TABLE]

where $r_{i}(t)$ and $x_{i}(t)$ are the 3D and 1D-axial positions of the atom $i$ at the simulation time $t$. Eq.(3) defines only the en-masse (vehicular) ion diffusion. The structural proton diffusion through Grotthuss hopping, which is not included in our simulation model, but which can be modeled using empirical valence bond (EVB) based water models [67, 85, 86], would increase the total diffusion rate of ions.

It should be noted that the pore diameters considered in Table 1 are typical for Nafion-like ionomers, and are larger than the size $\approx$0.4–1.1 nm of narrow carbon nanotubes. Thus, for the run parameters listed in Table 1 we do not enter the regime of single-file diffusion [87, 88, 89]. Moreover, the presence of the ions and grafted chains with hydrophilic ends is a suppressing factor to any occasional movement of hydrogen-bonded water clusters along the pore axis [90, 91].

4.1 Case A simulations for ($\overline{L}_{s}$,$d$,$\lambda$)

The runs A1-A7 from Table 1 are consistent with the ionomer pore swelling when excess water is absorbed by the PEM membrane. That is why we consider shorter sidechains with $L_{s}$=4 monomers which correspond to the sidechain protrusion length about 3–4 Å  in the pores of Nafion-like membranes [14].

We first plot the pair distribution functions $g_{OO}(r)$ of the terminal group oxygens in Figure 4a. The pair distribution functions are defined as

[TABLE]

where $\vec{r}_{i}^{(j)}$ denotes the position of the particle $i$ of the type $j$, $j$=O for the oxygen groups of sulfonates, $\delta(r)$ is the Dirac function, and $\theta(r)$ is the Heaviside step function.

The two maxima of the $g_{OO}(r)$ at separation distances a_{1}$$\approx0.45 nm and a_{2}$$\approx0.7 nm correspond to the first and second nearest neighbor sulfonate groups, respectively. The neighboring sulfonate groups share their hydronium ions in the sulfonate clusters of size 2 nm. In swollen pores $g_{OO}(a_{1})$ decreases whereas $g_{OO}(a_{2})$ increases, which is different from the case of the bulk ionomer, where both these quantities decrease at high solvation levels [37, 92]. This behavior is a result of the increase in the sidechain anchoring distance $r_{ss}$ which makes it harder for the terminal groups to come closer and share a hydronium ion at a separation distance $a_{1}$. As a consequence, in swollen pores the sidechain tips mostly share their counterions at the separation distance $a_{2}$. For larger $\lambda$ the discrepancy between $r_{ss}$ and $a_{2}$ increases, which indicates that the sidechains stretch to form sulfonate clusters.

The water-water pair distribution functions $g_{ww}(r)$, defined by the Eq.(4) with $j$=$w$ for the water oxygens, are shown in Figure 4b. For the low water contents \lambda$$\leq3 the most expected water-water separation distance is about 0.5 nm, which, more likely, characterizes the separation distance between water molecules that share the same ion or sulfonate group. At high water contents $\lambda$, hydration shells are formed around the charged entities together with the water clusters being formed in the pore center. Both these effects appear as the formation of a huge maximum at the touching water-water separation distance $a_{1}$= 0.35 nm for \lambda$$\geq6 in Figure 4b. The position of the second maximum of $g_{ww}(r)$ at a_{2}$$\approx0.6 nm is related to the radius of the second solvation layer in water clusters.

The ion-sulfonate association parameter, measured as the height of the SO${}_{3}^{-}$-ion pair distribution function $g_{SO_{3}-H}(r)$ at $r$=$\sigma$, is plotted in Figure 5a. This parameter holds information on how strongly the electrostatic field of sulfonate groups disturbs the distribution of hydronium ions in the pore. We see that such disturbance increases in swollen pores. This tendency, discussed in Refs. [13, 65], can be justified by the enhancement of the sulfonate-ion interaction

[TABLE]

in wide pores. In Eq.(5) $\epsilon_{s}$ is the permittivity of the medium and includes polarization effects of water in swollen pores, $r_{D}$ is the Debye screening length for the sulfonate groups and hydronium ions,

[TABLE]

where $\epsilon_{0}$ is the permittivity of free space, $n_{e}$ is the concentration of all charges (SO${}_{3}^{-}$ and HsO+) in the system and is defined as $n_{e}$=$2N_{s}/(\pi(d/2)^{2}L)$=$8N_{s}/(\pi d^{2}L)$=$8n_{s}/d$. Using Eq.(6) we get $r_{D1}$=0.028 nm in the unswollen pore of run A1, and $r_{D2}$=0.115 nm in the fully swollen pore of run A7 by assuming that $\epsilon_{s}$=5 in the first, and $\epsilon_{s}$=25 in the second cases. Then, for the relative change of the association strength $U_{2}(\sigma)/U_{1}(\sigma)$ we get U_{2}/U_{1}$$\gg1 which is a clear indication of the fact that in swollen pores the distribution of hydroniums is strongly disturbed by the electric field of terminal groups. This is exactly what is observed in Figure 5a, which, however, should not be interpreted as the localization effect of H3O+ ions around the sulfonate groups.

The localization or delocalization effects of the ions around their host terminal groups can be evaluated through their coordination number (CN)

[TABLE]

shown in Figure 5b. Here $\rho_{H}^{0}$=$N_{S}/V$ is the average sulfonate density in the pore of volume $V=\pi d^{2}L/4$, and $R_{m}$ is the position of the first minimum of $g_{SO_{3}-H}(r)$. The decrease of CN in swollen pores, as seen from Figure 5b, corresponds to the expected ion delocalization effect due to the water screening effects.

Figure 5 also shows the weakening of the ion-sulfonate association $g_{SO_{3}-H}(\sigma)$ when the pore surface charge density is increased. This is seen, for example, for $d$=2.3 nm, where the $g_{SO_{3}-H}(\sigma)$ is about 38 at $\lambda$=10 and $N_{s}$=100, compared to its value about 28 at $\lambda$=3 and $N_{s}$=200. This ion-sulfonate association weakening can be explained, again, using the screening arguments given in Eq.(6). In the pores with higher $N_{s}$ and lower $\lambda$, the parameters $r_{ss}$ and $\epsilon_{s}$, and therefore the screening length $r_{D}$, become smaller. At reduced $r_{D}$ the charges of all ions are strongly screened, and thus the opposite charges in the pore are loosely associated. As a result of this, the CN of the sulfonate groups also decreases in pores with lower surface charge densities, as seen in Figure 5b. In total, the competition between the increased influence of the sulfonate groups on the distribution of the hydronium ions and the increased delocalization of the hydronium ions in swollen pores determines the optimal pore diameter when the hydronium diffusion along the pore axis reaches its maximum value.

The radial distribution functions, defined as

[TABLE]

are shown in Figure 6 for the sulfonate oxygens ($j$=$O$),

hydronium

ions ($j$=$H$), and the water molecules ($j$=$w$). In Eq.(8) $r$ is the radial distance from the pore axis ($r$=0 corresponds to the pore center), and $r_{i}^{(j)}$ is the radial distance of the particle $i$ of the sort $j$. In Figure 6 the radial distributions are normalized by the average density of particles $j$ in the pore. The ion profiles $\rho_{H}(r)$ appears to be completely restricted to the pore boundary regions and no free ions penetrate to the pore center. This result differs from the ion radial distributions calculated using approaches based on Poisson-Boltzmann theory in charged pores [13, 25, 26]. In the latter the ion profile is non-zero in the whole pore volume.

The exclusion of ions from the pore center in Figure 6 effectively means that in charged pores with shorter grafted sidechains, the Grotthuss mechanism of structural proton transfer is restricted to the pore boundary areas.

The radial distribution of water molecules $\rho_{w}(r)$ shown in Figure 6 displays strong water structuring effects in the pore. First, in narrow pores with \lambda$$\leq6, the water mostly accumulates in the hydration shells of ions and sulfonates where it is strongly polarized. This is a consequence of the high surface-to-volume ratio effect in narrow pores, which works against the water clustering in the pore center. In the larger pores the water shows bulk-like properties in the pore central area.

The increase in the height of the maxima in $\rho_{O}(r)$ in Figure 6 in larger pores is related to the frustration of sidechains. When $r_{ss}$ becomes larger in swollen pores, the sidechains start to stretch out so that their tips can share free ions. Also, in larger pores the sidechain-pore wall interaction becomes stronger: while the hydrophobic part of the sidechain becomes more attracted to the pore wall, the hydrophilic tip segment of the sidechain becomes more repelled from the pore wall. An increased frustration of sidechains in larger pores, together with the enhanced sulfonate-ion association effect shown in Figure 5a, brings the ions to the vicinity of the pore walls. This is seen from the increase in the height of the maxima in $\rho_{H}(r)$ in Figure 6. As a result of these frustration-localization effects the flexibility of sidechains in larger pores degrades.

The orientation of the dipole moment of the SO${}_{3}^{-}$ head group, characterized by the function

[TABLE]

is shown in Figure 7, as a function of the water content $\lambda$. In Eq.(9) $\theta$ is the angle between the SO${}_{3}^{-}$ dipole moment and the long axis x of the pore, and $F$=–0.5 is expected for a limiting case when all the SO${}_{3}^{-}$ dipoles point towards the pore center. The orientation parameter $F$ shows a steep decrease when $\lambda$ is increased from 1 to 6, which is followed by a weak change in $F(\lambda)$ for \lambda$$\geq6. Such behavior can be explained by assuming that at $\lambda$=6 all ions develop a compact hydration shell which is strongly repelled from the walls of the larger pores. This results in the reorientation of the sulfonate dipoles towards the pore center. On the other hand, the accumulation of bulk water in the pore center has only a negligible effect on the sulfonate dipole reorientation.

The ion self-diffusion coefficients $D$ for 3D ion diffusion in the pore, and $D_{x}$ for 1D ion diffusion along the pore axis $x$, defined by Eq.(3), are shown in Figure 8 as a function of the pore size $d$. For all the hydration levels for the runs A1-A7, simulation results show $D_{x}>D$, showing that the ion diffusion along the pore axis is higher than the total ion diffusion in the pore [93]. While the 3D diffusion has almost a monotonic dependence on $\lambda$, the 1D diffusion reaches a maximum at $\lambda$=6. Such nonlinear behavior of the $D_{x}(d)$ is a consequence of the following two factors: First, as the pore swells, the increase in the grafting parameter $r_{ss}$ and the decrease in the sidechain flexibility constrain the ion diffusion along the pore. Second, as the pore shrinks, the low amount of water, insufficient to hydrate the ions, becomes highly structured and polarized by the charges. As a consequence, the ions become immobile in the structured water environment. The maximum of the ion diffusion, hence, happens in between these two extreme too-large and too-narrow pore cases. Figure 8 also reveals that the surface charges of the pores facilitate ion diffusion in the pore. This is because of the electrostatic delocalization effect of ions due to the weakening of the ion-sulfonate association at high $n_{s}$. In total, in the pores with shorter sidechains the ion diffusion is limited to the pore wall area, and its component along the pore axis is maximal if the electrostatically delocalized ion possesses a hydration shell with about six water molecules.

It is known that at sufficiently longer times the diffusion coefficient of the ions in the radial direction of the pore becomes very small, D_{r}$$\to0, thus the full diffusion coefficient

[TABLE]

has an asymptotic limit $D_{x}/D$=3. In our simulations this limit has been reached in the narrow pores of runs A1 and A2 for $N_{s}$=100, and in the narrow pore of the run A1 for $N_{s}$=200. In the larger pores of runs A3-A7 the ratio $D_{x}/D$ decreases from 3 to around 1.6–1.8. Therefore, runs much longer than 50 ns are needed to approach the asymptotic limit $D_{x}/D$=3.

4.2 Case B simulations for ($L_{s}$,$\overline{d}$,$\lambda$)

In this section we focus on the role of the sidechain protrusion length $L_{s}$ in providing fast ion diffusion along the pore axis. The runs B1-B7 from Table 1 were carried out at the constant pore size $d$=2.89 nm. The chain protrusion length $L_{s}$ was increased from 2 monomers to 16 monomers with a simultaneous decrease of the water content $\lambda$ from about 9 to 1. The shortest $L_{s}$=2 monomers corresponds to a completely hydrophilic sidechain consisting of only the two charged entities of the terminal group SO${}_{3}^{-}$. The tips of the longest sidechain with $L_{s}$=16 monomers can reach the pore center and thus are capable of forming hydrophilic sulfonate clusters there. The results of the previous section indicated that the ions mostly occupy the pore wall areas in the case of shorter sidechains $L_{s}$. For the longer $L_{s}$ it is expected that the ions would reach the pore center by following their host sulfonate groups. In both cases the ion diffusion depends on the state of water in the pore. In other words, a priori, it is not evident at which $L_{s}$ and $\lambda$ the ion diffusion will be maximal.

The sulfonate-sulfonate pair distribution functions $g_{00}(r)$ shown in Figure 9a clearly indicate an enhancement of the sulfonate clustering as the sidechain protrusion $L_{s}$ is increased. This result is fully expected, because long sidechains have more flexibility to form clusters in the pore volume. As the clusters grow in size, more ions and water molecules are shared between the neighboring sulfonate groups. This explains the shift of the position of the second maximum of $g_{00}(r)$, $a_{2}$, to the larger sulfonate-sulfonate separations. Since long sidechain protrusion lengths $L_{s}$ also mean lower values for $\lambda$ in the pore, it is obvious that sulfonate clustering features are essentially regulated by the protrusion length $L_{s}$ rather than by the water content $\lambda$.

The water pair distribution functions $g_{ww}(r)$ also show higher water cluster formations for the long sidechains, as seen in Figure 9b. Even if the long sidechains demand less water to fill the pore, it seems that all the available water is collected in the pore central area where it easily merges into larger clusters. This conclusion is supported by the radial distribution profiles of ions, water molecules and sulfonates, shown in Figure 10. Indeed, as the protrusion length $L_{s}$ increases, most of the water molecules cluster in the pore center, perhaps by being attracted there by the sulfonate tips of the protruding sidechains and the ions. For the sidechain monomer number $L_{s}\leq$8 the radial distributions in Figure 10a and Figure 10b look qualitatively the same as the densities shown in Figure 6 for \lambda$$\geq6: the water mostly occupies the pore central area, and the ions mostly lie between the water and sulfonate clusters. However, for the long chains with the number of monomers $L_{s}\geq$12, when the sulfonate groups reach the pore center, the ions have their maximal density on the $x$-axis. As seen from Figure 10c, in this region the difference between the radial densities of the sulfonate oxygens $\rho_{O}(r)$ and ions $\rho_{H}(r)$ gives rise to the charge separation in the pore. The central part of the pore has a positive charge whereas the rest of the pore is charged negatively. This result is the direct consequence of the long sidechains: a pore with the long sidechains, when the chains partly occupy the pore wall area, behaves like a narrow pore with shorter sidechains.

In narrow pores, when their radius is comparable with the sidechain length, the counterions enter the pore center generating a charge separation between the central and wall areas of the pore.

Their diffusion along the pore axis, however, depends on the sulfonate clustering structure in that area.

The $x$-axis profile of sulfonate clustering, defined as

[TABLE]

where $x_{i}$ is the $x$-coordinate of sulfonates, is shown in Figure 11 for the run B6 with $L_{s}$=12 monomers (the case of radial charge separation in the pore). The sulfonate clusters self-organize into compact formations along the pore axis with the average cluster size r_{c}$$\approx2 nm and a separation distance between them of $R_{c}$=2 nm. For $N_{S}$=200 sidechains along the $L$=16 nm long pore, a rough estimate gives about N_{c}$$\approx40–50 sulfonates per cluster formation.

Calculated ion diffusion coefficients $D$ and $D_{x}$ are shown in Figure 12 for the pore size $d$=2.89 nm and the sidechain anchoring distance $r_{ss}$=0.8 nm. Both diffusion coefficients show a nonmonotonic dependence on the sidechain protrusion length $L_{s}$. The maximal ion diffusion is observed for the run B4 with $\lambda$=5.7 and $L_{s}$=8 monomers. According to the findings of Paddison [14], who calculated ion diffusion in the framework of non-equilibrium statistical mechanics, the ion diffusion should decrease when the protrusion length of the anchored sidechains is increased. We believe that the difference between our findings and the results of Ref. [14] is based on the fact that we explicitly treat the dynamics and clustering behavior of the water, ions, and the sulfonates, which are not included in the statistical mechanical theory of Ref. [14]. The fact that the ion diffusion for the run B7 with $L_{s}$=16 monomers is very small compared to its value for the run B4 with $L_{s}$=8 monomers, implies that the clustering of sulfonates and ions in the pore center is a degrading factor for the free ion diffusing in the pore. Therefore, for getting faster ion diffusion, it seems that the central pore area should be free from ions and sulfonates. At the same time, similar to the case A runs in the previous section, the ions should be partly delocalized from their host sulfonates and moderately hydrated.

4.3 Case C simulations for ($L_{s}$,$d$,$\overline{\lambda}$)

In the case A and B simulations discussed in previous sections, the variation of $L_{s}$ and $d$ were accompanied with change in the water content $\lambda$. In the current section we fix the water content to $\lambda$=8, the average hydration number for Nafion-like ionomers, and focus on finding an enhanced ion diffusion in the charged pores. We choose the run B2 as a reference run for the case C simulations, and vary the sidechain protrusion length $L_{s}$ from 4 monomers to 16 monomers. Simulation results for the radial density distribution $\rho_{j}(r)$ of the ions ($j$=H), sulfonate oxygens ($j$=O), and the water molecules ($j$=w), shown in Figure 13, resemble the radial distributions for the case B runs in Figure 10 [94]. Here, again, in the case of long sidechains with $L_{s}$=16 monomers, the ions are electrostatically delocalized in the central pore area and generate a radial charge separation in the pore. A remarkable fact is the flattening of the ion and water profiles in the central pore area in Figure 13c at $\lambda$=8, in contrast to the steep gradients of the $\rho_{H}(r)$ and $\rho_{w}(r)$ in the same area as seen in Figure 10c for $\lambda$=2.8. At the same time, no such flattening of the sulfonate profile in Figure 13c from the increased water content is visible. As a result of these observations, we conclude that the radial charge separation in the pore decreases if more water is absorbed in the pore.

From the similarity of the ion and sulfonate radial distribution profiles in Figure 13b and Figure 10b it is more likely to expect that the run C3 with $L_{s}$=10 monomers and $d$=3.3 nm would provide maximal ion diffusion in the pore. The 3D distribution of ions for some of the case C runs are plotted in Figure 14.

The averaging of the ion positions in these plots was done during the 50 fs long simulations in fully equilibrated system to get instantaneous ion distribution profiles in the pore.

It is seen that for a shorter sidechains with L_{s}$$\leq10 monomers, the ions mostly accumulate near the pore walls, where they form continuous hollow-cylinder-like ionic pathways. The hollow area of this structure is filled with water molecules. For the sidechains with bigger monomer numbers $L_{s}$, the ions accumulate at the pore center by forming larger clusters with disrupted connections along the pore axis. These disrupted and hollow areas are filled either with water molecules or with sulfonates, with the latter hindering smooth ion passage through the pore.

Calculated diffusion coefficients for the ions as a function of the sidechain protrusion length $L_{s}$ for the runs C1-C7 are plotted in Figure 15. The axial diffusion coefficient $D_{x}$ for the entire range of the sidechain protrusion length $L_{s}$ variation is more than twice as large as the full diffusion coefficient $D$ of the ions, a result of the strongly suppressed ion diffusion in the radial direction, perpendicular to the pore walls. A maximum in the $D_{x}$ is seen for the sidechain protrusion length $L_{s}$=10 monomers (run C3), as is expected from the radial distributions of ions and sulfonates in Figure 14. As will be discussed in the next section, it is indeed possible the predict the pore parameters at which maximal ion diffusion along the pore axis is expected.

4.4 Maximal ion diffusion in charged pores

In this section we collect and discuss the maximal ion diffusion coefficients $D_{x}^{max}$ along the pore axis, detected for the runs A3, B4, and C3. The dependence of the $D_{x}^{max}$ on the sidechain protrusion length $L_{s}$ is plotted in Figure 16a. It appears that $D_{x}^{max}$ is a linear function of the sidechain protrusion length $L_{s}$. Taking into account the fact that long sidechains in larger pores behave like short sidechains in narrow pores, it seems plausible to find a scaling rule for the pore parameters $d$ and $L_{s}$ at which the ion diffusion is maximal along the pore axis. For the average sidechain extension into the pore center it is reasonable to introduce the gyration radius of the sidechain

[TABLE]

with a phenomenological Flory-like exponent $\nu$. Here $b$ is the segment length of the sidechain, $b$=0.154 nm for Nafion-like ionomers, For the polymer chains with excluded-volume interactions the exponent varies between 1/2 and 3/5 [95], whereas for attractive coil-globule collapsing polymers the Flory exponent is 1/3. A lower value $\nu$=1/4 has been reported for dilute branched polymers [96, 97]. Twice the gyration radius can be considered as the protrusion length of the sidechain into the pore center, $l_{s}$=$2R_{g}$.

In Figure 16b we test a phenomenological fit

[TABLE]

for the pore parameters $d$ and $L_{s}$ at which the ion diffusion is maximal. The scaling parameter $\Gamma$ measures the ratio of the sidechain protrusion length $l_{s}$ to the pore radius $d/2$. We found that for the fitting exponent $\nu$=1/4 the pore scaling parameter is almost constant at \Gamma$$\approx0.33 for $L_{s}$ changed from 4 monomers to 10 monomers, and $d$ changed from 2.7 nm to 3.3 nm. The fitting value $\nu$=1/4 is the same as the Flory exponent for branched polymers. This is a reasonable conclusion, because the sidechains in the pore create clusters which indeed can be viewed as a branched polymer. The result \Gamma$$\approx0.33 implies that the sidechain extension $l_{s}$ into the pore center should be about 1/3 of the pore radius for the pore to provide maximal ion diffusion along the pore axis.

An interesting question is how the optimal protrusion length $l_{s}$ depends on the diffusing ion size $\sigma_{i}$. Apparently there is no direct link between these two quantities. Under the condition that the free volume for the ions in the pore is kept constant, the change in $\sigma_{i}$ will affect the ion diffusion rates in all simulated cases A, B, and C in the same manner. Smaller ions will be more attracted to the oppositely charged sulfonates which will slow down their diffusion, whereas larger ions will be more readily released from the sulfonates and thus have higher diffusion rates. Therefore, we think that the diffusion lines shown in Figures 8, 12, 15, and 16a will retain their shape but shift either up or down depending on the decrease or increase in $\sigma_{i}$ respectively.

The maximal diffusion coefficients $D_{x}^{max}$ shown in Figure 16 are in the range D_{x}^{max}$$\approx(10–14)$\times$10*-6cm2/sec for the sidechain lengths L_{s}$$\approx8–10 monomers. These values are about half the size of the self-diffusion of water molecules in bulk water D_{H_{2}O}$$\approx23$\times$10-6cm2/sec [63, 64, 93, 98, 99, 100], and about 30–40% smaller than the simulated vehicular (en-masse) diffusion values for hydronium ion in bulk water D_{H_{3}O}$$\approx(17–20)$\times$10-6cm2/sec [59, 64]. Obviously, it is more relevant to compare our data for $D_{x}^{max}$ with the water and hydronium diffusion coefficients in ordinary ionomers with L_{s}$$\approx8 monomers like the Nafion membrane. In this case, however, one has to distinguish between the local $D^{l}$ and long-range $D^{lr}$ diffusion coefficients, with the $D^{l}$ corresponding to the particle diffusion in a single cluster at shorter time scales, and $D^{lr}$ corresponding to the particle diffusion at longer time scales when the particle has enough time to pass through several hydrophilic clusters in the ionomer. Generally, $D^{l}\gg D^{lr}$ because of the tortuousity and narrow bridges between neighboring hydrophilic clusters in the ionomer [101]. For example, in a fully hydrated ionomer with \lambda$$\approx15–17 the water molecule has a local self-diffusion coefficient D_{H_{2}O}^{l}$$\approx13$\times$10-6cm2/sec which is about twice as large as its long range self-diffusion coefficient D_{H_{2}O}^{lr}$$\approx6$\times$10-6cm2/sec [102], and much higher than its long range self-diffusion coefficient $D_{H_{2}O}^{lr}$(\lambda$$\approx5–10)$\approx$(1–5)$\times$10-6cm2/sec in partly hydrated Nafion [42, 103]. Additionally, one has to bear in mind that experimentally determined diffusion coefficients for hydronium ions include the Grotthuss proton hopping mechanism, which is not accounted for in our simulations. For example, $D_{H_{3}O}^{lr}$ in fully hydrated Nafion with the Grotthuss hopping accounted for reaches $D_{H_{3}O}^{lr}$(Grotthuss + en-masse)$\approx$(12–14)$\times$10-6cm2/sec [59] which is more than twice as large as the en-masse diffusion $D_{H_{3}O}^{lr}$(en-masse)$\approx$(4–5)$\times$10-6cm2/sec [59, 104]. The same is valid in partly hydrated Nafion with \lambda$$\approx6–10, where $D_{H_{3}O}^{lr}$(Grotthuss + en-masse)$\approx$(2–5)$\times$10-6cm2/sec is larger than $D_{H_{3}O}^{lr}$(en-masse)$\approx$(1–2)$\times$10-6*cm2/sec [63, 103].

Our data for $D_{x}^{max}$ corresponds to the long-range diffusion $D^{lr}$ of hydronium ions without the Grotthuss contribution. Therefore, we get a diffusion rate 2–3 times larger than $D_{H_{3}O}^{lr}$(en-masse) in fully hydrated Nafion [59, 104], and about 5–7 times larger in partly hydrated Nafion [63, 103].

5 Conclusions

In this work we analyzed ion diffusion phenomena in charged nanopores grafted with ionomer sidechains. Our aim was to determine the optimal pore parameters for obtaining enhanced ion diffusion.

We found that in the case of short sidechains, the hydronium ions mostly occupy the pore wall area and their distribution is strongly disturbed by the electric field of sulfonates in swollen pores. Additionally, the hydronium ions delocalize from their host sulfonate groups in wide pores. However, an increase in the anchoring distance $r_{ss}$ and associated with it a stretching of sidechains restricts the ion diffusion in swollen pores.

In narrow pores, on the other hand, the structuring and polarization effects of the water molecules hinder the free movement of the ions in the pore. Consequently, the ion diffusion in the charged pores grafted with short sidechains becomes a nonlinear function of the pore diameter. It attains a maximum in the pores with less structured water, high surface charge density, and moderately delocalized ions occupying only the pore wall area.

In the charged pores with a fixed pore size $d$, according to our simulations, the clustering of the sulfonates is strongly regulated by their protrusion length $L_{s}$ rather than by the water content $\lambda$. For the longer $L_{s}$, when the sidechain tips reach the pore center, a radial charge separation occurs in the pore. The pore center with excess ions is charged positively, while the pore wall area with excess sulfonates is charged negatively. Such charge separation, which is associated with the sulfonate clustering in the pore central area, suppresses the ion diffusion along the pore axis. An enhanced ion diffusion was found in the pores grafted with medium-size sidechains provided that the ions do not enter the central pore area, and the water is less structured around the ions. A similar conclusion is also made for the pores with a fixed water content. The medium-sized sidechain with tips not entering the pore central area allows the ions to form hollow-cylinder type hydrophilic pathways in the channel. The existence of such cylindrical shells with smooth and uninterrupted ion pathways is the necessary condition for getting an easy ion passage along the pore axis.

By collecting and analyzing the pore parameters for which our simulations have detected maximal ion diffusion along the pore axis, we found that the ion diffusion $D_{x}^{max}$ has a linear dependence on the sidechain protrusion length $L_{s}$. The long chains, provided that they do not reach the pore center, have very flexible tips which assist the ion diffusion along the pore axis. We also proposed a simple scaling rule for the pore parameters $d$ and $L_{s}$ with a Flory-like exponent 1/4. It appears that a maximal ion diffusion along the pore axis is possible if the effective length of the sidechain extension into the pore center, measured as twice the gyration radius of the sidechain with the Flory-like exponent 1/4, is about 1/3 of the pore radius $d/2$.

The simulated axial diffusion coefficients have maximal values in the range D_{x}^{max}$$\approx(10–14)$\times$10*-6cm2/sec for water content \lambda$$\approx6–8, and sidechain length L_{s}$$\approx8–10. These diffusion rates are about 2–3 times larger than the diffusion rate $D_{H_{3}O}^{lr}$(en-masse) in fully hydrated Nafion [59, 104], and about 5–7 times larger than the diffusion rate $D_{H_{3}O}^{lr}$(en-masse) in partly hydrated Nafion [63, 69, 103]. We believe that experimental realizations of our grafted pore set-up will produce even higher ion diffusion rates because of the Grotthuss structural diffusion contribution to $D_{x}^{max}$. For example, the Grotthuss proton diffusion rate in bulk water is $D_{p}$(Grotthuss)=70$\times$10-6cm2/sec, and with the water self diffusion $D_{p}$(en-masse)=23$\times$10-6cm2/sec added to it, the total proton diffusion rate roughly increases up to $D_{p}$(Grotthuss + en-masse)=93$\times$10-6*cm2/sec [64, 100]. In Molecular Dynamics simulations the Grotthuss hopping can be accounted for by implementing the empirical valence bond method [13, 15, 85].

These results obtained for ion diffusion in charged pores may be helpful in the further development of new emerging technologies, such as energy storage in nanoporous metals [6], or a new class of soft fluid actuators [105]. Our results also might contribute to the development of a new class of elastomer actuators called metallic muscles [106]. In these applications the pores of the metal inclusions are polymer coated and charged by doping with sulfuric acid. Nanocomposite elastomers with these nanoporous metal inclusions, besides generating electrostrictive strains [107, 108, 109], will also generate additional strain from the interaction of the electrostatic double layers between these highly charged objects. Then the total actuation response of the nanocomposite to the applied field will be greatly enhanced.

We also propose that the concept of the self-assembling of amphiphilic diblock copolymers in cylindrical confinements [110, 111] can be used to obtain grafted pores. In this case the free hydrophilic ends of the polymer will form brush conformations [112, 113]. The functionalization and hydration of these brushes will turn the pore into an ion-conducting channel.

Acknowledgments

E. Allahyarov is thankful for support of this work by the Deutsche Forschungsgemeinschaft (DFG) through the grant AL 2058/1-1 on nanocomposite elastomer actuators. H. Löwen thanks the Deutsche Forschungsgemeinschaft for support of the work through the SPP 1681 on magnetic hybrid materials. P. Taylor acknowledges support by the American Chemical Society Petroleum research Fund under Grant PRF#51995-ND7.