Insight into hydrogen production through molecular simulation of an electrode-ionomer electrolyte system
Reese E. Jones William C. Tucker Matthew J. L. Mills Sanjeev Mukerjee


rjones@sandia.gov 
Sandia National Laboratories, Livermore, CA 94551, USA

Sandia National Laboratories, Livermore, CA 94551, USA

Sandia National Laboratories, Livermore, CA 94551, USA

Northeastern University, Boston, MA 02115, USA









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 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.


Insight into hydrogen production through molecular simulation of an electrode-ionomer electrolyte system
    S. Mukerjee
    December 30, 2023
========================================================================================================






§ 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 <cit.>.
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 <cit.>.
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 <cit.> 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 <cit.> 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 NiO_x has an affinity for adsorbed OH^-, in-line with Markovic's theory of enhanced hydrogen evolution reaction on composite metal/metal-oxide surfaces <cit.>.


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 calculations such as filhol2006elucidation,jinnouchi2008electronic,skulason2010modeling,skulason2007density,karlberg2007cyclic,taylor2006first,rossmeisl2006calculated,janik2009first,otani2008electrode, 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 <cit.>.
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 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 (NiO_x) 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 thiel1987interaction,henderson2002interaction,hodgson2009water,carrasco2012molecular,stuve2012ionization.
The examination of the water-surface interactions has been pursued in great detail by combinations of microscopy and simulation (typically in vacuum), as in <ref>, 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:theory 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 <cit.>.
In sec:method 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:results and are discussed in sec:discussion in light of the hypotheses of Bates <cit.>



§ 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 <cit.>, Gouy  <cit.>, Chapman <cit.>, Stern <cit.> and others, and is still a topic of current research <cit.>.
Here we will briefly review the relevant theory to assist in interpreting the molecular results in sec:results.

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:

    -· (ϵϕ ) = ρ .

Here, ϕ is the electric potential and ϵ is the electric permittivity.
The net charge density ρ is given by:

    ρ = e ∑_a z_a c_a  ,

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, e z_a  ϕ + 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 ϕ:

    c_a = c̅_a exp(e z_a ϕ/k_B T)  ,

where c̅_a is the far-field/bulk concentration of species a.
Substituting eq:rho, and eq:Boltzmann in eq:Poisson results in the Poisson-Boltzmann (PB) equation.
The characteristic thickness of the diffuse layer is given by the Debye length 

    λ_D = √(k_B T ϵ/e^2 ∑_a z_a^2 c̅_a) ,

which is the similarity parameter in the solution of the linearized PB equations.
For more details see <ref>(Sec.13.3), <ref> (Chap.12), and <ref> (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 <cit.>.
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 <cit.> is:

    ϕ(x) = -2 k_B T/eln(1 + x/ℓ) + ϕ_0

and 

    c(x) = c_0/(1 + x/ℓ)^2 ,

where ϕ_0 is the potential at x=0, ℓ = 2 ϵ k_B T/e σ and σ is the surface charge of the electrode.
Note that the dielectric of the electrolyte, ϵ, is assumed to be spatially uniform in this derivation.
This diffuse layer solution has an electric field = - ϕ that is zero far away from the electrode and has the magnitude σ/ϵ at x=0, which is consistent with Gauss's law applied to the electrode (this discussed in more detail in the next section).




§ 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 <cit.>.

In the molecular model, the electrode-electrolyte system atoms and molecules interact with each other via the well-known CHARMM empirical potential <cit.>.
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:

    Φ(_α) 
       =  ∑_α <β 4ε_ab( (σ_ab/r_αβ)^12 - (σ_ab/r_αβ)^6  )_van der Waals
    + 1/4 πϵ_0∑_α <βq_α q_β/r_αβ_Coulomb
       + ∑_I [ ∑_α,β∈ M_I  k_ab r^2_αβ
    + ∑_α,β,γ∈ M_I  k_abcθ^2_αβγ
    + ∑_α,β,γ,μ∈ M_I  k_abcdϕ^2_αβγμ]_covalent

where ε_ab and σ_ab are the usual Lennard-Jones (LJ) pair parameters for species a and b,  q_α is the charge of atom α, r_αβ = _α - _β is the distance between atoms α and β, and a,b,c,d index atom types.
We employ a particle-particle particle-mesh (PPPM) solver <cit.> to efficiently compute the long-range Coulomb interactions with a short/long cutoff of 12 Å .
Here, the permittivity of free space is ε_0= 0.00552635 e/V-Å .
The intra-molecular bonds are effected by harmonic potentials based on pair distances r_αβ, 3 atom angles θ_αβγ, and 4 atom dihedral angles ϕ_αβγμ, where M_I is a set of like molecules.



 §.§ Electrode
 
Given the form of the interatomic potential, eq:md_potential, which includes Coulomb forces,  we need the point charges q_α for the classical representation of the electrode via eq:md_potential and the response of the electrode point charges to electrical bias/external potential V.
We use DFT with PBE/GGA level of theory <cit.> to obtain relaxed surface structures, and compute point charges via a Bader analysis of the charge density field <cit.>. 

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×8×1 Γ 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'_α(σ):

    q'_α(σ) = q_α - σ/e η Δ q_α

where q_α are the baseline point charges (such that ∑_α q_α = 0), Δ q_α are the perturbed point charges (corrected for the homogeneous background charge and normalized such that Δ q_α = 1 e for the surface atoms), σ is the target surface charge density, and η 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 ≈ 0–2.5 V in the experimental system.
Lastly, the scaling of eq:qperturbed with external voltage V follows V≈λ_D/ε_0σ; however, we recover the actual voltage in the MD systems via a Gauss box method described in sec:method.

In addition to the point charges required for Coulomb interactions, short-range parameters ε_ab and σ_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, 
σ_NiNi = 2.274 Å and ϵ_NiNi= 5.65 kcal/mol
<ref>(Table 1):
and for NiO,  
σ_OO   = 1.292 Å and ϵ_OO=  35.62 kcal/mol 
<ref>(Table 1).
[<ref> 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.



 §.§ 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:units, 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/cm^3 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:system.
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:system, 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:results 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 ϕ from the point charges q_α. 
Briefly, starting with Gauss's law for a quasi-one dimensional system, with electric field ≡ -ϕ = E _3 and cross-sectional area A, we obtain:

    ϕ(z)    = 
    ϕ(0)
    -1/ε_0 A∫_0^z < ∑_α | z_α∈[0,z̃] q_α> ^Q(z̃) dz̃ ,

where Q(z̃) is the total/net charge in z∈[0,z̃] and < ∙> 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_α, atomic dipoles q_α_α, and species type.
For instance, the water dipole density is given by:

    μ(z_I) = 1/V_I< ∑_α  q_α z_α w(z_I-z_α) > - 1/V_I q_I z_I  ,

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 

    q_I =  < ∑_α  q_α w(z_I-z_α) >  .

In  eq:dipole, 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.



§ 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:surface_species and fig:surface_orientation show the qualitative differences in surface coverage of the relevant species and water orientation for the three cases.
fig:surface_species and fig:surface_orientation will be discussed in more detail in the following sections.







 §.§ Electrode charge distributions

Using first principles methods described in sec:electrode we obtain the relaxed configuration and electron density of the electrode with and without an oxide layer.
fig:electrode_charge_density 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:electrode_charge_density 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 ≈2 Å.  
With a sequence of related simulations we examined where the excess electrons (induced by electrode charging) reside.
fig:electrode_point_charges 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:method 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 ≈ 80% of the excess electrons reside near the O, ≈ 10% near the Ni in the oxide layer and ≈ 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.







 §.§ Bare metal region of the electrode

With classical MD, we simulated the electrode-electrolyte systems over a sequence of surface charges: Δ q = { 0.0, 0.05, 0.10, 0.15} e per surface atom, corresponding to σ = { 0.007, 0.014, 0.021, 0.028} e/Å^2.
These surface charges were related to voltage through the Gauss box technique described in sec:method.

fig:potential-Ni shows that the voltages at the electrode surface (z=0) range from 0 to ≈ 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:PB and obtain ℓ≈ 0.5 Å for q = -0.1 e/atom.
This allowed us to estimate the effective dielectric in the diffuse layer ϵ = e σ/2 k_B Tℓ≈ 24.4 ϵ_0 which is considerably less than pure water (ϵ≈ 80 ϵ_0).


In fig:potential-Ni 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:surface_species and fig:surface_orientation (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 <cit.>.)

fig:profiles-Ni shows that at zero bias, H_2O, 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 ≈ 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, H_2O and OH^-, respond as mentioned.
From fig:surface_species, 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:electrode_charge_densitya.
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:rdf_profiles-Ni 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 H_2O 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:profiles-Ni).
Note no OH^- is present at 0–6Å from the electrode at q=-0.10 e/atom.
The coordination of N^+ with H_2O is uniform at zero bias and relatively unaffected at high bias which implies the OHP defined by H_2O 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 H_2O.
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 ≈ 20 Å at the highest bias, as fig:profiles-Ni indicates, due to the dissociation of the OH^- from the N^+ at high bias.











 §.§ Oxide covered region of the electrode

As observed in fig:surface_species and fig:surface_orientation, 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:surface_species also illustrates that there is distinctly increased water surface concentration at the NiO boundary.
The side view of the electrode in fig:surface_species and the closeup in fig:surface_orientation 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:electrode_point_charges).
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:potential-NiNiO.
There are some similarities and noticeable differences in the profiles compared to the bare Ni case shown in fig:potential-Ni.
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:potential-NiNiO compared to fig:potential-Ni 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:profiles-NiNiO are comparable to fig:profiles-Ni but also display distinct differences.
As in fig:profiles-Ni, 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:surface_species).
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 <cit.>).
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:profiles-Ni, 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:rdf_profiles-NiNiO are qualitatively similar to fig:rdf_profiles-Ni, implying that the electrode surface is not radically changing the coordination structures.
There is a slight variation in OH^- and H_2O concentration at low bias which affects OH^-:H_2O and N^+:H_2O coordination that is not observed in fig:rdf_profiles-Ni.
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:rdf_profiles-Ni 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.











§ DISCUSSION
 

In response to the hypotheses posed in Bates <cit.> 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:potential-Ni and fig:potential-NiNiO. 
A significant concentration of positive charge from the ionomer lies near the surface of the electrode, as can be inferred from fig:profiles-Ni and fig:profiles-NiNiO.

  * The ionomer N^+ are coordinated with both OH^- and H_2O.
fig:rdf_profiles-Ni and fig:rdf_profiles-NiNiO 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:surface_species, 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>.
From fig:profiles-Ni and fig:profiles-NiNiO 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:surface_species and fig:surface_orientation.

  * 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:surface_species and fig:surface_orientation.
This qualitative observation may corroborate how rough versus smooth electrodes exhibit qualitatively different reactivity, see <ref> (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:profiles-Ni and fig:profiles-NiNiO.

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 H_2 produced at the HER electrode is  small and mobile, and not likely to affect the components of the electrolyte.
If positive ions are present 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.




§ 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 H_2 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 Cr_2O_3, on the HER electrode.
Resolution of the dielectric field <cit.> 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> 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, ), 
LAMMPS (Sandia, ), and 
Bader (Univ. Texas ).
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.


56
fxundefined [1]
 ifx#1

fnum [1]
 #1firstoftwo
 secondoftwo
 
fx [1]
 #1firstoftwo
 secondoftwo
 





noop [0]secondoftwo

ref[1]@startlink#1@href
href[1]#1@endlink
anitize@url [0]`
12`$12`&12`#12`1̂2`_12`%12
startlink[1]
endlink[0]

rl [1]href #1












@bib@innerbibempty

[Rahman and Andrews(2006)]rahman2006special
  author author S. Rahman and author C. J. Andrews, @noop  journal journal
  Proceedings of the IEEE volume 94, pages 1781 (year 2006)NoStop
[Council(2004)]nrc2004hydrogen
  author author N. R. Council, @noop  title The Hydrogen Economy:
  Opportunities, Costs, Barriers, and R&D Needs (publisher
  National Academies Press, address
  https://www.nap.edu/read/10922/chapter/10, year
  2004)NoStop
[Navlani-García et al.(2018)Navlani-García, Mori,
  Kuwahara, and Yamashita]navlani2018recent
  author author M. Navlani-García, author K. Mori, author Y. Kuwahara,   and author H. Yamashita,  @noop  journal journal NPG Asia
  Materials , pages 1 (year 2018)NoStop
[Bates et al.(2015)Bates,
  Jia, Ramaswamy, Allen, and  Mukerjee]bates2015composite
  author author M. K. Bates, author Q. Jia, author N. Ramaswamy, author
  R. J. Allen,  and author
  S. Mukerjee, @noop  journal journal The Journal of Physical Chemistry C  volume 119, pages 5467 (year
  2015)NoStop
[Ünlü et al.(2011)Ünlü, Abbott, Ramaswamy,
  Ren, Mukerjee, and Kohl]unlu2011analysis
  author author M. Ünlü, author D. Abbott, author N. Ramaswamy,
  author X. Ren, author
  S. Mukerjee,  and author
  P. A. Kohl, @noop  journal journal Journal of The Electrochemical Society  volume 158, pages B1423 (year 2011)NoStop
[Subbaraman et al.(2011)Subbaraman, Tripkovic, Strmcnik,
  Chang, Uchimura, Paulikas,
  Stamenkovic, and Markovic]subbaraman2011enhancing
  author author R. Subbaraman, author D. Tripkovic, author D. Strmcnik,
  author K.-C. Chang, author M. Uchimura, author
  A. P. Paulikas, author
  V. Stamenkovic,  and author
  N. M. Markovic, @noop  journal journal Science volume
  334, pages 1256 (year 2011)NoStop
[Strmcnik et al.(2013)Strmcnik, Uchimura, Wang, Subbaraman, Danilovic, Van Der Vliet,
  Paulikas, Stamenkovic, and Markovic]strmcnik2013improving
  author author D. Strmcnik, author M. Uchimura,
  author C. Wang, author
  R. Subbaraman, author
  N. Danilovic, author
  D. Van Der Vliet, author
  A. P. Paulikas, author
  V. R. Stamenkovic,  and author
  N. M. Markovic, @noop  journal journal Nature chemistry volume 5, pages 300 (year 2013)NoStop
[Danilovic et al.(2012)Danilovic, Subbaraman, Strmcnik,
  Chang, Paulikas, Stamenkovic, and Markovic]danilovic2012enhancing
  author author N. Danilovic, author R. Subbaraman, author D. Strmcnik, author K.-C. Chang,
  author A. Paulikas, author V. Stamenkovic,  and author N. M. Markovic, @noop 
  journal journal Angewandte Chemie volume 124, pages 12663 (year
  2012)NoStop
[Filhol and Neurock(2006)]filhol2006elucidation
  author author J.-S. Filhol and author M. Neurock, @noop  journal journal
  Angewandte Chemie International Edition volume
  45, pages 402 (year 2006)NoStop
[Jinnouchi and Anderson(2008)]jinnouchi2008electronic
  author author R. Jinnouchi and author A. B. Anderson, @noop  journal journal
  Physical Review B volume 77, pages
  245417 (year 2008)NoStop
[Skúlason et al.(2010)Skúlason, Tripkovic, Björketun,
  Gudmundsdottir, Karlberg, Rossmeisl, Bligaard, Jónsson, and  Nørskov]skulason2010modeling
  author author E. Skúlason, author V. Tripkovic, author M. E. Björketun, author S. Gudmundsdottir, author G. Karlberg, author J. Rossmeisl,
  author T. Bligaard, author H. Jónsson,  and author J. K. Nørskov, @noop 
  journal journal The Journal of Physical Chemistry
  C volume 114, pages 18182 (year 2010)NoStop
[Skúlason et al.(2007)Skúlason, Karlberg, Rossmeisl,
  Bligaard, Greeley, Jónsson, and Nørskov]skulason2007density
  author author E. Skúlason, author G. S. Karlberg, author J. Rossmeisl,
  author T. Bligaard, author J. Greeley, author
  H. Jónsson,  and author
  J. K. Nørskov, @noop 
  journal journal Physical Chemistry Chemical
  Physics volume 9, pages 3241
  (year 2007)NoStop
[Karlberg et al.(2007)Karlberg, Jaramillo, Skulason,
  Rossmeisl, Bligaard, and Nørskov]karlberg2007cyclic
  author author G. Karlberg, author T. Jaramillo,
  author E. Skulason, author J. Rossmeisl, author
  T. Bligaard,  and author
  J. K. Nørskov, @noop 
  journal journal Physical review letters  volume 99, pages 126101 (year 2007)NoStop
[Taylor et al.(2006)Taylor,
  Wasileski, Filhol, and Neurock]taylor2006first
  author author C. D. Taylor, author S. A. Wasileski, author J.-S. Filhol,  and author M. Neurock, @noop  journal journal
  Physical Review B volume 73, pages
  165402 (year 2006)NoStop
[Rossmeisl et al.(2006)Rossmeisl, Nørskov, Taylor,
  Janik, and Neurock]rossmeisl2006calculated
  author author J. Rossmeisl, author J. K. Nørskov, author C. D. Taylor, author M. J. Janik,   and author M. Neurock,  @noop  journal journal The Journal of
  Physical Chemistry B volume 110, pages
  21833 (year 2006)NoStop
[Janik, Taylor, and Neurock(2009)]janik2009first
  author author M. J. Janik, author C. D. Taylor,   and author M. Neurock,  @noop  journal journal Journal of the
  Electrochemical Society volume 156, pages B126 (year 2009)NoStop
[Otani et al.(2008)Otani,
  Hamada, Sugino, Morikawa,
  Okamoto, and Ikeshoji]otani2008electrode
  author author M. Otani, author I. Hamada,
  author O. Sugino, author Y. Morikawa, author
  Y. Okamoto,  and author
  T. Ikeshoji, @noop  journal journal Journal of the Physical Society of Japan  volume 77, pages 024802 (year 2008)NoStop
[Spohr(1999)]spohr1999molecular
  author author E. Spohr, @noop  journal journal
  Electrochimica Acta volume 44, pages
  1697 (year 1999)NoStop
[Willard et al.(2009)Willard, Reed, Madden, and Chandler]willard2009water
  author author A. P. Willard, author S. K. Reed,
  author P. A. Madden,  and  author D. Chandler, @noop
   journal journal Faraday Discussions volume 141, pages 423 (year
  2009)NoStop
[Limmer et al.(2013)Limmer,
  Willard, Madden, and Chandler]limmer2013hydration
  author author D. T. Limmer, author A. P. Willard,
  author P. Madden,  and author D. Chandler, @noop 
  journal journal Proceedings of the National
  Academy of Sciences volume 110, pages
  4200 (year 2013)NoStop
[Yeh, Janik, and Maranas(2013)]yeh2013molecular
  author author K.-Y. Yeh, author M. J. Janik,   and author J. K. Maranas,  @noop  journal journal Electrochimica
  Acta volume 101, pages 308 (year 2013)NoStop
[Hughes and Walsh(2014)]hughes2014structure
  author author Z. E. Hughes and author T. R. Walsh, @noop  journal journal Journal
  of colloid and interface science volume 436,  pages 99 (year 2014)NoStop
[Dewan et al.(2014)Dewan,
  Carnevale, Bankura, Eftekhari-Bafrooei, Fiorin, Klein, and  Borguet]dewan2014structure
  author author S. Dewan, author V. Carnevale,
  author A. Bankura, author A. Eftekhari-Bafrooei, author G. Fiorin, author
  M. L. Klein,  and author
  E. Borguet, @noop  journal journal Langmuir volume
  30, pages 8056 (year 2014)NoStop
[Limmer et al.(2015)Limmer,
  Willard, Madden, and Chandler]limmer2015water
  author author D. T. Limmer, author A. P. Willard,
  author P. A. Madden,  and  author D. Chandler, @noop
   journal journal The Journal of Physical
  Chemistry C volume 119, pages 24016
  (year 2015)NoStop
[Takae and Onuki(2015)]takae2015fluctuations
  author author K. Takae and author A. Onuki,  @noop  journal journal The Journal of
  Chemical Physics volume 143, pages
  154503 (year 2015)NoStop
[Foroutan, Darvishi, and  Fatemi(2017)]foroutan2017structural
  author author M. Foroutan, author M. Darvishi,
   and author S. M. Fatemi,  @noop  journal journal Physical Review
  E volume 96, pages 033312 (year 2017)NoStop
[Thiel and Madey(1987)]thiel1987interaction
  author author P. A. Thiel and author T. E. Madey, @noop  journal journal Surface
  Science Reports volume 7, pages 211
  (year 1987)NoStop
[Henderson(2002)]henderson2002interaction
  author author M. A. Henderson, @noop  journal journal
  Surface Science Reports volume 46, pages 1 (year 2002)NoStop
[Hodgson and Haq(2009)]hodgson2009water
  author author A. Hodgson and author S. Haq,  @noop  journal journal Surface Science
  Reports volume 64, pages 381
  (year 2009)NoStop
[Carrasco, Hodgson, and  Michaelides(2012)]carrasco2012molecular
  author author J. Carrasco, author A. Hodgson,   and author A. Michaelides,  @noop  journal journal Nature materials  volume 11, pages 667 (year
  2012)NoStop
[Stuve(2012)]stuve2012ionization
  author author E. M. Stuve, @noop  title Ionization of water in
  interfacial electric fields: An electrochemical view,  (year
  2012)NoStop
[Nie et al.(2010)Nie,
  Feibelman, Bartelt, and Thürmer]nie2010pentagons
  author author S. Nie, author P. J. Feibelman,
  author N. Bartelt,  and  author K. Thürmer,  @noop  journal journal Physical review
  letters volume 105, pages 026102
  (year 2010)NoStop
[Helmholtz(1853)]helmholtz1853ueber
  author author H. v. Helmholtz, @noop  journal journal
  Annalen der Physik volume 165, pages
  353 (year 1853)NoStop
[Gouy(1910)]gouy1910constitution
  author author M. Gouy, @noop  journal journal J. Phys.
  Theor. Appl. volume 9, pages 457
  (year 1910)NoStop
[Chapman(1913)]chapman1913li
  author author D. L. Chapman, @noop  journal journal The
  London, Edinburgh, and Dublin philosophical magazine and journal of science  volume 25, pages 475 (year
  1913)NoStop
[Stern(1924)]stern1924theorie
  author author O. Stern, @noop  journal journal
  Zeitschrift für Elektrochemie und angewandte physikalische Chemie  volume 30, pages 508 (year
  1924)NoStop
[Merlet et al.(2014)Merlet,
  Limmer, Salanne, Van Roij,
  Madden, Chandler, and Rotenberg]merlet2014electric
  author author C. Merlet, author D. T. Limmer,
  author M. Salanne, author R. Van Roij, author
  P. A. Madden, author
  D. Chandler,  and author
  B. Rotenberg, @noop  journal journal The Journal of Physical Chemistry C  volume 118, pages 18291 (year 2014)NoStop
[Bard et al.(1980)Bard,
  Faulkner, Leddy, and Zoski]BardFaulkner
  author author A. J. Bard, author L. R. Faulkner,
  author J. Leddy,  and author C. G. Zoski, @noop 
  title Electrochemical methods: fundamentals and
  applications, Vol. volume 2 (publisher wiley
  New York, year 1980)NoStop
[Israelachvili(2011)]Israelachvili
  author author J. N. Israelachvili, @noop  title Intermolecular
  and surface forces (publisher Academic press, year 2011)NoStop
[Newman and Thomas-Alyea(2012)]Newman
  author author J. Newman and author K. E. Thomas-Alyea, @noop  title Electrochemical
  systems (publisher John Wiley & Sons, year
  2012)NoStop
[Kornyshev(2007)]kornyshev2007double
  author author A. A. Kornyshev, @noop  journal journal
  Journal of Physical Chemistry B  (year 2007)NoStop
[Oldham(2008)]oldham2008gouy
  author author K. B. Oldham, @noop  journal journal
  Journal of Electroanalytical Chemistry volume
  613, pages 131 (year 2008)NoStop
[Bazant, Storey, and Kornyshev(2011)]bazant2011double
  author author M. Z. Bazant, author B. D. Storey,
   and author A. A. Kornyshev,  @noop  journal journal Physical Review
  Letters volume 106, pages 046102
  (year 2011)NoStop
[Gray and Stiles(2018)]gray2018nonlinear
  author author C. Gray and author P. J. Stiles, @noop  journal journal
  European Journal of Physics  (year 2018)NoStop
[Lee et al.(2012)Lee,
  Nilson, Templeton, Griffiths,
  Kung, and Wong]lee2012comparison
  author author J. W. Lee, author R. H. Nilson,
  author J. A. Templeton, author S. K. Griffiths, author
  A. Kung,  and author
  B. M. Wong, @noop  journal journal Journal of chemical theory and computation  volume 8, pages 2012 (year
  2012)NoStop
[Lee et al.(2013)Lee,
  Templeton, Mandadapu, and Zimmerman]lee2013comparison
  author author J. W. Lee, author J. A. Templeton,
  author K. K. Mandadapu,  and  author J. A. Zimmerman,  @noop  journal journal Journal of
  Chemical Theory and Computation volume 9, pages 3051 (year 2013)NoStop
[Lee, Mani, and Templeton(2015)]lee2015atomistic
  author author J. W. Lee, author A. Mani,  and  author J. A. Templeton,  @noop  journal journal Langmuir  volume 31, pages 7496 (year
  2015)NoStop
[MacKerell Jr(2001)]mackerell2001atomistic
  author author A. D. MacKerell Jr, @noop  journal journal
  Computational Biochemistry and Biophysics , pages 7
  (year 2001)NoStop
[Hockney and Eastwood(1988)]hockney1988computer
  author author R. W. Hockney and author J. W. Eastwood, @noop  title Computer simulation
  using particles (publisher CRC Press, year
  1988)NoStop
[Perdew, Burke, and Ernzerhof(1996)]perdew1996generalized
  author author J. P. Perdew, author K. Burke,   and author M. Ernzerhof,  @noop  journal journal Physical review
  letters volume 77, pages 3865
  (year 1996)NoStop
[Tang, Sanville, and Henkelman(2009)]tang2009grid
  author author W. Tang, author E. Sanville,   and author G. Henkelman,  @noop  journal journal Journal of
  Physics: Condensed Matter volume 21, pages 084204 (year 2009)NoStop
[Heinz et al.(2008)Heinz,
  Vaia, Farmer, and Naik]heinz2008accurate
  author author H. Heinz, author R. Vaia,
  author B. Farmer,  and author R. Naik, @noop  journal journal The Journal of Physical Chemistry C  volume 112, pages 17281 (year 2008)NoStop
[Oliver, Watson, and Parker(1995)]oliver1995molecular
  author author P. Oliver, author G. Watson,   and author S. Parker,  @noop  journal journal Physical Review
  B volume 52, pages 5323 (year 1995)NoStop
[Note1()]Note1
  note Ref.+.1667emoliver1995molecular reports parameters for a Buckingham potential which we
  convert to a LJ parameterization by matching the potential well location and
  stiffness.Stop
[Mandadapu, Templeton, and  Lee(2013)]mandadapu2013polarization
  author author K. K. Mandadapu, author J. A. Templeton,  and author J. W. Lee, @noop  journal journal The
  Journal of chemical physics volume 139, pages 054115 (year 2013)NoStop
[Rossmeisl et al.(2016)Rossmeisl, Chan, Skúlason,
  Björketun, and Tripkovic]rossmeisl2016ph
  author author J. Rossmeisl, author K. Chan,
  author E. Skúlason, author M. E. Björketun,  and  author V. Tripkovic, @noop
   journal journal Catalysis Today volume 262, pages 36 (year
  2016)NoStop


