Kinetic Monte Carlo simulations of proton conductivity

The kinetic Monte Carlo method is used to model the dynamic properties of proton diffusion in anhydrous proton conductors. The results have been discussed with reference to a two-step process called the Grotthuss mechanism. There is a widespread belief that this mechanism is responsible for fast proton mobility. We showed in detail that the relative frequency of reorientation and diffusion processes is crucial for the conductivity. Moreover, the current dependence on proton concentration has been analyzed. In order to test our microscopic model the proton transport in polymer electrolyte membranes based on benzimidazole ${\mathrm{C}}_7{\mathrm{H}}_6{\mathrm{N}}_2$ molecules is studied.

Figure 1

A possible distribution of protons (solid dots) and allowed movements in the chain. Upper arrows represent rod rotations by angle $\pi$ that may but do not have to lead to a different configuration (this happens if the rotating rod is not occupied or occupied by two protons). Lower arrows represent acceptable hoppings of a proton between the neighboring rods for the particular configuration in the picture. The periodic boundary conditions permit a hop from the rightmost rod to the leftmost one.

Figure 2

A schematic description of the new configuration choice in the KMC algorithm. For the situation presented in the picture, provided that a number $u_1$ will be drawn the system is transformed to configuration number 2.

Figure 3

There are only two configurations requiring the additional factor representing effective Coulomb forces: one presented in the picture above and its mirror reflection. The dashed lines represent H-bond potentials. For the initial configuration above with one proton in each H-bond before the rod's rotation there is one proton and one vacancy in each H-bond which is energetically favorable. After the rotation two protons meet in one H-bond and two vacancies in another.

Figure 4

The log-log dependence of the proton current on the relative frequency ${\gamma }_T/{\gamma }_R$ for the half-filling case $c=0.5$. The individual curves are parametrized by the Coulomb potential $V_{\mathrm{Coul}}$.

Figure 5

The current dependence on the proton concentration including the presence of the Coulomb repulsion for different values of $V_{\mathrm{Coul}}$. The temperature and the external electric field are fixed.

Figure 6

The shape of the potential $V\left(x\right)$ in the presence of an external electric field. The dark shaded area (strips) shows the contribution to integral (15) determined by the value of $E$, the brighter one shows the range of energies contributing to Eq. (13). The hopping from the lower minimum, $V_{-}$, to the upper one, $V_{+}$, introduces the factor $exp[-(V_{+}-V_{-})/\left(k_{B}T\right)]$ to overcome the physically forbidden region for a proton with energy less than $V_{+}$. The energy in Eqs. (13, 14, 15) is measured from the upper minimum, i.e., $V_{+}=0$.

Figure 7

The current flow stabilizes after less than $5\times {10}^{-5}$ s, thus the switch-on effect may be neglected. Inset: the exemplary results for different system sizes converging quickly with the system size (the first two points are for $N=50$ and 100).

Figure 8

Comparison between the measured and simulated data for the benzimidazole.