Theoretical Analysis of Electrochemical Formation and Phase

Research Article www.acsami.org

Theoretical Analysis of Electrochemical Formation and Phase Transition of Oxygenated Adsorbates on Pt(111) Junxiang Chen, Siwei Luo, Yuwen Liu, and Shengli Chen* Hubei Key Laboratory of Electrochemical Power Sources, Key Laboratory of Analytical Chemistry for Biology and Medicine (Ministry of Education), Department of Chemistry, Wuhan University, Wuhan 430072, China

ACS Appl. Mater. Interfaces 2016.8:20448-20458. Downloaded from pubs.acs.org by UNIV OF GOTHENBURG on 01/25/19. For personal use only.

S Supporting Information *

ABSTRACT: The electrochemical oxygenation processes of Pt(111) surface are investigated by combining density functional theory (DFT) calculations and Monto Carlo (MC) simulations. DFT calculations are performed to construct force-field parameters for computing the energy of (√3 × √3)R30°-structured OH*-H2O* hydrogen-bonding networks (differently dissociated water bilayer) on the Pt(111) surface, with which MC simulations are conducted to probe the reversible H2O* ↔ OH* conversion in OH*-H2O* networks. The simulated isotherm (relation between electrode potential and OH* coverage) agrees well with that predicted by the experimental cyclic voltammetry (CV) in the potential region of 0.55−0.85 V (vs RHE). It is suggested that the butterfly shape of CV in this region is due to different variation trends of Pt-H2O* distance in low and high OH* coverages. DFT calculation results indicate that the oxidative voltammetry in the potential region from 0.85 V to ca. 1.07 V is associated with the dissociation of OH* to O*, which yields surface structures consisting of OH*-H2O* networks and (√3 × √3)-structured O* clusters. The high stability of the half-dissociated water bilayer (OH*-H2O* hydrogen-bonding network with equal OH* and H2O* coverages) formed in the butterfly region makes OH* dissociation initially very difficult in energetics, but become facile once starts due to the destabilization of OH* by the formed O* nearby. This explains the experimentally observed nucleation and growth behavior of O* phase formation and the high asymmetry of oxidation−reduction voltammetry in this potential region. KEYWORDS: oxygenated adsorbates, Pt(111), density functional theory, Monto Carlo simulations, phase transition, nucleation and growth

■

INTRODUCTION Platinum (Pt) is a versatile electrocatalyst for reactions relevant to energy conversion, for example, oxygen reduction, oxygen evolution, and oxidation of organic molecules. The oxygenated adsorbates play crucial roles in these electrocatalytic processes, by acting as reaction intermediates, surface spectators, and/or surface oxidants.1−8 Besides, they are also related to oxidation, reconstruction, and dissolution of surface atoms, which are vital to the stability of Pt electrocatalysts.4,7,9 Therefore, knowledge on the formation and structural evolution of oxygenated adsorbates on the Pt surface is of great significance in mechanistic understanding of various electrocatalytic reactions and in activity and stability optimization of electrocatalysts. Due to the limited accessibility of electrode/electrolyte interfaces to the ultravacuum (UHV) and low-temperature conditions, direct observation of electrochemically formed adsorbates on metal surfaces using modern surface-science techniques has been difficult. Even for the synchrotron-radiation based spectroscopies that could be used at solid/solution interfaces, it remains challenging to distinguish various oxygenated adsorbates on electrode surface due to their similar structures to each other and to the solvation water molecule. Therefore, current insights © 2016 American Chemical Society

into the adsorption and phase structures of oxygenated species on Pt electrode surface have been mainly gained from cyclic voltammetry (CV)1,2,4,10−13 and theoretical simulations.14−22 Especially, CV has been most widely used to probe the surface electrochemical processes of Pt. Figure 1a displays the well-established CV associated with the formation and reduction processes of oxygenated adsorbates on Pt(111) in 0.1 M HClO4 (a none-specifically electrolyte).11 In this study, we are mainly interested in the surface oxidation/reduction processes involved in regions I and II. Bulk oxidation of Pt may take place at potentials beyond region II through the so-called place-exchanging processes.14,23 The voltammetric response in region I is most likely associated with the formation/reduction of adsorbed hydroxyl groups (OH*, the postfix * is used to indicate an adsorbed species in this paper) from/to water molecules (H2O).2,4,10−19,23 It shows a unique butterfly shape, consisting of a pair of sharp wing oxidation/reduction peaks around 0.8 V and a pair of body Received: April 17, 2016 Accepted: July 5, 2016 Published: July 5, 2016 20448

DOI: 10.1021/acsami.6b04545 ACS Appl. Mater. Interfaces 2016, 8, 20448−20458

Research Article

ACS Applied Materials & Interfaces

Figure 1. (a) CV responses of Pt(111) in O2-free 0.1 M HClO4 measured with a potential scanning rate of 2 mV/s (replotted according to data in ref 11); (b) (√3 × √3)R30°-structured H2O*-H2O* network, i.e., the so-called water bilayer (left panel), and H2O*-H2O* network with equal coverage (1/3 monolayer) of OH* and H2O* (right panel) on Pt(111). The water bilayer consists of parallel and vertical H2O*s (H2O*p and H2O*v) respectively. The OH*-H2O* network in the right panel can be derived from the H2O*-H2O* network through dissociation of H2O*vs.

peak,9 which, however, were most likely associated with the defects (e.g., steps) on the prepared Pt(111) electrode surface.1,10,11 In addition, the simulated adsorption isotherms (coverage-potential relation) in these DFT+DMC studies failed to predict the butterfly shape of voltammetry, especially the steep OH* coverage change associated with the sharp wing peak (e.g., Figure 3 in ref 16). The deviation of these simulation results from the experiment voltammetry might be due to the inaccuracy of the kinetic parameters determined by DFT calculations. So far, it remains a great challenge to calculate the kinetic parameters of electrochemical proton/electron transfer reactions using DFT. The DMC simulations by Koper et al.19 have suggested that the butterfly voltammetry is associated with a disorder−order phase transition process of surface adsorption structures, which leads to formation of a (√3 × √3)R30°-structured OH* layer on the Pt(111) surface. In addition, the disorder−order transition model has been shown to reasonably explain the butterfly voltammetry observed for a variety of surface electrochemical processes on single crystal electrodes, which generally involve formation of an ordered adsorption layer.20−22 One of the particular example is the adsorption of Br− on Ag(100), which has been experimentally evidenced by surface X-ray scatting (SXS) experiments to undergo a disorder−order transition processes,29 The MC simulations on Br−/Ag(100) system agreed excellently with the experimental observation.22 These simulations have also revealed that the phase transition process (therefore the butterfly voltammetry) can be modulated through the distance-dependent interaction between adsorbates. For OH* adsorption, the coadsorption of H2O* would bring about some additional contributions, such as OH*H2O* and Pt- H2O* interactions, to the system energy, which have not been properly considered in these simulations. As will be shown in the present study, the Pt- H2O* interactions could dominate the energetics for OH* formation and reduction. The CV responses in region II should be mainly due to the formation and reduction of O*s.2,4,10−14 Recent electrochemical impedance analysis combined with DFT calculations by Bondarenko et al.14 has suggested that the oxidation process in this region yields an ordered phase of 1/3 ML O*s. As indicated by the highly unsymmetrical shape and positions of the oxidation and reduction peaks, the processes involved in this region are very irreversible and may have considerable kinetic limitation. Voltammetric measurements by Feliu and coworkers10,11 have suggested that the surface oxidation process in region II undergoes a nucleation and growth (N&G)

peaks spanning potentials from ca. 0.55 to 0.78 V. Following the sharp wing peak in the anodic branch, there is a wide flat region where only the double-layer charging current can be seen.1,10−12 When the potential is swept back in this region, the reduction branch exhibits a mirrored flat response.10,11 If assuming one electron transfer for per OH*, the voltammetric charge involved in the region I corresponds to the formation of ca. 1/3 monolayer (ML) OH*.1,10,11,16,19,25 Therefore, it has been generally believed that a √3 × √3R30°structured hexagonal coadsorption network, in which OH* and H2O* alternately reside atop Pt atoms and both have a coverage of 1/3 ML (Figure 1b, right panel), forms at the end of region I.10,11,14,16,18,25 This particular OH*-H2O* network, which will be denoted as the half-dissociated water bilayer in the following, can be derived through the dissociation of the vertical water molecules (H2O*v) in the H2O*-H2O* network of the same hexagonal structure, which will be denoted as the nondissociated water bilayer in the following (Figure 1b, left panel). Due to the simultaneous presence of extensive hydrogen bonds and Pt−OH bonds, the formation of the half-dissociated water bilayer should be energetically very favorable.15,24−27 The high stability of the half-dissociated water bilayer against further oxidation should be the origin of the wide flat region following the sharp wing peak in the anodic branch. Berna et al.28 have attributed the body and wing regions in the butterfly voltammetry of Pt(111) electrode to the dissociation of two types of water molecules, i.e., the free water molecules in the solution and the solvation water molecules bonded with anions, which were assumed to have different reaction energies. The experimental voltammetry thus has been numerically fitted based on empirical adsorption isotherms with adjustable parameters. Viswanathan et al.16 and Rai et al.18 have investigated the potential-dependent OH* and O* formation on Pt(111) with dynamic Monte Carlo (DMC) simulations based on kinetic and thermodynamic parameters computed by density functional theory (DFT) calculations. The results predicted the onset potential of OH* formation in good agreement with the experiment voltammetry. It was also shown that OH* adsorption is accompanied by H2O*, and they form (√3 × √3)R30°-structured network. However, these simulations suggested an onset potential lower than 0.9 V (vs RHE) for OH* dissociation, which seemed to conflict with the nearly vanished faraday current following the sharp wing peak as shown in Figure 1a. In some earlier studies, the measured CVs did show slight oxidation current after the sharp wing 20449


Research Article


Figure 2. Histograms of ni ≈ εi (red), ni′-εi (black) and ni″-εi (blue) spectra obtained through normal metropolis sampling, modified metropolis sampling with α = 0 and modified metropolis sampling with α = 0.5 respectively for θOH* = 0.25 ML (a) and θOH* = 0.17 ML (b). Green dash lines indicate the overlapping regions used to calculate λ (see the main text for details).

reaction, reaction R1 produces faraday currents. As implied by the highly symmetric oxidation and reduction responses and the nearly invariant shape and position of voltammetric peaks with potential scanning rates,10,11 reaction R1 is highly reversible. Therefore, it can be investigated with equilibrium MC simulations. Under equilibrium approximation of reaction R1, one has the following U ≈ θOH relation

mechanism, which may explain why the region II starts with a relatively sharp rise of oxidation current after a wide region in which nearly no oxidation current occurs, and why the oxidation and reduction voltammetry in this region is highly asymmetric. However, atom-level evidence has been rarely presented so far to reveal the nature of such a N&G-like behavior. Besides, it has been also disputed that whether only OH*s or both of OH*s and H2O*s are dissociated in this region.10 In this paper, we attempt to gain atom-level insights into the electrochemical oxygenation processes of Pt(111) surface through combined DFT calculations and Monto Carlo (MC) simulations. Considering the reversibility of the butterfly voltammetry and the possible formation of the half-dissociated water bilayer structure, the equilibrium MC simulations are used to investigate the conversion between H2O* and OH* in water bilayer networks, with the DFT-calculated energies of various water bilayer structures being used as inputs. The results produce a butterfly voltammetry similar to that observed experimentally, and suggest that the distance change between the water bilayer and Pt along with H2O* ↔ OH* conversion is responsible for the butterfly shape of the voltammetric responses. DFT calculations results indicate that the oxidative voltammetry in region II of CV is associated with the dissociation of OH* to O*, which yields surface structures consisting of OH*-H2O* networks and (√3 × √3)-structured O* clusters. The calculations also explain the experimentally observed N&G behavior of O* phase formation and the high asymmetry of oxidation−reduction voltammetry in this potential region.

■

eU = dE /d(NθOH) − T dS /d(NθOH) + 1/2μH2

where e is the electron charge, and μH2 is the standard chemical potential of the gaseous H2 molecule, which can be obtained from the DFT-calculated energy, vibrational entropy, and zeropoint energy (ZPE) of the H2 molecule. Here, the 1/2 μH2 is used to substitute for the chemical potential of H+ in solution (μH+) since they are the same in value as the potential is referenced to that of the reversible hydrogen electrode (RHE).15,30 The first term in the right side of eq 1, namely, the derivative of system energy (E) against the number of OH* at the surface (NθOH), represents the energy change due to the conversion of an H2O* to OH* in water bilayers at OH* coverage of θOH, which can be considered the standard adsorption energy of OH* at θOH. The term of dS/d(NθOH) in the right side of eq 1 is the change of configurational entropy (S) due to the conversion of an H2O* to OH* at θOH. The ways to determine E and S through MC simulations at a given θOH are described in the following. Under a given θOH, the OH*-H2O* network could have a variety of different configurations, with total number W = (N/ 3)!/((N/3-NθOH)!(NθOH)!). Each configuration represents a microscopic state of the system. Within a macroscopic time, the system undergoes various configurations with different probabilities, which can constitute a canonical ensemble. The thermodynamic mean energy of the canonical ensemble would give the macroscopic energy of the system, E, that is,

METHODS AND COMPUTATIONS

Monte Carlo Simulations. The canonical ensemble-based equilibrium Monte Carlo simulations were used here to investigate the conversion between H2O*v and OH* in the (√3 × √3)R30°-structured OH*-H2O* networks, namely,

E=

∑ pw εw = ∑ pεi εi w

H 2O*N(2/3 − θOH)OH*NθOH ↔ H 2O*N(2/3 − θOH) − 1OH*NθOH + 1 + H+ + e−

(1)

(2)

i

where pw refers to the probability that the system has the configuration w, with which the system has an energy of εw; pεi refers to the probability that the system has the energy of εi. According to the Boltzmann distribution law

(R1)

where N and θOH are the number of surface Pt atoms and OH* coverage. The surface structure would be the nondissociated water bilayer as θOH = 0 and the half-dissociated water bilayer as θOH = 1/3 ML. As a coupled proton/electron transfer

pεi = ni /Ñ = gi exp( −εi /kT )/∑ gi exp( −εi /kT ) i

20450

(3)


Research Article

ACS Applied Materials & Interfaces where ni is the number of configurations with energy εi in the assemble, Ñ = ∑ni is the total configuration number in the assemble, and gi is the degeneracy of the system at the energy level εi, i.e., the number of possible configurations (microscopic states) with energy εi; therefore, ∑gi = W. The values of ni depend on Ñ . Combining eqs 3 and 2 gives

E=

∑ niεi /Ñ

At medium OH* coverages at which the ni ≈ εi, and ni′ ≈ εi spectra have very weak overlapping, a sampling with α = 0.5 has to be performed. The resulted ni″ ≈ εi spectrum (Figure 2b, blue histograms) overlaps with both ni ≈ εi, and ni′ ≈ εi spectra, and λ can be calculated according to (the derivation is similar to that for eq 8) λ=

(4)

i

To determine the value of E for a given θOH, we used the Metropolis sampling algorithm to construct the assemble. The details are given in Supporting Information (SI, section S1). The simulations were performed on a 20√3 × 20√3-sized Pt (111) surface (N = 1200) by using a homemade code. The red histograms in Figure 2 show the ni ≈ εi distribution for two typical θOH* obtained through Metropolis MC sampling with Ñ = 106. In the Metropolis sampling process, the sampling probability of configurations well follows the Boltzmann distribution law (eq 3) as Ñ is large enough.31,32 Therefore, the obtained ni ≈ εi distributions can be used to calculate the system energy according to eq 4. For canonical ensemble, the entropy of the system can be determined with S = k ln ∑ gi exp[(E − εi)/kT ]

∑i ni″ exp(0.5εi /kT ) W ′ ∑j nj″ exp(0.5εj /kT ) ∑i ni′ W

NθOH

εw = E(N , θOH , w) − E(N , 0) =

To obtain the entropy, one has to know the gi-εi spectra of the system. According to eq 3, we have ⎛ε ⎞ gi = ni exp⎜ i ⎟ /λ ⎝ kT ⎠

(7)

Depending on the weighting factor α (0 ≤ α ≤ 1), the sampled assemble would be dominated by configurations spanning different energy ranges. The normal metropolis sampling (α = 1) gives ni-εi spectra located in the lowest energy region (Figure 2, red histograms). When α = 0, various configurations will be sampled with the same probability. This uniform sampling yields ni′-εi spectra located in relatively high energy region (Figure 2, black histograms), with the population distributions of configurations being very similar to that in gi ≈ εi spectrum. That is, ni′/W′ = gi/W, where W′ is the total number of the sampled configurations (here it is 106). In the cases when the ni ≈ εi, and ni′ ≈ εi spectra have a considerable overlapping region (Figure 2a), which occurs at relatively low and high OH* coverages (θOH* < 0.05or θOH* > 0.25 ML), The values of λ can be calculated according to (section S2 in SI for derivation details). λ = (W ′ ∑ ni exp(εi /kT ))/(W ∑ ni′) i

i EOH *

(10)

where E(N,θOH,w) refers to the energy for the studying system; E(N,0) refers to the energy of the system consisting of a nondissociated water bilayer (θOH = 0) on the same Pt(111) slab; EiOH* represents the formation energy of the ith OH*; and i ranges over all OH*s in the H2O*-OH* network. Through this definition, we have assumed that the difference between E(N,θOH,w) and E(N,0) can be expanded into the sum of the energy change due to the formation of all the individual OH*s in the network, i.e., EiOH*. For convenience, EiOH* will be denoted as the energy of the ith OH*. The values of EiOH* should depend on the neighboring environments of the corresponding OH*, namely, the numbers and distributions of the nearby OH*s. According to the numbers (k) and distribution modes (m) of the first nearest OH* neighbors, we can distinguish the surface OH*s into 13 types (SI, section 3: Figure S1), the energy of which correspond to 13 force-filed parameters (Ek,m). We also considered the effect of the second nearest OH* neighbors on the EiOH* by introducing a correcting term El to Ek,m, which is a summation of the contributions from 6 types of second nearest OH* neighbor environments (SI, section 4: eq S1 and Figure S2). These second-nearest-neighbor contributions correspond to 6 additional parameters. The definition of the first and second nearest sites is also illustrated in Figure 3. The energy of a given OH* can then be expressed by

where λ is a constant (SI, section S2). Thus, the gi-εi spectrum for a given θOH can be derived from the corresponding ni-εi spectrum obtained with Metropolis MC sampling if λ is known. In order to get λ, we also performed MC simulations by modifying the Metropolis sampling algorithm (SI, section S1), so that the appearing probability of configurations with energy εi is i

∑ i=1

(6)

pεi ′ = gi exp( −αεi /kT )/∑ gi exp( −αεi /kT )

(9)

in which j ranges over the overlapping region between ni ≈ εi and ni″ ≈ εi spectra (region A in Figure 2b) and i ranges over the overlapping region between ni′ ≈ εi and ni″ ≈ εi spectra (region B in Figure 2b). Force-Field Parameterization of System Energy. To calculate the macroscopic mean energy E and configurational entropy of the system at θOH using eqs 3 and 4, it is necessary to know the system energy at various configurations of an OH*-H2O* network, namely εw (see eq 2). For systems with very large sizes used in MC simulations (here N = 1200), direct calculation of εw using first-principle method is impossible. Therefore, we used a set of force-field parameters derived from DFT calculations on Pt(111) slabs with relatively smaller sizes to calculate the values of εw of large slabs. For simplicity, the value of εw is given as

(5)

i

i

∑j nj exp(εj /kT )

i EOH * = Ek , m + E l

(11)

Using DFT calculations, we obtained εw values of 39 systems which consist of Pt slabs with relatively small sizes (N = 9−27) and OH*-H2O* networks with different θOH and different configurations w. These DFT-calculated εw values are linear expressions of the 19 force field parameters, which yield a set of 39 linear equations containing 19 unknown variables. By solving these equations with the least-square method, the values

(8)

in which i ranges over the overlapping region (Figure 2a). 20451


Research Article


Figure 3. Illustrations of first (left panel) and second (right panel) nearest OH* neighbors (yellow) for a considered OH* (green). Figure 4. Comparison of the U ≈ θOH relation obtained from the MC simulations and that predicted by CV.

of the force field parameters can be determined (SI, section 5). These force-field parameters are then used to calculate εw in MC simulations according to eqs 10 and 11. DFT Calculations. Spin-polarized DFT calculations were performed with periodic supercells under the generalized gradient approximation (GGA) using the Perdew−Burke− Ernzerhof (PBE) functional for exchange-correlation and the ultrasoft pseudopotentials for nuclei and core electrons. The Kohn−Sham orbitals were expanded in a plane-wave basis set with a kinetic energy cutoff of 30 Ry and the charge-density cutoff of 300 Ry. The Fermi-surface effects have been treated by the smearing technique of Methfessel and Paxton, using a smearing parameter of 0.02 Ry. Periodically repeated four-layer slabs of face centered cube (fcc) structure were used to model the Pt(111) electrodes, with the bottom two layers fixed to DFT-optimized equilibrium lattice constant of 4.00 Å and the top two layers were allowed to relax with the adsorbates during the calculations until the Cartesian force components acting on each atom were below 10−3 Ry/Bohr and the total energy converged to within 10−5 Ry. Slabs of a variety of sizes, ranging from (3 × 3) to (3√3 × 3√3) were used. For slabs of (2√3 × 2√3) or smaller sizes, the Brillouin-zones were sampled with a 5 × 5 × 1 k-point mesh, while a 4 × 4 × 1 k-point mesh was used for those of larger sizes. The PWSCF codes contained in the Quantum ESPRESSO distribution33 were used to implement most of the calculations. The molecular entropies and zero-point-energies (ZPE) of H2O and H2, and the ZPE of O* were included in calculating the free energy of various surface reactions. The data are listed Table S2 (SI, section S5), in which the ZPE values were calculated using the Phonon module in ESPRESSO and the molecular entropies were taken from ref 34. The molecular entropies of adsorbed species were neglected. A temperature (T) of 300 K was used in the free energy calculations.

curve) given in the inset of Figure 4, which is related to the faraday current according to the following equation j = dQ /dt = d(eNθOH))/d(U /s) = seN (dθOH/dU ) (12)

Where s stands for the scanning rate. As compared with the experimental voltammetry, the simulated wing peak appears at some lower potential and is a bit broadened. This deviation might be due to the inaccuracy of the DFT-based force-field parameters. In the present DFT calculations, we have neglected some complexities that could affect the system energy change associated with the conversion between H2O* and OH*, for example, the changes in the zero-point-energy and charge density on the surface. These factors should be θOH-dependent. Better agreement between the simulated and experimental voltammetry could be achieved by introducing modification parameters to account for these factors. However, the results shown in Figure 4 should be able to reasonably assign the butterfly voltammetry to the reversible conversion between OH* and H2O* in the (√3 × √3)R30°-structured OH*H2O* networks. According to eq 12, the current at a given potential depends on the variation rates of θOH with U (or U with θOH). As known from eq 1, the U ≈ θOH relation is determined by the variation of the system energy E and entropy S with θOH. To gain deeper insight into the formation mechanism of the butterfly voltammetry through OH* ↔ H2O* conversion, the variation of the system energy E and entropy S with θOH should be investigated. Figure 5a shows the simulated variation of S with θOH (eq 5) and that predicted by the mean-filed approximation (MFA). The MFA assumes equal accessibility (probability) of different configurations of an OH*-H2O* network, regardless of the energy difference. In this case, S = k ln W = k ln ∑gi. In reality, the appearing probability of a configuration should depend on the system energy under this configuration according to the Boltzmann distribution law (eq 3). One can see that the actual configurational entropy deviates from the MFA prediction at coverages from ca. 0.05 ML to ca. 0.25 ML; while they are nearly the same at lower and higher coverages (Figure 5a). The reasons for this will be discussed later on. Despite the difference in the absolute entropy values, the S ≈ θOH curves given by the MC simulations and MFA have very similar shapes. In results, the derivative of the MFA and MC entropies have nearly the same dependence on θOH (Figure 5b), showing a monotonous (approximately linear) increase in θOH region of about 0.03−0.3

■

RESULTS AND DISCUSSION Conversion between OH* and H2O* in OH*-H2O* Networks. As stated in Monte Carlo simulations, reaction R1 may be responsible for the butterfly voltammetry in CV. To verify this, we compare the theoretical U ≈ θOH curve derived through MC simulations according to eq 1 and that obtained by integrating the butterfly region of CV (Figure 4). It can be seen that the MC results reasonably agree with the experiment. Both curves exhibit a turning of the variation rate of U with θOH (or θOH with U) at coverages around 0.2 ML (or potential around 0.75−0.78 V), which should correspond to the body to wing transition in the butterfly voltammetry. This is more clearly manifested by the derivative of theoretical isotherm (θOH ≈ U 20452


Research Article


Figure 5. (a) Entropies and (b) entropy derivatives obtained from MC simulations and MFA, and (c) energy derivatives from MC simulations. Under MFA, −dS/d(NθOH) = kB ln(θOH/(1/3 − θOH)).

Figure 6. (a) Top panels: top-view illustration of OH*s on Pt(111) obtained by removing all the H2O* molecules in the OH*-H2O* networks with different θOH; bottom panels: side-view illustration of distance variation between Pt and water bilayer. (b) DFT-calculated EOH (red dots) as a function of θOH (bottom x-axis) using models shown in the upper panels of Figure 6a, and EPt−H2O (black dots) as a function of distance between Pt and nondissociated water bilayer (top x-axis). (c) Average distance between Pt and H2O*s in OH*-H2O* networks as a function of θOH.

change the numbers of hydrogen bonds in the OH*-H2O* networks. Thus, we may conclude that the dissociation of a H2O* into OH* should result in negligible change in the interactions within the OH*-H2O* networks. There is another factor that may modulating the system energy E accompanying the OH* ↔ H2O* conversion, namely, the interaction between the water bilayer and Pt surface, which depends on their separation. For nondissociated water bilayer on Pt(111) surface, the optimized Pt−H2O* distance is ca. 3.29 Å, which is much longer than the optimized Pt−OH* distance (ca. 2.12 Å). Thus, the conversion of an H2O* to OH* would drag the adjacent H2O* molecules downward to maintain the hydrogen-bonding network, which may result in repulsion between Pt and water bilayer. To verify this, we estimated the Pt−H2O* repulsion energy by calculating the system energy of nondissociated water bilayer on Pt(111) surface (EPt−H2O, see section S7 in SI) under different Pt−H2O* distances (dPt−H2O, Figure 6a, bottom panels). As seen in Figure 6b (black data points), the calculated EPt−H2O increases rapidly as dPt−H2O is below 3 Å, indicating increased repulsion interaction between Pt and water bilayer with decreasing the dPt−H2O. In Figure 6c, the average Pt−H2O* distances in OH*-H2O* networks are plotted as a function of θOH. One can clearly see a turning at θOH ≈ 0.2 ML, below which the Pt-H2O* distance exhibits a continuous decrease with increasing θOH, while beyond which the Pt-H2O* distance varies little with θOH because that H2O*s are mostly in their appropriate positions to maintain the OH*H2O* network. Thus, we may conclude that the butterfly voltammetry in the region I of CV is mainly due to the variation of Pt-H2O* distance accompanying the conversion between H2O* and OH* in the OH*-H2O* networks. Although the change in the Pt-H2O distance does not directly contribute to

ML, with no turning point that may correspond to the bodywing transition in the butterfly voltammetry. As shown in Figure 5c, the dE/d(NθOH) ≈ θOH plot obtained from MC simulations shows an obvious turning at θOH ≈ 0.2 ML, below which the dE/d(NθOH) continuously increase with θOH, while above which the dE/d(NθOH) slightly decreases with increasing the θOH. The value of dE/d(NθOH) can be considered the standard adsorption energy of OH*. Usually, the standard adsorption energy of an adsorbate increases (goes positively) or decreases (goes negatively) with coverage if there is a repellent or attractive interaction between the adsorbate molecules. If the interaction between adsorbate molecules is weak (Langmuir adsorption), the standard adsorption energy would vary little with coverages. The variation of dE/d(NθOH) shown in Figure 5c thus seemingly suggests a transition from repellent interaction between OH*s to a weak attractive interaction as coverage goes above ca. 0.2 ML. This is somewhat unrealistic since one would expect more pronounced repelling between adsorbates at higher coverages. Considering that the OH*s in the OH*-H2O* network do not immediately neighbor with each other (separated by H2O*, Figure 3), the interaction between OH*s should be weak. To verify this, we calculated the integrated adsorption energy of OH* (EOH, section S7 in the SI) by removing all H2O* molecules in OH*-H2O* networks with different θOH (Figure 6a, top panels). As shown in Figure 6b (red data points), the calculated EOH values vary little with θOH, which suggested that OH*s interact rather weakly in the separated arrangement in OH*-H2O* networks. Similarly, one would not expect significant change in the interaction between H2O*s upon conversion of an H2O* into OH* in the OH*-H2O* networks. In addition, the conversion of a H2O* into OH* does not 20453


Research Article

ACS Applied Materials & Interfaces current production, it affects the relation of system energy E with θOH, and the U ≈ θOH relation accordingly. The later determine the current according to eq 12. It should be mentioned that there have been arguments on whether water molecules are arranged into well-defined bilayer structures on metal surfaces at room temperature,35 although the nondissociated water bilayer on Pt(111) has been implicated by results of the low-temperature spectroscopic measurements in UHV24,36 and DFT calculations.26 Considering the high stability of the OH*-H2O* hexagon network structure, it is expected that the OH*s and H2O*s would form extended networks as the OH* coverage reaches a critical value; even at beginning of the OH* formation (low OH* coverages), at least locally ordered OH*-H2O* hexagon networks would form on the surface. Although the surface sites away from the local OH*-H2O* networks may not have a well-defined water bilayer structure at very low OH* coverages, the calculated dE/d(NθOH) ≈ θOH curve shown in Figure 5c, and accordingly the adsorption isotherm in Figure 4 should be insignificantly affected because the influence of the background water structures on the calculation of E could be largely eliminated through the differentiation. It is also worth mentioning that the MC simulations by Koper et al.19 have shown that the body to wing transition could also be reasonably explained in terms of the disorder− order phase transition of OH* adsorption layer, which results in a large entropy decrease. In these simulations, the disorder− order phase transition crucially depends on the interaction between OH*, which was modeled by summation of distancedependent two-body pair potentials, with the simultaneous adsorption at the immediate nearest neighbor sites being excluded. The possible interaction between OH* and H2O* and that between Pt and H2O* are ignored in these models. Based on the present calculation results (Figure 6b), the variation of OH*-OH* interaction might be less pronounced than that of Pt-H2O* interaction accompanying OH* formation. It should be emphasized that the present results do not deny the order−disorder transition mechanism for the butterfly voltammetry; rather, they offer an alternate explanation. Considering that surface OH*s and H2O*s could interact with solvent molecules in the electric-double-layer, the OH*H2O* network may not fully represent the solvation structures at electrode surface. To investigate the extended solvation effect on the theoretical isotherms, we calculated the solvation energies of various OH*-H2O* networks using the continuum solvation model, namely the COSMO,37−39 in which the solvent environment out of the studying system (the slab and the OH*-H2O* network) is treated as a dielectric continuum of permittivity. We used a dielectric constant of 78.54 for H2O. The calculated solvation energies for the nondissociated and half-dissociated water bilayers were found to differ only ca. 0.008 eV, which is negligibly small as compared with the potential range (ca. 0.55−0.83 V) within which the OH* coverage changes from 0 to 1/3 ML. Now we explain why the configurational entropies given by MC are nearly the same as that predicted by the MFA at low and high θOH but deviate at the medium θOH. The reason can be found in Figure 7, in which the gie(E‑εi)/kT ≈ εi spectra and the corresponding gi ≈ εi spectra are compared for a low, medium and high OH* coverage, respectively. The distribution of gie(E‑εi)/kT determines the actual entropy of the system (eq 5); while the distribution of gi gives the MFA entropy. At a low

Figure 7. Spectra of gi ≈ εi (black) and gie(E‑εi)/kT ≈ εi (red) for different θOH of (a) 0.01, (b) 0.17, and (c) 0.32 ML. The values of gie(E‑εi)/kT in (a−c) have been multiplied by a factor of 103, 1012, and 10, respectively.

or high OH* coverage (Figures 6a and 6c), the system has a relatively small number of possible configurations with which the system energy does not differ significantly; accordingly, the gie(E‑εi)/kT ≈ εi and gi ≈ εi spectra are distributed in the similar narrow energy ranges. In this case, the entropy given by the gi ∼ εi spectra (MFA) would deviate insignificantly from that derived from the gie(E‑εi)/kT ≈ εi spectra (eq 5). At a medium coverage, the system generally has a large number of possible configurations which give a wide range of system energies; therefore, the major energy ranges where the gie(E‑εi)/kT ≈ εi and gi ≈ εi spectra span are fairly apart, making the two spectra give different entropy values. Transition between the OH*-H2O* Network Phase and O* Phase. The further oxidation of the half-dissociated water bilayer may be initialized through the dissociation of an H2O* (reaction R2) or an OH* (reaction R3). Because O* prefers to occupy the 3-fold fcc hollow site while OH* occupies the atop site, the dissociation of OH* to O* results in a vacancy in the hydrogen bonding network, which should be filled by a water molecule from the solution (H2Ol) to maintain the hydrogen bonding network (reaction R3) since an incomplete hydrogen bonding network is energetically much less stable than a complete one. Our calculation shows ca. 0.4 eV increase in the reaction free energy when keeping the vacancy unfilled by the H2Ol. H 2O*1/3N OH*1/3N ↔ H 2O*1/3N − 1OH*1/3N + 1 + H+ + e−

(R2)

H 2O*1/3N OH*1/3N + H 2Ol ↔ H 2O*1/3N + 1OH*1/3N − 1O* + H+ + e−

(R3)

To investigate preferred oxidation mode of the halfdissociated water bilayer, we calculated the standard reaction free energies of the reactions R2 and R3 (ΔG0R2 and ΔG0R3) using DFT with eq S5 and S6 (SI, section S8). In order to minimize the effect of the coverage change in DFT calculations, 20454


Research Article


Figure 8. Standard reaction free energies of reactions (a) R2 and (b) R3 (at 0 V) calculated using slabs of different sizes.

Figure 9. Optimized structures for (a) H2O*1/3N+1OH*1/3N‑1O* formed through reaction R3, (b) H2O*1/3N+2OH*1/3N‑2O*2 in which the second O* forms through dissociation of an OH* located in close proximity to the first O*, (c) H2O*1/3N+2OH*1/3N‑2O*2 in which the second O* forms through dissociation of an OH* located apart from the first O*, and (d) O* clusters and OH*-H2O* networks which coexist on the surface (the dashed lines mark the boundary). In the region of OH*-H2O* network phases, H2O*s reside atop Pt atoms; while H2O*s are located at the fcc sites in the region of O* cluster phase, with those above O* having a H-down configuration. The distance between the adsorbed oxygen atom and water molecule above it is ca. 2.77 Å, closed to the value of O···O distance for hydrogen bonding (2.75 Å). The distortion of water bilayer hexagons can be seen in boundary regions. The adsorbed oxygen atoms (O*) and that in H2O*s are given in yellow and red respectively for better illustration.

potential of region II in CV. This seemingly also confirms the R3 pathway for the further oxidation of the half-dissociated water bilayer. Equation 13 represents a MFA treatment of O* formation. Since this process is considerably irreversible, we do not attempt to use equilibrium MC simulation to treat it. Ideally, the kinetic MC (KMC) simulation should be used. This requires systematic determination of kinetic parameters (e.g., rate constants) for conversion between OH* and O* at different configurations, which is a nontrivial and timeconsuming task. Moreover, DFT calculations on the kinetics of proton-couple electron transfer reactions at electrochemical interfaces remain a great challenge. We are performing KMC simulations on the OH* ↔ O* transformation and hope to get meaningful results in the future. Therefore, we use the MFA to gain some qualitative insights into the O* formation process in present study. Fortunately, as has been implied by above MC simulations results, the MFA should cause negligible deviation of system entropy from the real values in the low or high ranges of adsorbate coverage (Figure 5a). Even though the MFA may give deviated configuration entropy at medium coverages, the differential configuration entropies, which directly determines the free energy of surface reactions, are negligibly affected by using the MFA. In addition, the relatively narrow distribution of system energy with configurations (Figure 7c) should also allow us to use DFT calculations to estimate the system energy at low and high coverages. Thus, the use of MFA should be justified at least at low and high ranges of O* coverages. Because the formation of an O* is always coupled with the conversion of an OH* to H2O* in the OH*-H2O* network (R3), the entropy

slabs as large as possible should be used to calculate the values of ΔG0R2 and ΔG0R3. We have used slabs of different sizes to perform calculations and found that the calculated values converged as N was larger than 21 (Figure 8). Since the values of ΔG0R2 and ΔG0R3 have the same dependence on the electrode potential (see eq S5 and S6), we compare them at zero potential (U = 0 V) for convenience. It can be seen that reaction R3 is ca. 0.2 eV more negative in reaction free energy than reaction R2, which suggests that the further oxidation of the half-dissociated water bilayer proceeds preferably through reaction R3. This is easy to be understood when considering that reaction R2 breaks the hydrogen bonding network by producing neighboring OH*s, which should be energetically disfavored. We can use the converged value of ΔG0R3 obtained at large slabs (ca. 1.03 eV, Figure 8b) to estimate the onset potential for the OH* dissociation in the half-dissociated water bilayer according to eq 13 by using an O* coverage value (θO) that is very small but could be visible in voltammetric responses. eU = ΔG 0 R3 + kBT ln(θO/(1/3 − θO))

(13)

By setting θO = 5 × 10−3, eq 13 gives a thermodynamic onset potential of ca. 0.93 V. Considering the possible kinetic barrier, this estimated onset potential for R3 agrees reasonably with the potentials where the region II starts in the oxidation direction of the CV (Figure 1a). Similarly, one can use the converged ΔG0R2 (ca. 1.24 eV, Figure 8a) to estimate the thermodynamic onset potential for reaction R2. The resulting values would be larger than 1.1 V, which is considerably higher than the onset 20455


Research Article

ACS Applied Materials & Interfaces change associated with an O* formation at relatively low θO would be similar to that for the reversed reaction of R1 at θOH close to 1/3 ML (θO = 1/3-θOH). Thus, the differential configuration entropy associated with an O* formation should be

energy barrier; whereas the formation of additional O* atoms around the already formed one requires lower energy barriers due to the destabilization of OH*s by the nearby O*. The present calculations also point out that the formation of O* phase is different from the usual N&G processes, in which a stable nucleus with a critical size (containing a critical number of atoms) has to form prior to the phase growth. Here, formation of one O* atom is enough to facilitate the phase growth, that is, an O* atom serves as a nucleus for phase growth. The N&G-like formation of O* phase is also indicated by the U-θO spectrum calculated using eq 15, in which the dE/d(NθO) is the differential of the system energy (E) at θO, kB ln(θO/(1/3 − θO)) represents the differential of configurational entropy, ZPEO* refers to the ZPE of O*, and μH2O is the chemical potential of water molecule under the partial pressure of 0.035 atm (the saturated vapor pressure of water at room temperature). It can be seen that the calculated O* formation potentials in a considerable range of O* coverages (θO < 0.15 ML) are less positive than the onset potential of OH* dissociation. The N&G-like formation of O* phase may be the origin of the high asymmetry (irreversibility) between the oxidation/reduction responses in the region II of CV. There might be other factors that contribute to the asymmetry between the oxidation and reduction responses in this region, for example, the chemical state change of surface Pt induced by O* formation. Surface Pt bonded with O* may be in higher valence state than that bonded with OH*. This makes Pt−O* bond stronger than Pt−OH* bond. Therefore, the formation and reduction of surface O* are different in energetics. As seen from Figure 1a, the reduction peak in the region II terminates at ca. 0.85 V, which indicates that the O* atoms are stable at potentials as low as ca. 0.85 V once formed. Interestingly, Figure 10 suggests the similar lowest potential

∂S /∂(NθO) = −∂S /∂(NθOH) = kB ln(θOH/(1/3 − θOH)) (14)

Following reaction R3, the formation of the second O* in a surface area may occur through the dissociation of an OH* located in close proximity to the existing O* (Figure 9, panel b), or an OH* fairly apart from the existing one (Figure 9, panel c). By comparing the energies of systems with surface structures of H2O*1/3N+2OH*1/3N‑2O*2 in which the two O*s are located in different places, we find that the second O* preferably forms through the dissociation of an OH* close to the existing O*. On slab larger with N ≥ 21, the formation of structure in panel b is about 0.14 eV more negative in reaction free energy than that in panel c. This suggests that the continuous dissociation of OH* would result in O* clusters rather than highly dispersed O* atoms on Pt (111) surface, yielding surface structures with O* clusters and OH*-H2O* networks coexisting (Figure 9, panel d). Because an O* formation is accompanied by transformation of an OH* to H2O*, the hydrogen bonding network in O* cluster areas would transform to H2O*-H2O* networks and shift positions so as to form hydrogen bonds with the underneath O*s. In surface regions covered by OH*-H2O* networks, OH*s and H2O*s remain located at top sites. Accordingly, the hexagon structure would be slightly deformed at the boundary between the O* clusters and OH*-H2O* networks (Figure 9, panel d). As O* coverage reaches 1/3 ML, the surface is purely covered by O* and there will form a global water bilayer with all H2O*s residing above O* (Figure S3). Another noticeable feature we have found is that the dissociation of OH*s nearby a formed O* becomes energetically easier than the initial OH* dissociation in the pristine halfdissociated water bilayer. For example, the formation of the second O* near the already formed O* shows a standard reaction free energy that is ca. 0.2 eV more negative than that for the first O* formation (reaction R3) on a slab with N = 21. This further implies that O* would prefer to form in close proximity to each other. It also suggests that the growth of an O* cluster is energetically less demanding than the initialization of O* formation. Considering that the growth of O* cluster and the first O* formation both occur through the OH* dissociation, the lower reaction free energy would correspond to lower active free energy according to the Brønsted−Evans− Polanyi (BEP) principle.40−43 Therefore, one can expect that the growth of O* phase would have lower free energy barrier, therefore more feasible kinetics, than the initial O* formation. Recently, Feliu and co-workers10,11 have concluded from systematic potentiostatic and potentiodynamic measurements that the voltammetric response in the early stage of region II is associated with a nucleation and growth (N&G) process. For instance, they observed chronoamperometric (current−time) curves exhibiting a maximum in a series potential-step experiments,10 which is the characteristic behavior of electrochemical phase formation with N&G mechanism. The present calculations thus provide an atom-level explanation on these experiment results. That is, the formation of a new O* atom through the OH* dissociation in the highly stable OH*-H2O* networks is difficult and has to overcome a considerably high

Figure 10. U ≈ θO dependence for O* formation calculated with eq 15.

for the OH* → O* transformation during the growth of O* phase, which seems to further justify the MFA treatment of O* formation at relatively low coverages. The O* formation potential becomes more positive than the onset value as the coverage becomes larger than ca. 0.15 ML. This should be due to the repulsion interaction between the neighboring O* clusters at high coverages. The complete dissociation of OH* would yield a √3 × √3-structured O* phase with coverage of 1/3 ML (Figure S3). According to the calculated U0-θ dependence, this occurs at a potential of ca. 1.07 V, which is also close to the end of the region II of CV (Figure 1a). eU = dE /d(NθO) + 1/2μH2 − μH2O + ZPEO * + TkB ln(θO/(1/3 − θO)) 20456

(15) DOI: 10.1021/acsami.6b04545 ACS Appl. Mater. Interfaces 2016, 8, 20448−20458


■

It should be pointed out that the U ∼ θO dependence shown in Figure 10 may have some inaccuracy at medium coverages, e.g., when 0.05 ML

Theoretical Analysis of Electrochemical Formation and Phase

Recommend Documents