Opening gates to oxygen reduction reactions on Cu(111) surface Aslihan Sumer and Santanu Chaudhuri
Citation: J. Chem. Phys. 142, 124703 (2015); doi: 10.1063/1.4914901 View online: http://dx.doi.org/10.1063/1.4914901 View Table of Contents: http://aip.scitation.org/toc/jcp/142/12 Published by the American Institute of Physics
THE JOURNAL OF CHEMICAL PHYSICS 142, 124703 (2015)
Opening gates to oxygen reduction reactions on Cu(111) surface Aslihan Sumer and Santanu Chaudhuria) Illinois Applied Research Institute, University of Illinois at Urbana-Champaign, Champaign, Illinois 61820, USA
(Received 29 October 2014; accepted 3 March 2015; published online 24 March 2015) Electrocatalytic reduction of oxygen is composed of multiple steps, including the diffusionadsorption-dissociation of molecular oxygen. This study explores the role of electrical double layer in aqueous medium in quantifying the rate of these coupled electrochemical processes at the electrode interface during oxygen reduction. The electronic, energetic, and configurational aspects of molecular oxygen diffusion and adsorption onto Cu(111) in water are identified through density functional theory based computations. The liquid phase on Cu(111) is modeled with hexagonal-ordered water bilayers, at two slightly different structures, with O–H bonds either facing the vacuum or the metal surface. The results indicate that the energetically preferred structure of water bilayers and adsorption configuration of O2 are different in cathodic and anodic potentials. The diffusion of O2 is found to be heavily hindered at the water/metal interface because of the ordering of water molecules in bilayers as compared to the bulk liquid. The unique correlations of diffusion and adsorption kinetics with water structure identified in this work can provide clues for improving oxygen reduction efficiency. C 2015 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4914901] I. INTRODUCTION
Single crystal metal surfaces render an important testing ground for improving our understanding of mechanistic steps in key reactions observed in complex polycrystalline samples. Single crystalline Cu(111) is one such important surface where catalytic, including those of electrocatalytic, processes are of major scientific and technological importance. Abundant and a less noble metal, Cu remains unusable in key technological areas, such as fuel cell applications, for which electrocatalytic oxygen reduction reactions (ORR) play a dominant role. Electrocatalytic ORR proceeds in aqueous environment. As established through many theoretical predictions and experimental evidence, first few layers of water on single crystal surfaces are highly ordered.1–7 The rate of reactions, in general, being controlled by charge and mass transfer, these water layers participate in ORR in two key ways: (a) they provide a way to transfer electrons to the surface adsorbate molecules and thereby, making some reactions energetically more probable compared to their gas-phase analogs, and (b) they limit or promote the kinetics of transport of these molecules between the catalyst surfaces and the bulk of the aqueous phase. The impact of the structural ordering in the water layers in creating potential barriers to the transport of the molecules or in controlling the direction and rate of the electron transfer is not yet fully explored. It is expected that evaluation of the structural, chemical, and electronic properties of the water layers under electrochemical conditions would open new paths for tilting the balance in a favorable direction for ORR and identifying viable strategies for operating ORR on cheaper catalysts.
a)Author to whom correspondence should be addressed. Electronic mail:
[email protected] 0021-9606/2015/142(12)/124703/8/$30.00
Current information on the physical and chemical properties of O2-Cu interaction is deducted from experiments conducted in vacuum phase, such as those based on EELS,8,9 XPS,10,11 UPS,12 LEED,12–14 one of the reasons being the difficulty of extracting information on the reactant species on the surface when surrounded by solvent molecules with the traditional spectroscopic methods. Theoretical modeling based on density functional theory (DFT) calculations offer an alternative method to explore the characteristics of adsorption, including those at the water/metal interface, with details of electronic structure.15–21 In this work, our DFT computations aim to answer some important questions regarding mainly the role of water layers in modifying the diffusion and adsorption of O2 molecules onto the Cu(111) surface, which is critical for improving ORR kinetics. In particular, we address the following important issues: (a) Does water/Cu interface show charge transfer and form a polarized double layer under neutral conditions? (b) Is there a preferential ordering of H2O molecules on Cu(111) under cathodic conditions? (c) How do the barrier to O2 diffusion and O2 adsorption configuration on Cu(111) depend on H2O ordering? and (d) How to identify a combination of cathodic potential and water ordering favorable for both transport and oxygen adsorption? II. COMPUTATIONAL DETAILS
DFT computations were performed using DMol322 with Double Numerical plus polarization (DNP) numerical basis set, DFT semi-core pseudopotentials with 4.4 Å cutoff and 3 × 3 × 1 k-point grid. The geometry optimizations were continued till a convergence criteria of 1 × 10−5 Ha in energy, 1 × 10−3 Ha/Å in force, and 1 × 10−3 Å in displacement are reached. SCF energy was converged with the direct inversion of the iterative subspace (DIIS) technique, and the iterations were run till the largest components of the DIIS density error
142, 124703-1
© 2015 AIP Publishing LLC
124703-2
A. Sumer and S. Chaudhuri
matrix dropped less than 1 × 10−6. Perdew-Burke-Ernzerhof functional (PBE) was chosen among the available exchangecorrelation functionals implemented in DMol3, as it calculated the lattice dimension of Cu crystal as 3.61 Å, very close to the experimental Cu crystal size of 3.63 Å,23 while the other functionals underestimated it. (Calculated lattice sizes according to the functionals are BLYP:3.60 Å; BOP, VWNBP and HCTH: 3.58 Å; PW91: 3.57 Å; RPBE: 3.54 Å; LDAVWN: 3.51 Å, and LDA-PWC: 3.50 Å.) In addition, PBE was previously found to be adequate in calculations of hydrogen bond strengths in water clusters.24 Combination of PBE with numerical basis sets are particularly accurate in representing the energetics of water clusters.25 Other higher levels of theories, i.e., hybrid and van der Waals inclusive functionals26 are not yet available in the code for periodic slab computations. The periodic slab models were composed of four layers of Cu(111), the bottom two of which were kept fixed at their bulk coordinates (according to an experimental Cu crystal size of 3.63 Å), while the top two were √ fully √ relaxed in geometry optimizations. The model had 2 3 × 2 3 surface periodicity and a distance of approximately 35 Å was kept between repeating metal slabs in z-direction. There were four layers of water on the top of the metal slab with each layer composed of eight H2O molecules arranged in an order of hexagonal rings. The water layers are labeled starting from the bottom one as Layer I and so on (Figure 1). The artificial contribution to the total energy due to the non-zero dipole moment of the model in zdirection was corrected with DMol3. The energy of the O2 in vacuum was calculated in a nonperiodic fashion. The spin state of the molecule was triplet. The internal bond length was found as 1.23 Å, very close to the experimental value of 1.21 Å.27 We explored the minimum energy paths for the diffusion and adsorption of O2 in H2O/Cu(111) with constrained geometry optimizations. We relaxed all degrees of freedom in the model during the optimizations, except for the z-coordinate of one of the O atoms of O2, which was gradually decreased by 1 Å, advancing the molecule closer to the Cu(111) surface. The distance between the partially constrained O atom and the average height of the topmost layer Cu atoms is called ‘distance’ throughout the text. The models were fully relaxed
FIG. 1. H2O/Cu(111) with (a) H-up and (b) H-down models; (c) pure H2O crystal (without metal) which has the same lattice dimensions with (a) and (b) in x and y and optimized in z.
J. Chem. Phys. 142, 124703 (2015)
at a distance of approximately 2 Å at the end of which O2 was adsorbed onto Cu(111).
III. RESULTS AND DISCUSSION A. Structural and electronic properties of water layers
The water layers in the model have H2O molecules within a hexagonal ordering parallel to the Cu(111) surface. Small water clusters on Cu(111) at low temperatures were previously identified with STM as in a similar hexagonal arrangement.28 We tested the stability of our model with ab-initio molecular dynamics (AIMD) simulation at 300 K and found that the ordering within the layers and the match between the layers and metal lattice were mainly preserved during the simulation for 5 ps.29 Another AIMD simulation at 330 K30 indicated a formation of a pentagonal network of H2O molecules within the layers above Cu(111) surface. We used two computational models of H2O/Cu(111) in the geometry optimizations, with slight differences in H2O orientation, aiming to alleviate the shortcomings of a static model in representing the properties of a system dynamic in nature. In H-up model (Figure 1(a)), O–H bonds point up towards the vacuum, while in H-down model (Figure 1(b)), they point towards the metal surface. The total energies of the two models are close: H-down model is more stable than the H-up model by only 0.004 eV/AA2 (0.03 eV/cell). Previous DFT computations had shown that for Au(111), Ag(111), and Pt(111); similarly, H-down model is energetically preferred over the H-up model.31 However, repetition of the computations of H2O/Cu(111) with the dispersion correction scheme of Tkatchenko-Scheffler (TS),32 as implemented in DMol3, to compensate the shortcomings of PBE in defining the hydrogen bond strengths in water, also showed that the energy preference is reversed in favor of H-up model over H-down model by 0.001 eV/Å.2 Carrasco et al. had previously shown the positive effect of dispersion forces on the adsorption strength of H2O molecules on transition metal surfaces,33 which can be the reason of the change in energy ordering of our two H2O/Cu(111) models with TS scheme. We assume that, due to the small energetic difference between the two models, the real water layers at Cu(111) interface are composed of a mixture of H2O molecules at both types of orientations, as observed in the AIMD simulation. The average interlayer distances between the four water layers above the Cu(111) surface are found to be around 3.7 Å for both of the computational models (Table I). If we remove the metal layers from the models and optimize the lattice dimension of this pure water crystal (Figure 1(c)) in z-direction, the energetically preferred interlayer distance between the water layers becomes 3.6 Å (Table I). This shows that there is a small expansion in water layers in z-direction, by approximately 3%, in our H2O/Cu(111) as compared to pure H2O. The distance in z-direction between the center of mass of eight H2O molecules in Layer I and the center of mass of the Cu atoms at the topmost metal layer as compared to the average interlayer distances of water is smaller in the case of H-up model but larger in the case of H–down model. Inclusion of the dispersion corrections into the computations decreases
124703-3
A. Sumer and S. Chaudhuri
J. Chem. Phys. 142, 124703 (2015)
TABLE I. Average interlayer distance between four water layers and distance between metal surface and water Layer I for H2O/Cu(111) with H-up and H-down models and pure water crystal. Calculations are done with and without dispersion corrections (TS—see text). For structures see Figures 1(a)–1(c).
dH2O-H2O(avg)/Å dH2O-H2O(avg)/Å (with TS) dH2O-Cu(111)/Å dH2O-Cu(111)/Å (with TS)
H-up model
H-down model
Pure water crystal
3.71 3.61 2.96 2.85
3.70 3.59 3.91 3.58
3.61 3.48 ... ...
both the average water interlayer distances and average watermetal distances. Dispersion effects also decrease the interlayer distance in the pure water crystal, which results in an increase in water density from 0.97 g/cm3 to 1.00 g/cm3. Partitioning the charge density with the Hirshfeld scheme, we analyzed the electronic structure of the water and metal layers (Table II) and evaluated the differences in the charge transfer direction between the two models. H-up model favors electron transfer from the metal surface to the water phase while the reverse is observed for H-down model. The biggest variations in the charge density of the two models are found at the vacuum/water and water/metal interfaces, i.e., topmost metal layer and water Layers I and IV. Although not shown here, the charge states of middle and bottom metal layers in H2O/Cu(111) are similar to the same layers in vacuum phase Cu(111), indicating the diminishing effect of H2O presence on the Cu layers far from the water/metal interface. B. Energetics of O2 solvation in water
To find the solvation energy of O2 within H2O/Cu(111), we placed O2 molecules between water Layers II and III and performed geometry optimizations (Figures 2(a) and 2(b)). It should be noted that there are two layers of water between O2 at this configuration and the vacuum, which are expected to diminish the effects of the water/vacuum intersection, e.g., eliminate the shortcomings of the model due to its finite size in z-direction. At the energetically preferred configurations found after full relaxations of O2 + H2O/Cu(111), O2 is found to be at a distance of around 8-9 Å for the two models (Table III), with a spin state of triplet. The solvation energy between O2 and water can be calculated using the formula, Esolvation = EO2+H2O/Cu(111) − EH2O/Cu(111) − EO2,gas,
(1)
where EO2+H2O/Cu(111) and EH2O/Cu(111) are the energies of water/metal models with and without O2, respectively, and EO2,gas is the energy of gas phase O2. Accordingly, the solvation energy of O2 was found as 0.02 and 0.48 eV, with the H-up and H-down models, respectively. We expect that the difference between the O2 solvation energies of the two models stems from the difference in the charge state of solvated O2.34 The solvation energy of O2 calculated in a similar fashion with pure water crystal (Figure 1(c)) is found to be 0.38 eV, i.e., close to the value obtained with H-down model, and the charge state of the molecule in water crystal is the same with the one in Hdown model. Both the water layers and metal topmost layer get involved in the donation of charges to O2 in H-up model, while mainly the surrounding water Layers II and III are affected by the presence of O2, in the case of H-down one (Table II). Although the water layers in our models are perfectly ordered and, thus, are structurally different than the experimental bulk water systems (which are known to have a tetrahedral hydrogen bonding network of H2O molecules with 80% of hydrogen bonds broken35), we calculated the zero-point energy and thermal corrections (obtained with harmonic vibrational frequencies of O2 in pure water crystal) to compare the solvation enthalpies of the two computational models with that of the experimentally measured value. The enthalpy and free energy of solvation at 298 K are found to be 0.20 and 0.40 eV lower, respectively, than the electronic solvation energy. Thus, the enthalpy of solvation is found to be −0.18 and 0.28 eV, for the H-up and H-down models, respectively. The free energies of solvation, calculated in a similar fashion, are −0.38 and 0.08 eV, for the two models. Thus, solvation is thermodynamically stable in H-up water model whereas slightly unstable in H-down model. The experimentally measured solvation enthalpy of O2 in bulk liquid water, which is −0.12 eV,36 is thus
TABLE II. Charges (in electronic units, e, from Hirshfeld scheme) of water and metal layers of H2O/Cu(111) and O2 + H2O/Cu(111) with O2 either between water layers II and III, or adsorbed on Cu(111). For structures, see Figures 1(a), 1(b), 2(a), 2(b), 2(e), and 2(f).
Layers Water layer IV Water layer III Water layer II Water layer I Top Cu layer Sum of two inner Cu layers Bottom Cu layer
H2O/Cu(111) H-up model
H2O/Cu(111) H-down model
O2,sol + H2O/Cu(111) H-up model
O2,sol + H2O/Cu(111) H-down model
O2,ads + H2O/Cu(111) H-up model
O2,ads + H2O/Cu(111) H-down model
−0.10 −0.01 0.00 −0.19 0.22 0.16
0.09 0.01 0.00 0.00 −0.17 0.15
0.00 −0.06 0.05 −0.15 0.27 0.15
0.08 −0.03 −0.08 0.03 −0.18 0.14
−0.07 −0.01 0.01 −0.22 0.64 0.15
0.12 −0.03 0.02 0.11 0.14 0.14
−0.08
−0.08
−0.08
−0.08
−0.09
−0.08
124703-4
A. Sumer and S. Chaudhuri
J. Chem. Phys. 142, 124703 (2015)
FIG. 2. Side and top views of O2 + H2O/Cu(111) configurations: Fully relaxed with O2 located between H2O molecules in Layers II and III, with (a) H-up model (b) H-down model; constrained with O2 kept at approximately 2.2 Å from Cu(111) with (c) H-up model (d) H-down model; and fully relaxed with O2 adsorbed on Cu(111) with (e) H-up model (f) H-down model.
between the solvation enthalpies of O2 in the two H2O/Cu(111) models. It should finally be noted that inclusion of the dispersion effects did not cause a major change in the O2 solvation energy in H2O/Cu(111) or pure water crystal (Table III), and thus, it is concluded that exclusion of these would cause no major effect on the relevant conclusions obtained from rest of the computations. C. Energetics of O2 diffusion in water towards Cu(111) surface
Starting from the coordinates of previously mentioned fully relaxed O2 + H2O/Cu(111) solvation system, we gradually decreased the distance of O2 to Cu(111) and performed constrained geometry optimizations, to find the minimum energy path for diffusion and adsorption of O2 at water/metal interface. In addition, two fully relaxed optimizations were performed in order to identify the stable interstitial positions of O2 between the H2O molecules in Layers I and II and the stable adsorption configuration on Cu(111) surface. Total energies obtained in this way as a function of O2 distance for H-up and H-down models are plotted in Figures 3(a) and 3(b), respectively, where the values are referenced to the initial total
energy of solvation system in order to make it easier to follow the energy trajectory for the O2 molecule approaching the metal surface. For both water models, the plots of energy have two peaks, which correspond to the configurations at which O2 is in the middle of hexagonal H2O rings, i.e., when they are “passing” through the water Layers II and I, in order. The energy barriers that are associated with the first of these peaks correspond to 0.36 and 0.32, for H-up and H-down models, respectively. Thus, the height of the first energy barriers is very close for the two models and in addition, they are close to the height of the diffusion barrier calculated with the pure water crystal (Figure 1(c)), at 0.34 eV, indicating the relatively small influence of water/metal interface on O2-H2O interaction at this distance. It should be noted here that the barriers of diffusion through ordered water layers in our models are higher than the experimentally measured barriers of O2 diffusion in bulk water, at 0.20 eV.37 Unlike the first ones, the second energy barriers, calculated as the energy difference of O2 + H2O/Cu(111) with O2 passing through water Layer I and relaxed at the interstitial position between Layers I and II, as 0.43 and 0.34 eV, for Hup and H-down models, respectively (Figures 3(a) and 3(b)), are found to be strongly dependent on the orientation of H2O with respect to the metal surface. Using these energy barriers
TABLE III. Solvation energy, charge (in electronic units, e, from Hirshfeld scheme) and distance of O2 to Cu(111) for fully relaxed O2 + H2O/Cu(111) with O2 between water layers II and III, with H-up and H-down models (Figures 2(a) and 2(b)) and in pure water crystal (Figure 1(c)). Calculations are done with and without dispersion corrections (TS—see text).
Esolvation/eV Esolvation/eV (with TS) q q (with TS) d/Å d/Å (with TS)
H-up model
H-down model
Water crystal
0.02 0.10 −0.19 −0.19 8.03 8.06
0.48 0.50 0.11 0.11 8.88 8.92
0.38 0.31 0.11 0.11 ... ...
124703-5
A. Sumer and S. Chaudhuri
J. Chem. Phys. 142, 124703 (2015)
become 1.29 × 105 and 4.09 × 106 s−1, respectively. Assuming the O2 concentration at the H2O/Cu(111) interface to be identical to the O2 solubility in bulk water at 298 K and 1 atm, i.e., 8.3 mg/l, the rate of adsorptions for H-up and Hdown models are found as 7.29 × 1018 and 2.30 × 1020 m−2 s−1, respectively. For comparison, the rate of adsorption on Cu(111) in vacuum, as calculated from r ads = √
P
S(θ),
2πmkT
(5)
using the partial pressure of O2 in air at 1 atm and the initial sticking coefficient, S(0), of 0.001,13 is 5.71 × 1023 m−2 s−1, which is higher than the adsorption rate of O2 in solution because of the lack of barriers to the adsorption in vacuum. It should also be noted that as an alternative path to molecular diffusion and then adsorption onto Cu(111), we also tested the possibility of chemical transformation of O2 to OOH in solution by exploring the minimum energy path and the energetics of the reaction.29 We identified an energy barrier of 0.50 eV for reduction, which is higher than the barrier for diffusion and adsorption. The reaction was also found to be endothermic, indicating an overall lack of thermodynamic and kinetic driving force for chemical transformation of O2 to OOH in water. D. Structural and energetic properties of O2 adsorption over Cu(111)
FIG. 3. Relative energies of constrained (•) and fully relaxed () O2 + H2O/Cu(111) configurations as a function of distance of O2 to Cu(111) with (a) H-up and (b) H-down models. Energy barriers of O2 diffusion through H2O molecules in Layer I are also shown.
as activation energies for adsorption, Eact, in an Arrhenius form of rate equation, k = k 0e
−E act kT
,
(2)
it can be shown that there is a difference in the adsorption rate constants, k, between the two models at 298 K, by H −up
H −down) −(E act −E act k H −up 1 kT =e = , (3) H −down 31 k meaning a higher rate of adsorption of O2 molecule on Cu(111) when the water/metal interface is composed of H-down oriented H2O molecules. It is possible to calculate the pre-exponential factor, k0, from initial state f (4) k 0 = transition state f
as 2.62 × 1012 with the vibrational frequencies, f, calculated, for simplicity, in pure water crystal. Then, adsorption rate constant, k, for H-up and H-down models (kH−up and kH−down)
Starting from the outcome of the final constrained optimization of O2 + H2O/Cu(111) with O2 at a distance of approximately 2.2 Å (Figures 2(c) and 2(d)), a fully relaxed optimization is performed during which O2 forms bonds with Cu atoms and gets adsorbed on Cu(111). During the relaxation and bond forming between O2 and Cu atoms, O2 changes its spin state from triplet to singlet. The crossing of the spin states happens at approximately 1.9-2.0 Å distance to the metal surface. It was previously found using DFT calculations that a similar decrease in the spin state of the adsorption complex occurs at Cun clusters with n as small as 5 atoms.15 The final adsorption configuration of O2 over Cu(111) is found to be dependent on the computational model used: It is obtained as three-fold with O-O molecular axis perpendicular to the Cu(111) surface, and three-fold with O-O axis parallel to Cu(111) surface with Hup and H-down models, respectively (Figures 2(e) and 2(f), Table IV). DFT computations on vacuum phase O2/Cu(111) gave the energetically most preferred adsorption configuration of O2 on Cu(111) in absence of H2O molecules as parallel oriented on four-fold site (Table IV), while a parallel oriented threefold adsorption configuration, which was previously identified as the preferred adsorption configuration in vacuum,18 is only 1 meV higher in energy. We calculated the adsorption energies for the final O2 + H2O/Cu(111) adsorption configurations identified with the two models with Eads = EO2,ads+H2O/Cu(111) − EH2O/Cu(111) − EO2,gas.
(6)
Adsorption energies of O2/Cu(111) in vacuum, for comparison, are calculated with
124703-6
A. Sumer and S. Chaudhuri
J. Chem. Phys. 142, 124703 (2015)
TABLE IV. Top views, adsorption sites, energies and charges (in electronic units, e, from Hirshfeld scheme) of O2 for O2 + H2O/Cu(111) and O2/Cu(111) (vacuum phase). Perpendicular O2
Parallel O2
Three-fold (fcc) −0.25 −0.34
Three-fold (hcp) −0.78 −0.42
Three-fold (fcc) −0.14 −0.32
Four-fold −0.79 −0.40
O2 + H2O/Cu(111)
Adsorption site Eads/eV q O2/Cu(111)
Adsorption site Eads/eV q
Eads = EO2/Cu(111) − ECu(111) − EO2,gas.
(7)
Water does not cause a major change in the stability of parallel oriented O2, while it decreases the adsorption energy of perpendicularly oriented one (Table IV). In both water and vacuum, parallel oriented configuration is more stable than perpendicular one, and the difference in the adsorption energies is larger in vacuum than in water. The identification of the adsorption configuration of O2 as perpendicular to the metal surface in the case of H-up model is a result of the interaction between O2 and H2O through hydrogen bonding, the degree of which is dependent on the
H2O orientations with respect to the metal surface. As it was mentioned before, the distance between the topmost layer Cu atoms and water Layer I is smaller at H-up model than at Hdown one. Consequently, when O2 gets adsorbed on Cu(111) in water, H2O’s in close vicinity at the H-up model hinder its full relaxation to an energetically more stable adsorption configuration, whereas the molecule is, relatively, under less influence of hydrogen bonding at the H-down model. The influence of surrounding H2O’s on perpendicularly adsorbed O2 can also be evaluated through a comparative analysis of local density of states (LDOS) graphs. LDOS at upper atom of perpendicular O2 and H atoms of H-up oriented H2O molecules surrounding it have common features between −23 and −20 eV (Figures 4(a) and 4(b)), indicating the presence of hydrogen bonds between the molecules. For comparison, we also optimized the adsorption geometry of perpendicular O2 on Cu(111) in water with the H-down model and also, in vacuum, and found that LDOS at O2 computed with these two last models are very similar to each other and thus, prove the lack of hydrogen bonds between H2O’s and O2 in H-down model (Figures 4(c) and 4(d)). As a final note, we also checked the influence of lower water density in our models on the final adsorption configuration and found no effect of increasing water density on the final conclusions.29 Lower energy of parallel oriented O2 as compared to the perpendicularly oriented one points out to an energetic driving force for reorientation of the molecule even if its initial adsorption configuration after diffusion is perpendicular. So, in a separate work, we computed the energetics of possible reorientation and chemical transformation paths for O2 over Cu(111) in water as modeled with H-up model and found that the energy barrier for reorientation of perpendicular O2 is higher than the energy barrier for the reduction of the molecule to OOH.38 These
FIG. 4. LDOS at (a) O of perpendicularly adsorbed O2 with H-up model; (b) H of H2O molecules neighboring O2 at model in (a); (c) O of perpendicularly adsorbed O2 with H-down model; (d) O of perpendicularly adsorbed O2 in vacuum.
124703-7
A. Sumer and S. Chaudhuri
findings overall indicate that on the Cu(111) surface, along with thermodynamically and kinetically more favorable threefold adsorbed parallel oriented O2 species, it is possible to observe a less stable perpendicularly adsorbed O2, and the reduction reaction may thus proceed via two separate paths starting from these two different adsorption configurations. E. Adsorption induced electronic changes and effects of an electric field
Through an integration of the charge along the cell length (Figures 5(a) and 5(b)), we calculated the dipole moment of H2O/Cu(111) with H-up and H-down models as 6.51 and −11.81 D (corresponding to 0.095 and −0.172 D/Å2), respectively. The difference in the directions of dipole causes an opposite response in the energy of the two models to variations in an external electric field,29 as was also shown previously for other metal surfaces.39–41 It can be speculated that gradual reorientation of the H2O molecules following the variations in electric field would eventually lead to formation of homogeneous H-up and H-down layers at anodic and cathodic potentials, respectively, assuming the hexagonal ordering in the water layers is preserved duringthe potential changes. Our
J. Chem. Phys. 142, 124703 (2015)
computations indicate that this structuring in water layers at the metal interface has an influence on the energetics of O2 diffusion. Thus, it is expected that as the magnitude of the electric field towards the metal surface increases, i.e., as potentials shift to cathodic values, the rate of O2 diffusion and the coverage of O2 on metal surface would also increase. In the presence of an increasing electric field towards the vacuum, i.e., at anodic electrode potentials, the total coverage of O2 on metal surface would decrease and the ratio of perpendicularly adsorbed O2 to parallel adsorbed one would increase. In parallel to the indirect influence of an electric field through water structuring at the metal interface, the charge transfer from the metal to the adsorbate and the consequent changes in the dipole of the system during adsorption also causes the adsorption energetics to respond to the variations in potential of an electrocatalytic system. The charge density comparison, based on Hirshfeld partitioning (Table II) and also plots of charge density difference (between optimized charge density of the system and the sum of the densities of separate atoms that make up system), averaged in x and y dimensions (Figures 5(a) and 5(b)), before and after O2 adsorption indicate an adsorption induced charge redistribution at the water layers, outside region of the water/metal interface where the electric field is directly influential in an electrocatalytic system. Thus, in order to explicitly evaluate the electric field induced changes in O-metal bond formation energetics, we calculated the dipole of the O2/Cu(111) system in vacuum. We found an adsorption induced dipole moment for O2/Cu(111), for both fourfold and perpendicular adsorption configurations, as −1.32 D, which hints a stabilization (destabilization) at cathodic (anodic) potentials. For O2 adsorbed on either Pt nanoclusters or Pt(111) surface, it was shown previously that the energetic response to an external homogeneous electric field (F) varying between ±0.5 V/Å follows a first order Stark equilibrium, ∆E = −µF,
(8)
where µ is the adsorption induced dipole moment without the field.42 For the O2/Cu(111) system, this equilibrium implies a change in adsorption energy of both fourfold and perpendicular configurations by ±0.14 eV in external electric fields of ±0.5 V/Å. Adsorption induced dipole of Cu(111) is higher than the dipole of O2/Pt(111), calculated previously as −0.36 D, due to the charge state of the adsorbate on Pt, which is as low as −0.09.43 Hirshfeld charge analysis of O2/Cu(111), on the other hand, indicates a relatively larger net charge on the adsorbate, ∼ −0.3-0.4 (Table IV), explaining the calculated relatively larger dipole of the system.
IV. CONCLUSIONS
FIG. 5. Plane (xy) averaged charge density difference at H2O/Cu(111) before (dashed line) and after (solid line) O2 adsorption. (a) Perpendicularly adsorbed O2 with H-up model, (b) parallel adsorbed O2 with H-down model.
We computed the energetics of diffusion in water and surface adsorption for O2 molecule at the water/Cu(111) interface. We found that adsorption rate and configuration of O2 and the structure of water network at the metal interface have some unique correlations which can be exploited using electrostatic potentials for favorable ORR kinetics. So-called H-up oriented H2O molecules cause higher barriers for the diffusion
124703-8
A. Sumer and S. Chaudhuri
of O2 towards Cu(111), as compared to opposite oriented Hdown H2O molecules, and thus, H-up ordering reduces the adsorption rate. As the potentials decrease from positive to negative, i.e., from anodic to cathodic values, H-down ordering becomes energetically more stable, causing the H-up oriented H2O molecules flip over. This indicates that O2 adsorption rate on Cu(111) will increase under cathodic conditions. Along with O2 surface coverage, the potential induced structuring in water network also influences the adsorption configuration of O2 on Cu(111). It is shown that while the thermodynamically most stable parallel adsorption configuration is also kinetically favorable at cathodic potentials, a perpendicular adsorption configuration is expected to increase its relative coverage at anodic potentials. Analysis of the local density of states indicated the role of hydrogen bonds in the adsorption of O2 at this configuration. The calculated dipole moment of the O2 on Cu(111) at both configurations was the same and higher than was found previously for O2/Pt(111). The fundamental barrier to efficient ORR on cheaper and more abundant metals such as Cu(111) is primarily due to lack of understanding of the double layer under electrostatic potentials. The complexity of the electrical double layer still needs further investigation such as a realistic polarization scheme for the water layer and ionic species adsorption, e.g., NaCl, on Cu-surface, which can affect the water structure at metal interface.44 Thermal effects on the structural, electronic, and chemical properties of the system would influence the interaction between water and Cu(111) surfaces, so are needed to be included with molecular dynamic simulations. For a water double layer without charged species, the DFT results presented in this work for a planar Cu(111) surface reveal important details on the role of water bilayer in controlling the path and adsorption kinetics of O2. Further spectroscopic measurements such as sum frequency generation spectroscopy will be important for identifying the configuration of the adsorbate in anodic/cathodic potentials identified in this work. ACKNOWLEDGMENTS
The authors acknowledge Joseph Osborne from Boeing Company for his thoughtful suggestions and interest in ORR processes on Cu surface. S.C. acknowledges funding from Joint Center for Aerospace Technology Innovation (JCATI) and Boeing Company. 1P. A. Thiel and T. E. Madey, Surf. Sci. Rep. 7, 211 (1987);
ibid. 46, 1 (2002). Meng, L. F. Xu, E. G. Wang, and S. W. Gao, Phys. Rev. Lett. 89, 176104 (2002). 3M. Morgenstern, J. Muller, T. Michely, and G. Comsa, Z. Phys. Chem. 198, 43 (1997).
2S.
J. Chem. Phys. 142, 124703 (2015) 4K.
Jacobi, K. Bedürftig, Y. Wang, and G. Ertl, Surf. Sci. 472, 9 (2001).
5A. Glebov, A. P. Graham, A. Menzel, and J. P. Toennies, J. Chem. Phys. 106,
9382 (1997). Haq, J. Harnett, and A. Hodgson, Surf. Sci. 505, 171 (2002). 7S. Standop, A. Redinger, M. Morgenstern, T. Michely, and C. Busse, Phys. Rev. B 82, 161412 (2010). 8T. Sueyoshi, T. Sasaki, and Y. Iwasawa, Surf. Sci. 365, 310 (1996). 9K. Prabhakaran, Surf. Sci. 225, L25 (1990). 10M. K. Rajumon, K. Prabhakaran, and C. N. R. Rao, Surf. Sci. 233, L237 (1990). 11G. R. Gruzalski, D. M. Zehner, and J. F. Wendelken, Surf. Sci. 159, 353 (1985). 12A. Spitzer and H. Lüth, Surf. Sci. 118, 136 (1982). 13F. H. P. M. Habraken, C. M. A. M. Mesters, and G. A. Bootsma, Surf. Sci. 97, 264 (1980). 14R. W. Judd, P. Hollins, and J. Pritchard, Surf. Sci. 171, 643 (1986). 15X. Yuan, L. Liu, X. Wang, M. Yang, K. A. Jackson, and J. Jellinek, J. Phys. Chem. A 115, 8705 (2011). 16A. Soon, M. Todorova, B. Delley, and C. Stamp, Phys. Rev. B 73, 165424 (2006). 17L. Padilla-Campos, J. Mol. Struct.: THEOCHEM 815, 63 (2007). 18Y. Xu and M. Mavrikakis, Surf. Sci. 494, 131 (2001). 19M. Alatalo, A. Puisto, H. Pitkänen, A. S. Foster, and K. Laasonen, Surf. Sci. 600, 1574 (2006). 20M. Neurock, R. A. Van santen, W. Biemolt, and A. P. J. Jansen, J. Am. Chem. Soc. 116, 6860 (1994). 21C. D. Taylor, R. G. Kelly, and M. Neurock, J. Electrochem. Soc. 154, F217 (2007). 22B. Delley, J. Chem. Phys. 113, 7756 (2000). 23J. T. Drotar, T.-M. Lu, and G.-C. Wang, J. Appl. Phys. 96, 7071 (2004). 24B. Santra, A. Michaelides, and M. Scheffler, J. Chem. Phys. 127, 184104 (2007). 25F. Li, L. Wang, J. Zhao, J. R. Xie, K. E. Riley, and Z. Chen, Theor. Chem. Acc. 130, 341 (2011). 26J. Carrasco, B. Santra, J. Klimes, and A. Michaelides, Phys. Rev. Lett. 106, 026101 (2011). 27CRC Handbook of Chemistry and Physics, edited by D. R. Lide (CRC Press, Boca Raton, FL, 1993). 28A. Michaelides and K. Morgenstern, Nat. Mater. 6, 597 (2007). 29See supplementary material at http://dx.doi.org/10.1063/1.4914901 for the static and dynamics computations on structure of water, effects of water density and electric field. 30R. Nadler and J. Fernandez Sanz, J. Mol. Model 18, 2433 (2012). 31S. Schnur and A. Groß, New J. Phys. 11, 125003 (2009). 32A. Tkatchenko and M. Scheffler, Phys. Rev. Lett. 102, 073005 (2009). 33J. Carrasco, J. Klimeš, and A. Michaelides, J. Chem. Phys. 138, 024708 (2013). 34R. M. Lynden-Bell and J. Rasaiah, J. Chem. Phys. 107, 1981 (1997). 35Z. Li, “Improving ab initio molecular dynamics of liquid water,” Ph.D. dissertation (Department of Physics, Princeton University, 2012). 36M. H. Abraham, J. Am. Chem. Soc. 104, 2085 (1982). 37C. F. Wong and W. Hayduk, Can. J. Chem. Eng. 68, 849 (1990). 38A. Sumer and S. Chaudhuri, “Modeling aqueous phase electro-reduction of O2 on Cu(111) using density functional based computations” (unpublished). 39C. G. Sanchez, Surf. Sci. 527, 1 (2003). 40A. Roudgar and A. Groß, Chem. Phys. Lett. 409, 157 (2005). 41G. L. Borges, J. G. Gordon, O. R. Melroy, D. G. Wiesler, D. Yee, and L. B. Sorensen, Nature 368, 444 (1994). 42M. P. Hyman and J. W. Medlin, J. Phys. Chem. B 109, 6304 (2005). 43L. Qi, X. Qian, and J. Li, Phys. Rev. Lett. 101, 146101 (2008). 44C. G. Guymon, R. L. Rowley, J. N. Harb, and D. R. Wheeler, Condens. Matter Phys. 8, 335 (2005). 6S.
Citation: J. Chem. Phys. 142, 124703 (2015); doi: 10.1063/1.4914901 View online: http://dx.doi.org/10.1063/1.4914901 View Table of Contents: http://aip.scitation.org/toc/jcp/142/12 Published by the American Institute of Physics
THE JOURNAL OF CHEMICAL PHYSICS 142, 124703 (2015)
Opening gates to oxygen reduction reactions on Cu(111) surface Aslihan Sumer and Santanu Chaudhuria) Illinois Applied Research Institute, University of Illinois at Urbana-Champaign, Champaign, Illinois 61820, USA
(Received 29 October 2014; accepted 3 March 2015; published online 24 March 2015) Electrocatalytic reduction of oxygen is composed of multiple steps, including the diffusionadsorption-dissociation of molecular oxygen. This study explores the role of electrical double layer in aqueous medium in quantifying the rate of these coupled electrochemical processes at the electrode interface during oxygen reduction. The electronic, energetic, and configurational aspects of molecular oxygen diffusion and adsorption onto Cu(111) in water are identified through density functional theory based computations. The liquid phase on Cu(111) is modeled with hexagonal-ordered water bilayers, at two slightly different structures, with O–H bonds either facing the vacuum or the metal surface. The results indicate that the energetically preferred structure of water bilayers and adsorption configuration of O2 are different in cathodic and anodic potentials. The diffusion of O2 is found to be heavily hindered at the water/metal interface because of the ordering of water molecules in bilayers as compared to the bulk liquid. The unique correlations of diffusion and adsorption kinetics with water structure identified in this work can provide clues for improving oxygen reduction efficiency. C 2015 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4914901] I. INTRODUCTION
Single crystal metal surfaces render an important testing ground for improving our understanding of mechanistic steps in key reactions observed in complex polycrystalline samples. Single crystalline Cu(111) is one such important surface where catalytic, including those of electrocatalytic, processes are of major scientific and technological importance. Abundant and a less noble metal, Cu remains unusable in key technological areas, such as fuel cell applications, for which electrocatalytic oxygen reduction reactions (ORR) play a dominant role. Electrocatalytic ORR proceeds in aqueous environment. As established through many theoretical predictions and experimental evidence, first few layers of water on single crystal surfaces are highly ordered.1–7 The rate of reactions, in general, being controlled by charge and mass transfer, these water layers participate in ORR in two key ways: (a) they provide a way to transfer electrons to the surface adsorbate molecules and thereby, making some reactions energetically more probable compared to their gas-phase analogs, and (b) they limit or promote the kinetics of transport of these molecules between the catalyst surfaces and the bulk of the aqueous phase. The impact of the structural ordering in the water layers in creating potential barriers to the transport of the molecules or in controlling the direction and rate of the electron transfer is not yet fully explored. It is expected that evaluation of the structural, chemical, and electronic properties of the water layers under electrochemical conditions would open new paths for tilting the balance in a favorable direction for ORR and identifying viable strategies for operating ORR on cheaper catalysts.
a)Author to whom correspondence should be addressed. Electronic mail:
[email protected] 0021-9606/2015/142(12)/124703/8/$30.00
Current information on the physical and chemical properties of O2-Cu interaction is deducted from experiments conducted in vacuum phase, such as those based on EELS,8,9 XPS,10,11 UPS,12 LEED,12–14 one of the reasons being the difficulty of extracting information on the reactant species on the surface when surrounded by solvent molecules with the traditional spectroscopic methods. Theoretical modeling based on density functional theory (DFT) calculations offer an alternative method to explore the characteristics of adsorption, including those at the water/metal interface, with details of electronic structure.15–21 In this work, our DFT computations aim to answer some important questions regarding mainly the role of water layers in modifying the diffusion and adsorption of O2 molecules onto the Cu(111) surface, which is critical for improving ORR kinetics. In particular, we address the following important issues: (a) Does water/Cu interface show charge transfer and form a polarized double layer under neutral conditions? (b) Is there a preferential ordering of H2O molecules on Cu(111) under cathodic conditions? (c) How do the barrier to O2 diffusion and O2 adsorption configuration on Cu(111) depend on H2O ordering? and (d) How to identify a combination of cathodic potential and water ordering favorable for both transport and oxygen adsorption? II. COMPUTATIONAL DETAILS
DFT computations were performed using DMol322 with Double Numerical plus polarization (DNP) numerical basis set, DFT semi-core pseudopotentials with 4.4 Å cutoff and 3 × 3 × 1 k-point grid. The geometry optimizations were continued till a convergence criteria of 1 × 10−5 Ha in energy, 1 × 10−3 Ha/Å in force, and 1 × 10−3 Å in displacement are reached. SCF energy was converged with the direct inversion of the iterative subspace (DIIS) technique, and the iterations were run till the largest components of the DIIS density error
142, 124703-1
© 2015 AIP Publishing LLC
124703-2
A. Sumer and S. Chaudhuri
matrix dropped less than 1 × 10−6. Perdew-Burke-Ernzerhof functional (PBE) was chosen among the available exchangecorrelation functionals implemented in DMol3, as it calculated the lattice dimension of Cu crystal as 3.61 Å, very close to the experimental Cu crystal size of 3.63 Å,23 while the other functionals underestimated it. (Calculated lattice sizes according to the functionals are BLYP:3.60 Å; BOP, VWNBP and HCTH: 3.58 Å; PW91: 3.57 Å; RPBE: 3.54 Å; LDAVWN: 3.51 Å, and LDA-PWC: 3.50 Å.) In addition, PBE was previously found to be adequate in calculations of hydrogen bond strengths in water clusters.24 Combination of PBE with numerical basis sets are particularly accurate in representing the energetics of water clusters.25 Other higher levels of theories, i.e., hybrid and van der Waals inclusive functionals26 are not yet available in the code for periodic slab computations. The periodic slab models were composed of four layers of Cu(111), the bottom two of which were kept fixed at their bulk coordinates (according to an experimental Cu crystal size of 3.63 Å), while the top two were √ fully √ relaxed in geometry optimizations. The model had 2 3 × 2 3 surface periodicity and a distance of approximately 35 Å was kept between repeating metal slabs in z-direction. There were four layers of water on the top of the metal slab with each layer composed of eight H2O molecules arranged in an order of hexagonal rings. The water layers are labeled starting from the bottom one as Layer I and so on (Figure 1). The artificial contribution to the total energy due to the non-zero dipole moment of the model in zdirection was corrected with DMol3. The energy of the O2 in vacuum was calculated in a nonperiodic fashion. The spin state of the molecule was triplet. The internal bond length was found as 1.23 Å, very close to the experimental value of 1.21 Å.27 We explored the minimum energy paths for the diffusion and adsorption of O2 in H2O/Cu(111) with constrained geometry optimizations. We relaxed all degrees of freedom in the model during the optimizations, except for the z-coordinate of one of the O atoms of O2, which was gradually decreased by 1 Å, advancing the molecule closer to the Cu(111) surface. The distance between the partially constrained O atom and the average height of the topmost layer Cu atoms is called ‘distance’ throughout the text. The models were fully relaxed
FIG. 1. H2O/Cu(111) with (a) H-up and (b) H-down models; (c) pure H2O crystal (without metal) which has the same lattice dimensions with (a) and (b) in x and y and optimized in z.
J. Chem. Phys. 142, 124703 (2015)
at a distance of approximately 2 Å at the end of which O2 was adsorbed onto Cu(111).
III. RESULTS AND DISCUSSION A. Structural and electronic properties of water layers
The water layers in the model have H2O molecules within a hexagonal ordering parallel to the Cu(111) surface. Small water clusters on Cu(111) at low temperatures were previously identified with STM as in a similar hexagonal arrangement.28 We tested the stability of our model with ab-initio molecular dynamics (AIMD) simulation at 300 K and found that the ordering within the layers and the match between the layers and metal lattice were mainly preserved during the simulation for 5 ps.29 Another AIMD simulation at 330 K30 indicated a formation of a pentagonal network of H2O molecules within the layers above Cu(111) surface. We used two computational models of H2O/Cu(111) in the geometry optimizations, with slight differences in H2O orientation, aiming to alleviate the shortcomings of a static model in representing the properties of a system dynamic in nature. In H-up model (Figure 1(a)), O–H bonds point up towards the vacuum, while in H-down model (Figure 1(b)), they point towards the metal surface. The total energies of the two models are close: H-down model is more stable than the H-up model by only 0.004 eV/AA2 (0.03 eV/cell). Previous DFT computations had shown that for Au(111), Ag(111), and Pt(111); similarly, H-down model is energetically preferred over the H-up model.31 However, repetition of the computations of H2O/Cu(111) with the dispersion correction scheme of Tkatchenko-Scheffler (TS),32 as implemented in DMol3, to compensate the shortcomings of PBE in defining the hydrogen bond strengths in water, also showed that the energy preference is reversed in favor of H-up model over H-down model by 0.001 eV/Å.2 Carrasco et al. had previously shown the positive effect of dispersion forces on the adsorption strength of H2O molecules on transition metal surfaces,33 which can be the reason of the change in energy ordering of our two H2O/Cu(111) models with TS scheme. We assume that, due to the small energetic difference between the two models, the real water layers at Cu(111) interface are composed of a mixture of H2O molecules at both types of orientations, as observed in the AIMD simulation. The average interlayer distances between the four water layers above the Cu(111) surface are found to be around 3.7 Å for both of the computational models (Table I). If we remove the metal layers from the models and optimize the lattice dimension of this pure water crystal (Figure 1(c)) in z-direction, the energetically preferred interlayer distance between the water layers becomes 3.6 Å (Table I). This shows that there is a small expansion in water layers in z-direction, by approximately 3%, in our H2O/Cu(111) as compared to pure H2O. The distance in z-direction between the center of mass of eight H2O molecules in Layer I and the center of mass of the Cu atoms at the topmost metal layer as compared to the average interlayer distances of water is smaller in the case of H-up model but larger in the case of H–down model. Inclusion of the dispersion corrections into the computations decreases
124703-3
A. Sumer and S. Chaudhuri
J. Chem. Phys. 142, 124703 (2015)
TABLE I. Average interlayer distance between four water layers and distance between metal surface and water Layer I for H2O/Cu(111) with H-up and H-down models and pure water crystal. Calculations are done with and without dispersion corrections (TS—see text). For structures see Figures 1(a)–1(c).
dH2O-H2O(avg)/Å dH2O-H2O(avg)/Å (with TS) dH2O-Cu(111)/Å dH2O-Cu(111)/Å (with TS)
H-up model
H-down model
Pure water crystal
3.71 3.61 2.96 2.85
3.70 3.59 3.91 3.58
3.61 3.48 ... ...
both the average water interlayer distances and average watermetal distances. Dispersion effects also decrease the interlayer distance in the pure water crystal, which results in an increase in water density from 0.97 g/cm3 to 1.00 g/cm3. Partitioning the charge density with the Hirshfeld scheme, we analyzed the electronic structure of the water and metal layers (Table II) and evaluated the differences in the charge transfer direction between the two models. H-up model favors electron transfer from the metal surface to the water phase while the reverse is observed for H-down model. The biggest variations in the charge density of the two models are found at the vacuum/water and water/metal interfaces, i.e., topmost metal layer and water Layers I and IV. Although not shown here, the charge states of middle and bottom metal layers in H2O/Cu(111) are similar to the same layers in vacuum phase Cu(111), indicating the diminishing effect of H2O presence on the Cu layers far from the water/metal interface. B. Energetics of O2 solvation in water
To find the solvation energy of O2 within H2O/Cu(111), we placed O2 molecules between water Layers II and III and performed geometry optimizations (Figures 2(a) and 2(b)). It should be noted that there are two layers of water between O2 at this configuration and the vacuum, which are expected to diminish the effects of the water/vacuum intersection, e.g., eliminate the shortcomings of the model due to its finite size in z-direction. At the energetically preferred configurations found after full relaxations of O2 + H2O/Cu(111), O2 is found to be at a distance of around 8-9 Å for the two models (Table III), with a spin state of triplet. The solvation energy between O2 and water can be calculated using the formula, Esolvation = EO2+H2O/Cu(111) − EH2O/Cu(111) − EO2,gas,
(1)
where EO2+H2O/Cu(111) and EH2O/Cu(111) are the energies of water/metal models with and without O2, respectively, and EO2,gas is the energy of gas phase O2. Accordingly, the solvation energy of O2 was found as 0.02 and 0.48 eV, with the H-up and H-down models, respectively. We expect that the difference between the O2 solvation energies of the two models stems from the difference in the charge state of solvated O2.34 The solvation energy of O2 calculated in a similar fashion with pure water crystal (Figure 1(c)) is found to be 0.38 eV, i.e., close to the value obtained with H-down model, and the charge state of the molecule in water crystal is the same with the one in Hdown model. Both the water layers and metal topmost layer get involved in the donation of charges to O2 in H-up model, while mainly the surrounding water Layers II and III are affected by the presence of O2, in the case of H-down one (Table II). Although the water layers in our models are perfectly ordered and, thus, are structurally different than the experimental bulk water systems (which are known to have a tetrahedral hydrogen bonding network of H2O molecules with 80% of hydrogen bonds broken35), we calculated the zero-point energy and thermal corrections (obtained with harmonic vibrational frequencies of O2 in pure water crystal) to compare the solvation enthalpies of the two computational models with that of the experimentally measured value. The enthalpy and free energy of solvation at 298 K are found to be 0.20 and 0.40 eV lower, respectively, than the electronic solvation energy. Thus, the enthalpy of solvation is found to be −0.18 and 0.28 eV, for the H-up and H-down models, respectively. The free energies of solvation, calculated in a similar fashion, are −0.38 and 0.08 eV, for the two models. Thus, solvation is thermodynamically stable in H-up water model whereas slightly unstable in H-down model. The experimentally measured solvation enthalpy of O2 in bulk liquid water, which is −0.12 eV,36 is thus
TABLE II. Charges (in electronic units, e, from Hirshfeld scheme) of water and metal layers of H2O/Cu(111) and O2 + H2O/Cu(111) with O2 either between water layers II and III, or adsorbed on Cu(111). For structures, see Figures 1(a), 1(b), 2(a), 2(b), 2(e), and 2(f).
Layers Water layer IV Water layer III Water layer II Water layer I Top Cu layer Sum of two inner Cu layers Bottom Cu layer
H2O/Cu(111) H-up model
H2O/Cu(111) H-down model
O2,sol + H2O/Cu(111) H-up model
O2,sol + H2O/Cu(111) H-down model
O2,ads + H2O/Cu(111) H-up model
O2,ads + H2O/Cu(111) H-down model
−0.10 −0.01 0.00 −0.19 0.22 0.16
0.09 0.01 0.00 0.00 −0.17 0.15
0.00 −0.06 0.05 −0.15 0.27 0.15
0.08 −0.03 −0.08 0.03 −0.18 0.14
−0.07 −0.01 0.01 −0.22 0.64 0.15
0.12 −0.03 0.02 0.11 0.14 0.14
−0.08
−0.08
−0.08
−0.08
−0.09
−0.08
124703-4
A. Sumer and S. Chaudhuri
J. Chem. Phys. 142, 124703 (2015)
FIG. 2. Side and top views of O2 + H2O/Cu(111) configurations: Fully relaxed with O2 located between H2O molecules in Layers II and III, with (a) H-up model (b) H-down model; constrained with O2 kept at approximately 2.2 Å from Cu(111) with (c) H-up model (d) H-down model; and fully relaxed with O2 adsorbed on Cu(111) with (e) H-up model (f) H-down model.
between the solvation enthalpies of O2 in the two H2O/Cu(111) models. It should finally be noted that inclusion of the dispersion effects did not cause a major change in the O2 solvation energy in H2O/Cu(111) or pure water crystal (Table III), and thus, it is concluded that exclusion of these would cause no major effect on the relevant conclusions obtained from rest of the computations. C. Energetics of O2 diffusion in water towards Cu(111) surface
Starting from the coordinates of previously mentioned fully relaxed O2 + H2O/Cu(111) solvation system, we gradually decreased the distance of O2 to Cu(111) and performed constrained geometry optimizations, to find the minimum energy path for diffusion and adsorption of O2 at water/metal interface. In addition, two fully relaxed optimizations were performed in order to identify the stable interstitial positions of O2 between the H2O molecules in Layers I and II and the stable adsorption configuration on Cu(111) surface. Total energies obtained in this way as a function of O2 distance for H-up and H-down models are plotted in Figures 3(a) and 3(b), respectively, where the values are referenced to the initial total
energy of solvation system in order to make it easier to follow the energy trajectory for the O2 molecule approaching the metal surface. For both water models, the plots of energy have two peaks, which correspond to the configurations at which O2 is in the middle of hexagonal H2O rings, i.e., when they are “passing” through the water Layers II and I, in order. The energy barriers that are associated with the first of these peaks correspond to 0.36 and 0.32, for H-up and H-down models, respectively. Thus, the height of the first energy barriers is very close for the two models and in addition, they are close to the height of the diffusion barrier calculated with the pure water crystal (Figure 1(c)), at 0.34 eV, indicating the relatively small influence of water/metal interface on O2-H2O interaction at this distance. It should be noted here that the barriers of diffusion through ordered water layers in our models are higher than the experimentally measured barriers of O2 diffusion in bulk water, at 0.20 eV.37 Unlike the first ones, the second energy barriers, calculated as the energy difference of O2 + H2O/Cu(111) with O2 passing through water Layer I and relaxed at the interstitial position between Layers I and II, as 0.43 and 0.34 eV, for Hup and H-down models, respectively (Figures 3(a) and 3(b)), are found to be strongly dependent on the orientation of H2O with respect to the metal surface. Using these energy barriers
TABLE III. Solvation energy, charge (in electronic units, e, from Hirshfeld scheme) and distance of O2 to Cu(111) for fully relaxed O2 + H2O/Cu(111) with O2 between water layers II and III, with H-up and H-down models (Figures 2(a) and 2(b)) and in pure water crystal (Figure 1(c)). Calculations are done with and without dispersion corrections (TS—see text).
Esolvation/eV Esolvation/eV (with TS) q q (with TS) d/Å d/Å (with TS)
H-up model
H-down model
Water crystal
0.02 0.10 −0.19 −0.19 8.03 8.06
0.48 0.50 0.11 0.11 8.88 8.92
0.38 0.31 0.11 0.11 ... ...
124703-5
A. Sumer and S. Chaudhuri
J. Chem. Phys. 142, 124703 (2015)
become 1.29 × 105 and 4.09 × 106 s−1, respectively. Assuming the O2 concentration at the H2O/Cu(111) interface to be identical to the O2 solubility in bulk water at 298 K and 1 atm, i.e., 8.3 mg/l, the rate of adsorptions for H-up and Hdown models are found as 7.29 × 1018 and 2.30 × 1020 m−2 s−1, respectively. For comparison, the rate of adsorption on Cu(111) in vacuum, as calculated from r ads = √
P
S(θ),
2πmkT
(5)
using the partial pressure of O2 in air at 1 atm and the initial sticking coefficient, S(0), of 0.001,13 is 5.71 × 1023 m−2 s−1, which is higher than the adsorption rate of O2 in solution because of the lack of barriers to the adsorption in vacuum. It should also be noted that as an alternative path to molecular diffusion and then adsorption onto Cu(111), we also tested the possibility of chemical transformation of O2 to OOH in solution by exploring the minimum energy path and the energetics of the reaction.29 We identified an energy barrier of 0.50 eV for reduction, which is higher than the barrier for diffusion and adsorption. The reaction was also found to be endothermic, indicating an overall lack of thermodynamic and kinetic driving force for chemical transformation of O2 to OOH in water. D. Structural and energetic properties of O2 adsorption over Cu(111)
FIG. 3. Relative energies of constrained (•) and fully relaxed () O2 + H2O/Cu(111) configurations as a function of distance of O2 to Cu(111) with (a) H-up and (b) H-down models. Energy barriers of O2 diffusion through H2O molecules in Layer I are also shown.
as activation energies for adsorption, Eact, in an Arrhenius form of rate equation, k = k 0e
−E act kT
,
(2)
it can be shown that there is a difference in the adsorption rate constants, k, between the two models at 298 K, by H −up
H −down) −(E act −E act k H −up 1 kT =e = , (3) H −down 31 k meaning a higher rate of adsorption of O2 molecule on Cu(111) when the water/metal interface is composed of H-down oriented H2O molecules. It is possible to calculate the pre-exponential factor, k0, from initial state f (4) k 0 = transition state f
as 2.62 × 1012 with the vibrational frequencies, f, calculated, for simplicity, in pure water crystal. Then, adsorption rate constant, k, for H-up and H-down models (kH−up and kH−down)
Starting from the outcome of the final constrained optimization of O2 + H2O/Cu(111) with O2 at a distance of approximately 2.2 Å (Figures 2(c) and 2(d)), a fully relaxed optimization is performed during which O2 forms bonds with Cu atoms and gets adsorbed on Cu(111). During the relaxation and bond forming between O2 and Cu atoms, O2 changes its spin state from triplet to singlet. The crossing of the spin states happens at approximately 1.9-2.0 Å distance to the metal surface. It was previously found using DFT calculations that a similar decrease in the spin state of the adsorption complex occurs at Cun clusters with n as small as 5 atoms.15 The final adsorption configuration of O2 over Cu(111) is found to be dependent on the computational model used: It is obtained as three-fold with O-O molecular axis perpendicular to the Cu(111) surface, and three-fold with O-O axis parallel to Cu(111) surface with Hup and H-down models, respectively (Figures 2(e) and 2(f), Table IV). DFT computations on vacuum phase O2/Cu(111) gave the energetically most preferred adsorption configuration of O2 on Cu(111) in absence of H2O molecules as parallel oriented on four-fold site (Table IV), while a parallel oriented threefold adsorption configuration, which was previously identified as the preferred adsorption configuration in vacuum,18 is only 1 meV higher in energy. We calculated the adsorption energies for the final O2 + H2O/Cu(111) adsorption configurations identified with the two models with Eads = EO2,ads+H2O/Cu(111) − EH2O/Cu(111) − EO2,gas.
(6)
Adsorption energies of O2/Cu(111) in vacuum, for comparison, are calculated with
124703-6
A. Sumer and S. Chaudhuri
J. Chem. Phys. 142, 124703 (2015)
TABLE IV. Top views, adsorption sites, energies and charges (in electronic units, e, from Hirshfeld scheme) of O2 for O2 + H2O/Cu(111) and O2/Cu(111) (vacuum phase). Perpendicular O2
Parallel O2
Three-fold (fcc) −0.25 −0.34
Three-fold (hcp) −0.78 −0.42
Three-fold (fcc) −0.14 −0.32
Four-fold −0.79 −0.40
O2 + H2O/Cu(111)
Adsorption site Eads/eV q O2/Cu(111)
Adsorption site Eads/eV q
Eads = EO2/Cu(111) − ECu(111) − EO2,gas.
(7)
Water does not cause a major change in the stability of parallel oriented O2, while it decreases the adsorption energy of perpendicularly oriented one (Table IV). In both water and vacuum, parallel oriented configuration is more stable than perpendicular one, and the difference in the adsorption energies is larger in vacuum than in water. The identification of the adsorption configuration of O2 as perpendicular to the metal surface in the case of H-up model is a result of the interaction between O2 and H2O through hydrogen bonding, the degree of which is dependent on the
H2O orientations with respect to the metal surface. As it was mentioned before, the distance between the topmost layer Cu atoms and water Layer I is smaller at H-up model than at Hdown one. Consequently, when O2 gets adsorbed on Cu(111) in water, H2O’s in close vicinity at the H-up model hinder its full relaxation to an energetically more stable adsorption configuration, whereas the molecule is, relatively, under less influence of hydrogen bonding at the H-down model. The influence of surrounding H2O’s on perpendicularly adsorbed O2 can also be evaluated through a comparative analysis of local density of states (LDOS) graphs. LDOS at upper atom of perpendicular O2 and H atoms of H-up oriented H2O molecules surrounding it have common features between −23 and −20 eV (Figures 4(a) and 4(b)), indicating the presence of hydrogen bonds between the molecules. For comparison, we also optimized the adsorption geometry of perpendicular O2 on Cu(111) in water with the H-down model and also, in vacuum, and found that LDOS at O2 computed with these two last models are very similar to each other and thus, prove the lack of hydrogen bonds between H2O’s and O2 in H-down model (Figures 4(c) and 4(d)). As a final note, we also checked the influence of lower water density in our models on the final adsorption configuration and found no effect of increasing water density on the final conclusions.29 Lower energy of parallel oriented O2 as compared to the perpendicularly oriented one points out to an energetic driving force for reorientation of the molecule even if its initial adsorption configuration after diffusion is perpendicular. So, in a separate work, we computed the energetics of possible reorientation and chemical transformation paths for O2 over Cu(111) in water as modeled with H-up model and found that the energy barrier for reorientation of perpendicular O2 is higher than the energy barrier for the reduction of the molecule to OOH.38 These
FIG. 4. LDOS at (a) O of perpendicularly adsorbed O2 with H-up model; (b) H of H2O molecules neighboring O2 at model in (a); (c) O of perpendicularly adsorbed O2 with H-down model; (d) O of perpendicularly adsorbed O2 in vacuum.
124703-7
A. Sumer and S. Chaudhuri
findings overall indicate that on the Cu(111) surface, along with thermodynamically and kinetically more favorable threefold adsorbed parallel oriented O2 species, it is possible to observe a less stable perpendicularly adsorbed O2, and the reduction reaction may thus proceed via two separate paths starting from these two different adsorption configurations. E. Adsorption induced electronic changes and effects of an electric field
Through an integration of the charge along the cell length (Figures 5(a) and 5(b)), we calculated the dipole moment of H2O/Cu(111) with H-up and H-down models as 6.51 and −11.81 D (corresponding to 0.095 and −0.172 D/Å2), respectively. The difference in the directions of dipole causes an opposite response in the energy of the two models to variations in an external electric field,29 as was also shown previously for other metal surfaces.39–41 It can be speculated that gradual reorientation of the H2O molecules following the variations in electric field would eventually lead to formation of homogeneous H-up and H-down layers at anodic and cathodic potentials, respectively, assuming the hexagonal ordering in the water layers is preserved duringthe potential changes. Our
J. Chem. Phys. 142, 124703 (2015)
computations indicate that this structuring in water layers at the metal interface has an influence on the energetics of O2 diffusion. Thus, it is expected that as the magnitude of the electric field towards the metal surface increases, i.e., as potentials shift to cathodic values, the rate of O2 diffusion and the coverage of O2 on metal surface would also increase. In the presence of an increasing electric field towards the vacuum, i.e., at anodic electrode potentials, the total coverage of O2 on metal surface would decrease and the ratio of perpendicularly adsorbed O2 to parallel adsorbed one would increase. In parallel to the indirect influence of an electric field through water structuring at the metal interface, the charge transfer from the metal to the adsorbate and the consequent changes in the dipole of the system during adsorption also causes the adsorption energetics to respond to the variations in potential of an electrocatalytic system. The charge density comparison, based on Hirshfeld partitioning (Table II) and also plots of charge density difference (between optimized charge density of the system and the sum of the densities of separate atoms that make up system), averaged in x and y dimensions (Figures 5(a) and 5(b)), before and after O2 adsorption indicate an adsorption induced charge redistribution at the water layers, outside region of the water/metal interface where the electric field is directly influential in an electrocatalytic system. Thus, in order to explicitly evaluate the electric field induced changes in O-metal bond formation energetics, we calculated the dipole of the O2/Cu(111) system in vacuum. We found an adsorption induced dipole moment for O2/Cu(111), for both fourfold and perpendicular adsorption configurations, as −1.32 D, which hints a stabilization (destabilization) at cathodic (anodic) potentials. For O2 adsorbed on either Pt nanoclusters or Pt(111) surface, it was shown previously that the energetic response to an external homogeneous electric field (F) varying between ±0.5 V/Å follows a first order Stark equilibrium, ∆E = −µF,
(8)
where µ is the adsorption induced dipole moment without the field.42 For the O2/Cu(111) system, this equilibrium implies a change in adsorption energy of both fourfold and perpendicular configurations by ±0.14 eV in external electric fields of ±0.5 V/Å. Adsorption induced dipole of Cu(111) is higher than the dipole of O2/Pt(111), calculated previously as −0.36 D, due to the charge state of the adsorbate on Pt, which is as low as −0.09.43 Hirshfeld charge analysis of O2/Cu(111), on the other hand, indicates a relatively larger net charge on the adsorbate, ∼ −0.3-0.4 (Table IV), explaining the calculated relatively larger dipole of the system.
IV. CONCLUSIONS
FIG. 5. Plane (xy) averaged charge density difference at H2O/Cu(111) before (dashed line) and after (solid line) O2 adsorption. (a) Perpendicularly adsorbed O2 with H-up model, (b) parallel adsorbed O2 with H-down model.
We computed the energetics of diffusion in water and surface adsorption for O2 molecule at the water/Cu(111) interface. We found that adsorption rate and configuration of O2 and the structure of water network at the metal interface have some unique correlations which can be exploited using electrostatic potentials for favorable ORR kinetics. So-called H-up oriented H2O molecules cause higher barriers for the diffusion
124703-8
A. Sumer and S. Chaudhuri
of O2 towards Cu(111), as compared to opposite oriented Hdown H2O molecules, and thus, H-up ordering reduces the adsorption rate. As the potentials decrease from positive to negative, i.e., from anodic to cathodic values, H-down ordering becomes energetically more stable, causing the H-up oriented H2O molecules flip over. This indicates that O2 adsorption rate on Cu(111) will increase under cathodic conditions. Along with O2 surface coverage, the potential induced structuring in water network also influences the adsorption configuration of O2 on Cu(111). It is shown that while the thermodynamically most stable parallel adsorption configuration is also kinetically favorable at cathodic potentials, a perpendicular adsorption configuration is expected to increase its relative coverage at anodic potentials. Analysis of the local density of states indicated the role of hydrogen bonds in the adsorption of O2 at this configuration. The calculated dipole moment of the O2 on Cu(111) at both configurations was the same and higher than was found previously for O2/Pt(111). The fundamental barrier to efficient ORR on cheaper and more abundant metals such as Cu(111) is primarily due to lack of understanding of the double layer under electrostatic potentials. The complexity of the electrical double layer still needs further investigation such as a realistic polarization scheme for the water layer and ionic species adsorption, e.g., NaCl, on Cu-surface, which can affect the water structure at metal interface.44 Thermal effects on the structural, electronic, and chemical properties of the system would influence the interaction between water and Cu(111) surfaces, so are needed to be included with molecular dynamic simulations. For a water double layer without charged species, the DFT results presented in this work for a planar Cu(111) surface reveal important details on the role of water bilayer in controlling the path and adsorption kinetics of O2. Further spectroscopic measurements such as sum frequency generation spectroscopy will be important for identifying the configuration of the adsorbate in anodic/cathodic potentials identified in this work. ACKNOWLEDGMENTS
The authors acknowledge Joseph Osborne from Boeing Company for his thoughtful suggestions and interest in ORR processes on Cu surface. S.C. acknowledges funding from Joint Center for Aerospace Technology Innovation (JCATI) and Boeing Company. 1P. A. Thiel and T. E. Madey, Surf. Sci. Rep. 7, 211 (1987);
ibid. 46, 1 (2002). Meng, L. F. Xu, E. G. Wang, and S. W. Gao, Phys. Rev. Lett. 89, 176104 (2002). 3M. Morgenstern, J. Muller, T. Michely, and G. Comsa, Z. Phys. Chem. 198, 43 (1997).
2S.
J. Chem. Phys. 142, 124703 (2015) 4K.
Jacobi, K. Bedürftig, Y. Wang, and G. Ertl, Surf. Sci. 472, 9 (2001).
5A. Glebov, A. P. Graham, A. Menzel, and J. P. Toennies, J. Chem. Phys. 106,
9382 (1997). Haq, J. Harnett, and A. Hodgson, Surf. Sci. 505, 171 (2002). 7S. Standop, A. Redinger, M. Morgenstern, T. Michely, and C. Busse, Phys. Rev. B 82, 161412 (2010). 8T. Sueyoshi, T. Sasaki, and Y. Iwasawa, Surf. Sci. 365, 310 (1996). 9K. Prabhakaran, Surf. Sci. 225, L25 (1990). 10M. K. Rajumon, K. Prabhakaran, and C. N. R. Rao, Surf. Sci. 233, L237 (1990). 11G. R. Gruzalski, D. M. Zehner, and J. F. Wendelken, Surf. Sci. 159, 353 (1985). 12A. Spitzer and H. Lüth, Surf. Sci. 118, 136 (1982). 13F. H. P. M. Habraken, C. M. A. M. Mesters, and G. A. Bootsma, Surf. Sci. 97, 264 (1980). 14R. W. Judd, P. Hollins, and J. Pritchard, Surf. Sci. 171, 643 (1986). 15X. Yuan, L. Liu, X. Wang, M. Yang, K. A. Jackson, and J. Jellinek, J. Phys. Chem. A 115, 8705 (2011). 16A. Soon, M. Todorova, B. Delley, and C. Stamp, Phys. Rev. B 73, 165424 (2006). 17L. Padilla-Campos, J. Mol. Struct.: THEOCHEM 815, 63 (2007). 18Y. Xu and M. Mavrikakis, Surf. Sci. 494, 131 (2001). 19M. Alatalo, A. Puisto, H. Pitkänen, A. S. Foster, and K. Laasonen, Surf. Sci. 600, 1574 (2006). 20M. Neurock, R. A. Van santen, W. Biemolt, and A. P. J. Jansen, J. Am. Chem. Soc. 116, 6860 (1994). 21C. D. Taylor, R. G. Kelly, and M. Neurock, J. Electrochem. Soc. 154, F217 (2007). 22B. Delley, J. Chem. Phys. 113, 7756 (2000). 23J. T. Drotar, T.-M. Lu, and G.-C. Wang, J. Appl. Phys. 96, 7071 (2004). 24B. Santra, A. Michaelides, and M. Scheffler, J. Chem. Phys. 127, 184104 (2007). 25F. Li, L. Wang, J. Zhao, J. R. Xie, K. E. Riley, and Z. Chen, Theor. Chem. Acc. 130, 341 (2011). 26J. Carrasco, B. Santra, J. Klimes, and A. Michaelides, Phys. Rev. Lett. 106, 026101 (2011). 27CRC Handbook of Chemistry and Physics, edited by D. R. Lide (CRC Press, Boca Raton, FL, 1993). 28A. Michaelides and K. Morgenstern, Nat. Mater. 6, 597 (2007). 29See supplementary material at http://dx.doi.org/10.1063/1.4914901 for the static and dynamics computations on structure of water, effects of water density and electric field. 30R. Nadler and J. Fernandez Sanz, J. Mol. Model 18, 2433 (2012). 31S. Schnur and A. Groß, New J. Phys. 11, 125003 (2009). 32A. Tkatchenko and M. Scheffler, Phys. Rev. Lett. 102, 073005 (2009). 33J. Carrasco, J. Klimeš, and A. Michaelides, J. Chem. Phys. 138, 024708 (2013). 34R. M. Lynden-Bell and J. Rasaiah, J. Chem. Phys. 107, 1981 (1997). 35Z. Li, “Improving ab initio molecular dynamics of liquid water,” Ph.D. dissertation (Department of Physics, Princeton University, 2012). 36M. H. Abraham, J. Am. Chem. Soc. 104, 2085 (1982). 37C. F. Wong and W. Hayduk, Can. J. Chem. Eng. 68, 849 (1990). 38A. Sumer and S. Chaudhuri, “Modeling aqueous phase electro-reduction of O2 on Cu(111) using density functional based computations” (unpublished). 39C. G. Sanchez, Surf. Sci. 527, 1 (2003). 40A. Roudgar and A. Groß, Chem. Phys. Lett. 409, 157 (2005). 41G. L. Borges, J. G. Gordon, O. R. Melroy, D. G. Wiesler, D. Yee, and L. B. Sorensen, Nature 368, 444 (1994). 42M. P. Hyman and J. W. Medlin, J. Phys. Chem. B 109, 6304 (2005). 43L. Qi, X. Qian, and J. Li, Phys. Rev. Lett. 101, 146101 (2008). 44C. G. Guymon, R. L. Rowley, J. N. Harb, and D. R. Wheeler, Condens. Matter Phys. 8, 335 (2005). 6S.
