Insight into hydrogen production through molecular simulation of an electrode-ionomer electrolyte system

TL;DR

This study uses molecular simulations combining first principles and classical dynamics to analyze electrode-ionomer systems at high voltage and pH, aiming to understand factors affecting hydrogen production efficiency.

Contribution

It provides detailed molecular insights into electrode-electrolyte interactions and tests hypotheses from prior experimental work on hydrogen evolution.

Findings

Potential profile changes with bias voltage and oxide coverage.

Alterations in solvation and species concentrations near the electrode.

Water orientation at reactive sites influences hydrogen evolution.

Abstract

In this work, we examine metal electrode-ionomer electrolyte systems at high voltage / negative surface charge and at high pH to assess factors that influence hydrogen production efficiency. We simulate the hydrogen evolution electrode interface investigated experimentally in Bates et al., Journal of Physical Chemistry C, 2015 using a combination of first principles calculations and classical molecular dynamics. With this detailed molecular information, we explore the hypotheses posed in Bates et al. In particular we examine the response of the system to increased bias voltage and oxide coverage in terms of the potential profile, changes in solvation and species concentrations away from the electrode, surface concentrations, and orientation of water at reactive surface sites. We discuss this response in the context of hydrogen production.

Equations

$$-\nabla \cdot (ϵ\nabla \varphi )=\rho .$$

$$\rho =e\sum _a{z}_a{c}_a,$$

$$c_a={\bar{c}}_a\exp \left(\frac{e{z}_a\varphi }{k_{B}T}\right),$$

$${\lambda }_D=\sqrt{\frac{k_{B}Tϵ}{e^2{\sum }_a{z}_a^2{\bar{c}}_a}},$$

$$\varphi (x)=-2\frac{k_{B}T}{e}\ln \left(1+\frac{x}{\ell }\right)+{\varphi }_0$$

$$c(x)=\frac{c_0}{{\left(1+\frac{x}{\ell }\right)}^2},$$

$$\Phi ({\mathbf{x}}_{\alpha })$$

$$+\underset{\text{covalent}}{\underbrace{\sum _I\left[\sum _{\alpha ,\beta \in M_I}{k}_{ab}{r}_{\alpha \beta }^2+\sum _{\alpha ,\beta ,\gamma \in M_I}{k}_{abc}{\theta }_{\alpha \beta \gamma }^2+\sum _{\alpha ,\beta ,\gamma ,\mu \in M_I}{k}_{abcd}{\varphi }_{\alpha \beta \gamma \mu }^2\right]}}$$

$$q_{\alpha }^{\prime }(\sigma )=q_{\alpha }-\frac{\sigma }{e\eta }\Delta q_{\alpha }$$

$$\phi(z)=\phi(0)-\frac{1}{\varepsilon_{0}A}\int_{0}^{z}\overbrace{\bigl{<}\sum_{\alpha|z_{\alpha}\in[0,\tilde{z}]}q_{\alpha}\bigr{>}}^{Q(\tilde{z})}\,\mathrm{d}\tilde{z}\ ,$$

$$\mu(z_{I})=\frac{1}{V_{I}}\Bigl{<}\sum_{\alpha}q_{\alpha}z_{\alpha}w(z_{I}-z_{\alpha})\Bigr{>}-\frac{1}{V_{I}}q_{I}z_{I}\ ,$$

$$q_{I}=\Bigl{<}\sum_{\alpha}q_{\alpha}w(z_{I}-z_{\alpha})\Bigr{>}\ .$$

Full text

Ionomer electrolyte, electrolysis, molecular dynamics

I Introduction

Hydrogen has the capacity to be an ecologically-friendly fuel since water is the primary by-product of its use. Many technological and economic challenges remain in realizing a viable hydrogen economy and energy system Rahman and Andrews (2006); Council (2004). The central issue is that molecular hydrogen gas does not occur naturally in abundance and must be produced industrially. Currently, the majority of hydrogen is generated by either high-temperature/high-energy methane reforming or the water-gas shift reaction which produces significant carbon dioxide, while the electrolysis of water accounts for a relatively minor proportion of its production Navlani-García et al. (2018). As electrolysis is not directly dependent on fossil fuels there is strong motivation to develop a low cost, low energy electrolytic process. Since electrolysis traditionally depends on expensive Pt-group catalysts, transition metal catalysts are being developed.

Mukerjee and co-workers have been developing Ni based catalysts with ionomer-based electrolytes Bates et al. (2015); Ünlü et al. (2011) which show promise but are still confronted by development challenges in part due to the need for a more fundamental understanding of the electrode-electrolyte interaction. In Bates et al. Bates et al. (2015) a number of hypotheses were put forward:

• •

The electrical potential is significantly altered by the ionomer. The ionomer extends the “electrified interface” near the catalyst. The potential at the inner Helmholtz plane dictates the rate of the hydrogen evolution reaction (HER). This potential is more positive than it would be in absence of ionomer and this affects the reactant water molecules at the inner plane.

• •

The water molecules inside the inner Helmholtz plane orient with dipoles pointing away from the electrode due to the presence of the ionomer, and this facilitates H–OH cleavage on the metal catalyst

• •

The ionomer itself straddles the inner and outer Helmholtz planes and is directly chemisorbed on the metal surface.

• •

The majority of the metal surface is coordinated with water molecules.

• •

Nanoscale heterogeneity provides a high density of adjacent metal/metal-oxide sites where metallic Ni has an affinity for H-bonding, and NiOx has an affinity for adsorbed OH-, in-line with Markovic’s theory of enhanced hydrogen evolution reaction on composite metal/metal-oxide surfaces Subbaraman et al. (2011); Strmcnik et al. (2013); Danilovic et al. (2012).

In this work we investigate these hypotheses via molecular simulation. Since a system encompassing the compact/diffuse layers of the long chain ionomer-based electrolyte is too large for ab initio calculations such as Refs. (9; 10; 11; 12; 13; 14; 15; 16; 17), which allow charge transfer and spontaneous dissociation, we employ classical molecular dynamics (MD) to model the system at the relevant length and time scales. Water-metal interfaces have been studied with MD for some time Spohr (1999); Willard et al. (2009); Limmer et al. (2013); Yeh, Janik, and Maranas (2013); Hughes and Walsh (2014); Dewan et al. (2014); Limmer et al. (2015); Takae and Onuki (2015); Foroutan, Darvishi, and Fatemi (2017). The absence of chemical reactions in these simulations is offset by the fact that these reactions are fast compared to diffusion timescales resolvable with MD and we pre-populate the simulation with the relevant chemical species. We assume simulations of the isothermal steady-state accounting for electrostatic and steric interactions with the experimentally observed species concentrations is informative of the transport-limited steady operation of the actual cell. Also, we apply ab initio methods to assess the near electrode charge environment and use this information in the MD model of the electrode. We focus on the hydrogen evolution reaction (HER) environment at the negatively charged electrode (cathode). Here the ionomer is in contact with bare metal (Ni) and metal oxide (NiOx) at the electrode surface. The ionomer-based electrolyte has characteristics that differ from the well-studied dissolved salt electrolytes. For instance, the ionomer is relatively immobile due to its polymeric structure with fixed charge centers (N+), while the counter ions (OH-) are mobile. The influence of the immobile charge in ionomer strands/chains on performance is central to our investigation.

Given its importance, the structure and chemistry of water near a metal surface has long and intense field of study in and of itself, which is reviewed in Refs. (27; 28; 29; 30; 31). The examination of the water-surface interactions has been pursued in great detail by combinations of microscopy and ab initio simulation (typically in vacuum), as in Ref. (32), and has lead to many postulated and observed intact, partially dissociated, and fully dissociated configurations of water at uncharged atomically flat and rough surfaces. Despite the limitations of classical MD, we attempt to interpret the results of this study in the context of this deep body of research.

In Sec.​ II we briefly review the relevant theory. With representative electrode-ionomer systems we simulate the response of these systems to a range of external electrical bias and characterize this response by a variety of means to address the hypotheses in Bates et al. Bates et al. (2015). In Sec.​ III we describe how we use spatial resolution of charge density, per-species radial distributions and other fields to quantify the location of the Helmholtz planes and significant concentrations of the reactive species. The results are given in Sec.​ IV and are discussed in Sec.​ V in light of the hypotheses of Bates et al. Bates et al. (2015)

II Theory

Due to high electric fields near electrodes, charged mobile species tend to pack and form characteristic structures and concentration profiles. This phenomena is modelled in classical theory by Helmholtz Helmholtz (1853), Gouy Gouy (1910), Chapman Chapman (1913), Stern Stern (1924) and others, and is still a topic of current research Merlet et al. (2014). Here we will briefly review the relevant theory to assist in interpreting the molecular results in Sec.​ IV.

The most commonly used theory is the Gouy-Champman-Stern (GCS) model, composed of a compact layer of ions next to the electrode and a diffuse layer beyond. The compact layer is bounded by the inner Helmholtz plane (IHP) at a layer of unsolvated ions adsorbed on the electrode surface and the outer Helmholtz plane (OHP) where ions are fully solvated. The interior of the compact layer is assumed to be charge free due to steric effects. The structure of the diffuse layer is governed by the Poisson-Boltzmann (PB) equation where the electrostatic interaction of the ions is given by the Poisson equation:

[TABLE]

Here, $\phi$ is the electric potential and $\epsilon$ is the electric permittivity. The net charge density $\rho$ is given by:

[TABLE]

where $z_{a}$ is the valence of species $a$, $c_{a}$ is its concentration, and $e$ is the unit of elemental charge. In equilibrium, the absence of fluxes implies that both the electrochemical potential, $ez_{a}\phi+k_{B}T\ln c_{a}$, and the temperature, $T$, are constant across the system domain. This condition leads to species concentrations, $c_{a}$, that vary with the local potential $\phi$:

[TABLE]

where $\bar{c}_{a}$ is the far-field/bulk concentration of species $a$. Substituting Eq. (2), and Eq. (3) in Eq. (1) results in the Poisson-Boltzmann (PB) equation. The characteristic thickness of the diffuse layer is given by the Debye length

[TABLE]

which is the similarity parameter in the solution of the linearized PB equations. For more details see Ref. (38)(Sec.13.3), Ref. (39) (Chap.12), and Ref. (40) (Sec.7.4).

Real electrode-electrolyte systems deviate from GCS for a number of reasons such as: finite ion size and correlation effects due to van der Waals interactions, and a non-uniform dielectric field due to varying concentration and discrete charges Kornyshev (2007); Oldham (2008); Bazant, Storey, and Kornyshev (2011). Also theoretical concepts such as full solvation and the location of Helmholtz planes become less well-defined. In the present case, an electrolyte with a dense ionomer, the charge centers of ionomer chains are effectively immobile and a simple model assumes these ionomer charges provide a background charge density which is fully solvated and screened far from the electrode. It follows that the PB equation governs the excess charge of the counter ion. This paradigm holds where counter ions stay bound to ionomer charge centers, and fails where the electric field is strong enough to dissociate the mobile species from the charge centers. The solution to the PB equation for the counter ions only Gray and Stiles (2018) is:

[TABLE]

and

[TABLE]

where $\phi_{0}$ is the potential at $x=0$, $\ell=2\frac{\epsilon k_{B}T}{e\sigma}$ and $\sigma$ is the surface charge of the electrode. Note that the dielectric of the electrolyte, $\epsilon$, is assumed to be spatially uniform in this derivation. This diffuse layer solution has an electric field $\mathbf{E}=-\boldsymbol{\nabla}\phi$ that is zero far away from the electrode and has the magnitude $\frac{\sigma}{\epsilon}$ at $x=0$, which is consistent with Gauss’s law applied to the electrode (this discussed in more detail in the next section).

III Method

To simulate the metal electrode-ionomer based electrolyte system, we employ a combination of density functional theory (DFT) and classical molecular dynamics (MD). We use DFT to determine the charge state of the metal HER electrode, which is our primary focus. We use MD to model the dynamics of the molecular species of the electrolyte through interplay of electrostatic, elastic, steric and diffusion forces in the overall system. By its explicit representation of atoms and atomic bonding, MD is known to capture deviations from classical theories such as GCS Lee et al. (2012, 2013); Lee, Mani, and Templeton (2015).

In the molecular model, the electrode-electrolyte system atoms and molecules interact with each other via the well-known CHARMM empirical potential MacKerell Jr (2001). This often employed potential depends on atomic charges and proximity. It is composed of short-range van der Waals interactions, long-range Coulomb interactions, and intramolecular covalent bonds:

[TABLE]

where $\varepsilon_{ab}$ and $\sigma_{ab}$ are the usual Lennard-Jones (LJ) pair parameters for species $a$ and $b$, $q_{\alpha}$ is the charge of atom $\alpha$, $r_{\alpha\beta}=\|\mathbf{x}_{\alpha}-\mathbf{x}_{\beta}\|$ is the distance between atoms $\alpha$ and $\beta$, and $a,b,c,d$ index atom types. We employ a particle-particle particle-mesh (PPPM) solver Hockney and Eastwood (1988) to efficiently compute the long-range Coulomb interactions with a short/long cutoff of 12 Å . Here, the permittivity of free space is $\varepsilon_{0}$= 0.00552635 e/V-Å . The intra-molecular bonds are effected by harmonic potentials based on pair distances $r_{\alpha\beta}$, 3 atom angles $\theta_{\alpha\beta\gamma}$, and 4 atom dihedral angles $\phi_{\alpha\beta\gamma\mu}$, where $M_{I}$ is a set of like molecules.

III.1 Electrode

Given the form of the interatomic potential, Eq. (7), which includes Coulomb forces, we need the point charges $q_{\alpha}$ for the classical representation of the electrode via Eq. (7) and the response of the electrode point charges to electrical bias/external potential $V$. We use DFT with PBE/GGA level of theory Perdew, Burke, and Ernzerhof (1996) to obtain relaxed surface structures, and compute point charges via a Bader analysis of the charge density field Tang, Sanville, and Henkelman (2009).

Since only studies of small systems are feasible with DFT, we examined specific domains of a partially oxidized Ni electrode: bare Ni, partially oxidized Ni, and Ni covered by an NiO layer. For each, we create a small, laterally periodic system to calculate charge density using a FCC unit cell with a 3.52 Å lattice constant for Ni regions and a cubic B1 unit cell with a 4.17 Å lattice constant for NiO regions. We select (100)-oriented (non-polar) surfaces for both Ni and NiO. The surfaces neighbor vacuum regions, not representative electrolytes, for simplicity and under the assumption that the proximity of the electrolyte evokes only perturbations of the charge density. We employ an energy cutoff of 500 eV for the metallic systems and 800 eV for the systems containing oxides, together with an 8$\times$8$\times$1 $\Gamma$ centered $k$-point grid. After computing a baseline charge density, we add excess electrons to emulate negative charging of the electrode of interest. Finally, we obtain the electrode point charges $q^{\prime}_{\alpha}(\sigma)$:

[TABLE]

where $q_{\alpha}$ are the baseline point charges (such that $\sum_{\alpha}q_{\alpha}=0$), $\Delta q_{\alpha}$ are the perturbed point charges (corrected for the homogeneous background charge and normalized such that $\Delta q_{\alpha}$ = 1 $e$ for the surface atoms), $\sigma$ is the target surface charge density, and $\eta$ is surface atom density. The number of excess electrons added the systems to obtain the perturbed charge field is on the order of 0.1 $e$ per surface atom and on par with the perturbation needed to achieve the voltage bias $V$ in the range $V\approx$ 0–2.5 V in the experimental system. Lastly, the scaling of Eq. (8) with external voltage $V$ follows $V\approx\frac{\lambda_{D}}{\varepsilon_{0}}\sigma$; however, we recover the actual voltage in the MD systems via a Gauss box method described in Sec.​ III.

In addition to the point charges required for Coulomb interactions, short-range parameters $\varepsilon_{ab}$ and $\sigma_{ab}$ for the electrode interactions with the ionomer are needed. By assuming traditional Lorentz-Berthelot mixing, only the LJ self-self pair parameters are required. We obtain these from published, surface specific parameterizations: for Ni, $\sigma_{\text{NiNi}}=$ 2.274 Å and $\epsilon_{\text{NiNi}}$= 5.65 kcal/mol Ref. (52)(Table 1): and for NiO, $\sigma_{\text{OO}}=$ 1.292 Å and $\epsilon_{\text{OO}}$= 35.62 kcal/mol Ref. (53)(Table 1).111Ref. (53) reports parameters for a Buckingham potential which we convert to a LJ parameterization by matching the potential well location and stiffness. For simplicity (since no additional parameters are needed) and efficiency we neglect thermal motion of the electrode and fix the locations of its Ni and O atoms.

III.2 Electrolyte

The selected ionomer (PAP-DP-60) is an amorphous material composed of long polymer chains with N+ charge centers charge-balanced by OH- together with water. Each chain is comprised of charged (c) and neutral units (n), see Fig. 1, in a ratio c:n = 60:40 with lengths of 30 to 40 units in a random sequence. The mass density of the ionomer is 1.1 g/cm3 with about 40% water by weight at room temperature and pH in the range 13 to 14+.

Using these experimental measurements, we created representative models of the ionomer electrolyte, see Fig. 2. First we created chains from the units using a random sequence that respected the experimental bounds on length and c:n ratio. Next we added molecular water and OH- ions. To achieve the experimentally measured density, we compressed the ionomer-electrolyte mixture then relaxed the compressed configuration at 2000 K via isothermal dynamics to relieve unphysical local configurations for 0.4 ns. Finally, we cooled the system to 300 K at 10 K/ps and let it equilibrate for 0.1 ns. To complete the electrode-electrolyte system seen in Fig. 2, we added the primary Ni electrode plus a soft wall and counter electrode to bound the system. The separation of the capping wall and the counter electrode was expedient due to different densities and lengths of the replicas while maintaining the same effective gradient in external potential. Due to computational cost we can only place the counter electrode 400 nm from the surface of the primary electrode. As stated in the introduction, our focus is on the Ni electrode where the HER occurs and, hence, we model it in detail and simplify the counter electrode. For our purposes this is sufficient as long as the diffuse layer near the primary electrode relaxes to the bulk, and the region away from the primary electrode is not depleted of mobile species. As we will see in Sec.​ IV both these conditions are satisfied. Also, in preliminary studies, the lateral dimensions of the systems where increased to the point where the charge density profiles and related measures converged to within the expected statistical noise.

Since the ionomer is amorphous, we generated a number of statistically equivalent replica systems to improve sampling and reduce finite size effects. The isothermal dynamics of the electrolyte were evolved with a Nosé-Hoover algorithm which accommodated the rigidity of the intramolecular bonds. All reported results are the average of 8 configurations time-averaged and simulated for 0.4 ns each with a 1 fs time-step.

Finally, we use a so-called Gauss box to recover the electric potential $\phi$ from the point charges $q_{\alpha}$. Briefly, starting with Gauss’s law for a quasi-one dimensional system, with electric field $\mathbf{E}\equiv-\boldsymbol{\nabla}\phi=E\mathbf{e}_{3}$ and cross-sectional area $A$, we obtain:

[TABLE]

where $Q(\tilde{z})$ is the total/net charge in $z\in[0,\tilde{z}]$ and $\left<\bullet\right>$ is the average over replica systems and steady-state trajectories in each system. In a similar fashion, we estimate charge density, dipole density, and species concentration profiles using coarse-grained point charges $q_{\alpha}$, atomic dipoles $q_{\alpha}\mathbf{x}_{\alpha}$, and species type. For instance, the water dipole density is given by:

[TABLE]

where $z_{I}$ is the location of points on the coarse-grained sampling grid, $V_{I}$ is the effective volume associated with $I$, $w$ is a partition-of-unity, tent-like kernel function, and

[TABLE]

In Eq. (10), the second term is a correction for net charge in a particular bin $I$. Also, we use a kernel width of 3 Å based on the size of water molecules, and cutoff the averaging kernel $w(z)$ at electrode surface so that the average only includes the time-dependent data in the electrode. The water dipole density is used to determine water orientation and the species concentrations indicate which species are present in the compact and diffuse layers. Also, we extract coarse-grained profiles of radial distribution functions (RDFs) of particular species to quantify changes in solvation with 3 Å resolution. This information allows us to estimate the location of the OHP and its response to bias.

IV Results

We compare the response of the electrolyte to electrode bias voltage using a variety of configurational metrics in order to provide insight and support to the experimental hypotheses described in the Introduction. We focus on (a) bare Ni, (b) fully NiO covered and (c) partially NiO covered regions of the Ni electrode. Fig. 3 and Fig. 4 show the qualitative differences in surface coverage of the relevant species and water orientation for the three cases. Fig. 3 and Fig. 4 will be discussed in more detail in the following sections.

IV.1 Electrode charge distributions

Using first principles methods described in Sec.​ III.1 we obtain the relaxed configuration and electron density of the electrode with and without an oxide layer. Fig. 5 shows a top down view of the charge density of a Ni surface partially covered by a NiO monolayer. The NiO island is in the lower left of Fig. 5 and the exposed portion of the underlying Ni surface can be seen in the remainder of the figure. The initially square island relaxes to a more diamond-like shape but the variation in charge density across the partial layer is minimal, with charge transfer from the O to Ni nuclei as expected. Also, the charge density of the exposed Ni metal is relatively unperturbed by the oxide layer. The mean separation between the Ni and the NiO layer is $\approx$2 Å. With a sequence of related simulations we examined where the excess electrons (induced by electrode charging) reside. Fig. 6 plots the point charges (in excess of normal valence: metal Ni 0, oxide Ni +2, O -2) for bare Ni and Ni covered by 1, 2, and 3 monolayers of NiO (the point charges for an overall neutral system and one with excess electrons are extracted with a Bader method described in Sec.​ III and then differenced). As expected, in the bare Ni, the excess electron density resides primarily in the outermost layer of Ni. For a single layer of NiO covering Ni $\approx$ 80% of the excess electrons reside near the O, $\approx$ 10% near the Ni in the oxide layer and $\approx$ 10% near the Ni in the underlying metal. This charge splitting induces a surface dipole moment. It is apparent that the charge distribution and the dipole effects become more complex as more oxide layers are added to the system.

IV.2 Bare metal region of the electrode

With classical MD, we simulated the electrode-electrolyte systems over a sequence of surface charges: $\Delta q=\{0.0,0.05,0.10,0.15\}$ $e$ per surface atom, corresponding to $\sigma=\{0.007,0.014,0.021,0.028\}$ $e/\AA^{2}$. These surface charges were related to voltage through the Gauss box technique described in Sec.​ III.

Fig. 7 shows that the voltages at the electrode surface ($z$=0) range from 0 to $\approx$ 2 V for these surface charges. Also, it is apparent from the charge density profiles that: (a) both a double layer next to the electrode and a diffuse layer grading to the bulk form exist, and (b) the decay to the bulk region with zero net charge occurs within 20Å from the outer layer of the electrode. We fit the potential profile in the diffuse region (approximately 5–20Å from the electrode surface) to the Poisson-Boltzmann solution in Eq. (5) and obtain $\ell\approx$ 0.5 Å for $q=$ -0.1 $e$/atom. This allowed us to estimate the effective dielectric in the diffuse layer $\epsilon=\frac{e\sigma}{2k_{B}T}\ell\approx 24.4\epsilon_{0}$ which is considerably less than pure water ($\epsilon\approx 80\epsilon_{0}$).

In Fig. 7 we also see that the region of significant dipole moment near the electrode decays to the bulk on the same scale as the potential that induces the orientation of the water molecules. Most significantly, the dipole moment is mostly negative and increases in magnitude with applied bias, which corresponds to the expected conformation of the H atoms of the water molecules oriented closest to the negative electrode, and this orientation becoming more dominant with increased bias. This finding is corroborated by direct observation, as in Fig. 3 and Fig. 4 (at $q=$ -0.1 $e$/atom), where the white, partially positively charged H atoms lie on the electrode surface with the red, partially negatively charged O atoms of the same water molecules pointing away from the surface. Either one or both hydrogens is directly associated with distinct metal surface atoms. In a second layer, we observe some water molecules bridging water on the surface with their hydrogens associated with the oxygens of the surface adsorbed waters. Even at zero bias it appears there is some orientational preference for the water molecules which must be induced by the particular LJ parameters for H and O atoms in the water molecules. (The preferential alignment of water with uncharged surfaces is widely observed Carrasco, Hodgson, and Michaelides (2012).)

Fig. 8 shows that at zero bias, H2O, OH- and the N+ of the ionomer are adsorbed to the electrode surface and maintain uniform bulk concentrations away from the electrode. At higher biases OH- is excluded near the electrode, and water fills the vacated region and forms structured layers. The spatial distribution of N+ (averaged over time and the replica systems) is effectively uniform in the electrolyte region except there is a significantly higher concentration of N+ at the interface (with corresponding neighboring region of depletion). This appears to be due to the fact that ends of the chains are relatively mobile, unentangled and attracted to surface of the negative electrode. The concentration of N+ at the electrode more than doubles with the applied bias of $\approx$ 2 V and the response appears to be nonlinear, most likely due to an energy barrier in straightening the ends of the ionomer chains. This can be seen as competition of elastic versus Coulomb forces for the relatively immobile N+. Note that these biases are high enough to dissociate water in the real system and the ionomer near the electrode might be likewise affected. The mobile species, H2O and OH-, respond as mentioned. From Fig. 3, it apparent that no ring-like water structures form and that, generally, the N+ of the ionomer appear in regions of low concentrations of surface adsorbed water. We presume that the near surface N+ of the ionomer displace water and disrupt the tendency for the water to form regular structures. The OH- transitions from significant surface adsorbed species at no bias to progressively vacated regions with increased voltage. Here again, a threshold (determined by the relative magnitude of the overall electric field separating OH- and N+ versus the LJ and Coulomb forces binding them) appears to be operating. For biases less than $q=$-0.05 $e$/atom (not shown) OH- remains bound to N+. Note that effectively no OH- is present near the electrode at the $q=$-0.10 $e$/atom surface charge as shown in Fig. 5a. Water, on the other hand, displays a surface adsorbed layer, a depletion zone, and then a uniform bulk concentration at no bias. With increasing voltage more structured water layers appear near the surface and the closest layer becomes more packed (higher surface coverage).

Fig. 9 shows the spatially binned RDFs for zero bias and $q=$-0.10 $e$/atom. First of all, we see that OH- is closely coordinated with H2O and this strong association does not change wherever OH- is present. At voltage the only effect is decreasing the magnitude (not the location) of the RDF peak due to decreased OH- concentration near the electrode surface (refer to Fig. 8). Note no OH- is present at 0–6Å from the electrode at $q=$-0.10 $e$/atom. The coordination of N+ with H2O is uniform at zero bias and relatively unaffected at high bias which implies the OHP defined by H2O solvation of N+ is within 3 Å of the outer layer of electrode surface. The effects of distance from the electrode and bias on the coordination of N+ with OH- are more complex due to changing local electric field and concentration. Given the location of the RDF peaks, the primary coordination shell of OH- appears to mix with that of H2O. Most significant of the changes is that N+ coordination with OH- tightens and sharpens at high bias and near the electrode, particularly in the 6–9Å bin at $q$= -0.1 $e$/atom. In the next bin, 9–12Å, we see two peaks corresponding to a two shell arrangement; and, by the 15–18Å bin, the RDF resembles that of the bulk in both the zero and high bias cases. If we use the coordination of N+ with OH- to indicate the location of the OHP, we see the OHP shift away from the electrode surface reaching $\approx$ 20 Å at the highest bias, as Fig. 8 indicates, due to the dissociation of the OH- from the N+ at high bias.

IV.3 Oxide covered region of the electrode

As observed in Fig. 3 and Fig. 4, the presence of a (single layer of) NiO qualitatively changes the response of the HER electrode to electrical bias. Comparing the surface coverage of a bare Ni region of the electrode with partially and fully covered NiO regions, we see approximately 40% increase of adsorbed water where there is oxide on the surface. Fig. 3 also illustrates that there is distinctly increased water surface concentration at the NiO boundary. The side view of the electrode in Fig. 3 and the closeup in Fig. 4 shows water molecules on the surface of the NiO covered electrode with one of the O-H legs of the water molecule flat on the surface and the other with H embedded in the oxide layer in direct association with an oxide O atom. These effects are, in large part, due to the differences in the unbiased point charges and induced charge distribution (refer to Fig. 6). Since oxides are dielectrics, they change local and long-range electric fields by screening due to induced and/or permanent dipoles.

At the same surface charges and comparable voltages to the simulation of the bare Ni region of the electrode discussed in the previous section, we compute the electrolyte’s response to electrical bias for a Ni electrode with a monolayer of NiO as shown in Fig. 10. There are some similarities and noticeable differences in the profiles compared to the bare Ni case shown in Fig. 7. Firstly, there is a larger counter charge density next to the electrode in the NiO-on-Ni case, which we can associate with the partially positively charged H in the increased density of surface water molecules. Examining the dipole density, the H-toward orientation is strong and barely shifts with increasing bias due to the strong, local polarization of the electrode itself. This near-surface dipole density is reflected in the electric potential which becomes constant but not zero in the far-field. This also leads to a lower effective voltage difference for the same surface charge, as expected with a dielectric layer. In addition, the more prominent positive dipole peak seen in Fig. 10 compared to Fig. 7 is apparently due the dipole created by neighboring waters to since no inversions of water molecules with the H toward the surface were observed.

Likewise the concentration profiles shown in Fig. 11 are comparable to Fig. 8 but also display distinct differences. As in Fig. 8, the water density is uniform in the far-field and shows distinct signs of structured layers near the electrode. On the other-hand, the (time-averaged) water density at the electrode surface is much higher than in the bare metal case (approximately 40%, as corroborated by the snapshots shown in Fig. 3). Also, this spike in concentration is essentially unaffected by bias, a fact corroborated by nearly constant dipole density at the surface. Remarkably, significant OH- remains associated with the electrode surface at high bias, but concentration decreases with increased bias (similar behavior has been observed in other systems Lee, Mani, and Templeton (2015)). Otherwise, the OH- is similar to that for bare Ni, with a depletion zone increasing with increased bias and a nearly constant bulk concentration. The N+ concentration is also similar to that for the bare Ni, shown in Fig. 8, but has a pronounced double peak near the interface. This feature is perhaps related to significant residual OH- concentration at the surface with bias.

The spatially binned RDFs shown in Fig. 12 are qualitatively similar to Fig. 9, implying that the electrode surface is not radically changing the coordination structures. There is a slight variation in OH- and H2O concentration at low bias which affects OH-:H2O and N+:H2O coordination that is not observed in Fig. 9. The presence of OH- near the N+ closest to the electrode in this case also appears in the N+:OH- RDF but it is omitted in Fig. 9 due to sampling noise and an effort to maintain clarity. Also, in this case the N+:OH- RDF clearly relaxes to bulk distribution at zero bias.

V Discussion

In response to the hypotheses posed in Bates et al. Bates et al. (2015) and reiterated in the Introduction, we find:

• •

The electrical potential is altered by the ionomer and the oxide layer primarily due to dielectric effects of the immobile charges and their interplay with the mobile ones, as observed in Fig. 7 and Fig. 10. A significant concentration of positive charge from the ionomer lies near the surface of the electrode, as can be inferred from Fig. 8 and Fig. 11.

• •

The ionomer N+ are coordinated with both OH- and H2O. Fig. 9 and Fig. 12 indicate that solvation and hence the location of the OHP is changing with electrode bias voltage and depending on the solvating species ranges from 3 Å off the electrode surface to 15 Å. A region depleted of OH- forms near the electrode with electrical bias and grows with increased bias. Since the N+ are relatively immobile, there are significant N+ atoms in the depletion region, but these N+ remain coordinated with water.

• •

On the surface of the electrode, refer to Fig. 3, the N+ in the ionomer apparently displaces water and disrupts the tendency of the water to form ordered surface structures like rings, as in Ref. (32). From Fig. 8 and Fig. 11 is it clear that significant concentration of the ionomer N+ lie near the electrode surface.

• •

The water molecules near the electrode surface orient in the manner dictated by the electric field (with O further away from the surface) in a mixture of either one or both hydrogens associated with surface atoms; however, in the surfaces with an oxide layer the primary association is with a surface oxygen and the configuration where one of the H legs of the water molecule lies on the electrode surface seems to dominate, as can be seen in Fig. 3 and Fig. 4.

• •

Approximately 40% of the bare metal surface is coordinated with water molecules and this surface coverage almost doubles in areas covered by an oxide layer.

• •

Nanoscale heterogeneity such as at the edge of a oxide region provides local configurations favorable for water adsorption, as evidenced by the regions of high water coverage at the edge of the partially oxidized electrode in Fig. 3 and Fig. 4. This qualitative observation may corroborate how rough versus smooth electrodes exhibit qualitatively different reactivity, see Ref. (27) (Sec. 3.1 and 4.2) and references therein.

• •

The presence of a NiO monolayer apparently allows OH- to remain near the electrode at high negative surface charges, compare Fig. 8 and Fig. 11.

Although we are presently unable to simulate the relevant reactions at the necessary time- (molecular diffusion) and length- (diffuse layer) scales, the fact that the reactions are fast compared with the processes MD represents well gives us some confidence that the results are representative of steady-state conditions. Also, the H2 produced at the HER electrode is small and mobile, and not likely to affect the components of the electrolyte. If positive ions are present e.g. K+ from dissolved KOH, they will likely absorb to the electrode surface and hence screen the remainder of the electrolyte, having stronger effects on performance.

VI Conclusion

We were able to simulate and examine the HER electrode-ionomer electrolyte interface with atomic detail using a combination of DFT and classical MD techniques. We observed configuration changes in response to external bias and the oxide coverage of the electrode. Information of this nature is relevant to the efficiency of the water splitting process. In particular the concentration of the reactants, water, and the electric field at the interface have a direct relation to H2 production.

In future work, we will explore in more detail the effect of surface roughness/nano-structuring and the presence of other phases, such as Cr2O3, on the HER electrode. Resolution of the dielectric field Mandadapu, Templeton, and Lee (2013) may also shed light on the operation of these types of water splitting cells and give rise to more accurate and informative theories. Also, the entropic changes due to changes in pH discussed in Ref. (56) are likely correlated with structural characterizations presented in this work, such as the spatial variation of radial distributions, and are another topic for potential future work.

Acknowledgments

We would like to thank Norman Bartelt and Jeremy Templeton (Sandia) for insightful discussions on this work. This work enabled by VASP (TU Vienna, https://www.vasp.at), LAMMPS (Sandia, https://lammps.sandia.gov), and Bader (Univ. Texas http://theory.cm.utexas.edu/henkelman/code/bader/). This work was supported by the U.S. Department of Energy, Office of Energy Efficiency and Renewable Energy (EERE), specifically the Fuel Cell Technologies Office Sandia National Laboratories is a multimission laboratory managed and operated by National Technology and Engineering Solutions of Sandia, LLC., a wholly owned subsidiary of Honeywell International, Inc., for the U.S. Department of Energy’s National Nuclear Security Administration under contract DE-NA0003525. The views expressed in the article do not necessarily represent the views of the U.S. Department of Energy or the United States Government.