Ab initio atomistic thermodynamics study of the early stages of Cu ( 100 ) oxidation

Wissam A. Saidi,1,* Minyoung Lee,2 Liang Li,3 Guangwen Zhou,3 and Alan J. H. McGaughey2,† 1Department of Chemical and Petroleum Engineering, University of Pittsburgh, Pittsburgh, Pennsylvania 15261, United States 2Department of Mechanical Engineering, Carnegie Mellon University, Pittsburgh, Pennsylvania 15213, United States 3Department of Mechanical Engineering, State University of New York, Binghamton, New York 13902, United States (Received 14 September 2012; revised manuscript received 28 November 2012; published 26 December 2012)


I. INTRODUCTION
Oxidation of metals and metallic alloys plays a crucial role in many technologically important materials, processes, and devices including high-temperature resistant coatings, heterogeneous catalysis, and microelectronics.In particular, oxidation of copper surfaces is of fundamental and practical importance.Copper, CuO, and Cu 2 O are candidates to replace rare and expensive noble metal catalysts for applications including methane synthesis, 1-3 catalytic conversion of nitrogen oxides, 4 water-gas shift, 5,6 and preventing CO poisoning in fuel cells. 7,8he mechanisms manifested in the early stages of copper oxidation are not well understood, especially under different environmental conditions.For example, using in situ ultrahigh vacuum transmission electron microscopy, Zhou and Yang found that the substrate temperature affects the morphology of the Cu 2 O islands that grow on the Cu(100) surface during oxidation, resulting in triangular, hut, rod, and pyramid shapes. 9][15] The early stages of Cu(100) oxidation have been investigated by experimental and computational studies.  The tial oxidation of clean Cu(100) starts when oxygen molecules dissociate and atomic oxygen adsorbs at the facecentered-cubic (fcc) hollow sites. 32Upon further oxidation, oxygen coverage increases until the c(2 × 2) phase is formed.This state then transitions into the (2 √ 2 × √ 2)R45 • missingrow reconstruction (MRR) associated with an ejection of every fourth copper atom from the surface. 19Further oxygen adsorption may lead to Cu 2 O formation.Using scanning tunneling microscopy (STM), Jensen et al. observed the MRR and inferred a "squeezing out" model to explain the release of every fourth row of copper atoms. 16Fujita et al. observed, using STM and low-energy electron diffraction, c(2 × 2) domains with a maximum size of 1 × 1.5 nm 2 . 19][22][23][24][25][26] Despite numerous investigations, the mechanism for the transformation between the MRR and Cu 2 O is not clear.Using x-ray photoelectron spectroscopy (XPS), x-ray induced Auger electron spectroscopy (XAES), and STM, Lampimaki • domains that are 1.8 Å higher than their surroundings, which they attributed to the presence of trapped copper adatoms. 29Lee and McGaughey showed using density function theory (DFT) calculations that the elevation can also be explained by the presence of subsurface oxygen. 31This observation is supported by the finding that subsurface oxygen leads to a large increase in the stability of the missing-row reconstructed Cu(100) surface (consistent with Kangas et al. 37 ) and encourages the formation of Cu 2 O-like structures. 33Experimentally, using XPS, XAES, and STM, Lahtonen et al. found that subsurface oxygen leads to the growth of disordered Cu 2 O. 30 Recently, Li et al. proposed a new precursor to the onset of bulk oxidation of Cu(100) that is defined by merged missing-row nanodomains. 35Using DFT calculations, they showed that boundaries formed from merged missing-row nanodomains mismatched by a half unit-cell lead to preferred oxygen adsorption at the subsurface tetrahedral sites.The Cu-O tetrahedrons along the domain boundary strikingly resemble Cu 2 O.This finding is particularly interesting as the oxygen atoms in Cu 2 O reside in the tetrahedral sites [(1/4,1/4,1/4) and (3/4,3/4,3/4)] of the fcc copper lattice, while previous studies found that subsurface oxygen atoms in the missing-row reconstruction prefer octahedral rather than tetrahedral sites. 33FT calculations correspond to zero-temperature and zeropressure conditions.To extend DFT results to the finite temperature and pressure conditions encountered in experiments, we herein apply ab initio atomistic thermodynamics.9][40][41][42][43][44][45][46] Duan et al. previously investigated the stability of atomic oxygen-covered Cu(100) surfaces at finite temperatures and pressures up to the MRR using ab initio thermodynamics. 32They showed that the MRR is the most stable surface structure before bulk Cu 2 O formation but did not investigate the transition from the MRR to Cu 2 O island formation.
One challenge in studying the MRR to Cu 2 O transition is that kinetic hindrance could limit the full oxidation of Cu 2 O, and thus, less thermodynamically stable phases might be manifested experimentally.Previously, a strong kinetic hindrance to the bulk oxide formation of Pd(100) was identified at temperatures as high as 675 K. 47 Additionally, it is believed that the MRR is kinetically limited to 0.5 monolayer (ML) coverage, as it is found to remain reactive towards O 2 dissociation. 30New surface structures leading to subsurface oxide formation can be obtained by increasing the oxygen impingement rate on the surface. 30The objective of this study is to integrate more than sixty structures using the ab initio thermodynamics approach in order to show potential precursors for the transformation from the MRR to Cu 2 O island formation at finite temperature and pressure conditions.In particular, we discuss the transition from the oxygen chemissorption layer to subsurface oxide formation, a topic that is largely unexplored.
All slab structures include five layers where the bottom layer is fixed.The slab is thick enough in our calculations, as evidenced by noting that surface restructuring happens in the first and second copper layers only, while the third and fourth layers do not show any noticeable structural changes from those of the bulk fcc copper phase.Periodic images along the direction perpendicular to the surface are separated by 11 Å of vacuum.In all of our calculations we looked for adsorption on one side of the Cu(100) slab only.We used an 8 × 8 × 1 kpoint mesh generated using the Monkhorst-Pack scheme. 54he density of the k-point mesh, the energy cutoff, and the number of copper layers are chosen based on convergence tests performed on the p(2 √ 2 × 2 √ 2) missing-row reconstructed surface, which was reported previously. 33For example, the energy difference in the surface energy compared to more strict options obtained using a 600 eV energy cutoff, a 12 × 12 × 1 k-point mesh, and eight copper layers are 5 meV/ Å2 for the energy cutoff and k-point mesh, and 0.8 meV/ Å2 for the number of layers.Finally, all of our calculations are spin-averaged except for those of the oxygen atom and dimer, which are done using spin-polarized calculations.
To study the finite temperature and pressure conditions between a Cu(100) surface and a gas-phase oxygen environment and to determine the most stable surface structure for each condition, we use the ab initio atomistic thermodynamics framework developed by Reuter et al. 42 The surface is assumed to be in thermodynamical equilibrium with the gas phase environment. 32,42The copper and oxygen reservoirs are assumed to be those of bulk fcc copper and molecular oxygen.
To compare the stabilities of different surface structures, we define the surface free energy with respect to the clean Cu(100) surface, γ (T ,p O 2 ), as Here, G Cu/O , N Cu , and N O are the Gibbs free energy and the number of copper and oxygen atoms of the oxidized surface, all measured with respect to the clean Cu(100) surface with the same surface area, A surf .μ Cu and μ O are the copper and oxygen chemical potentials, respectively, which measure the energy penalty for exchanging particles with the atomic reservoirs.
The system temperature and oxygen partial pressure are T and p O 2 , respectively.We approximate the Gibbs free energy G by the total DFT energy, as typically done in ab initio thermodynamic analysis. 42This approximation is justifiable because solids are incompressible so the pressure contribution pV to the Gibbs free energy is negligible.Additionally, the contributions from vibrational energy and configurational disorder are also negligible.To verify this for our systems, we determined the vibrational modes using an Einstein model for two of the most stable structures that we investigated (MRR and 2hol7, see Fig. 4).We found that their vibrational contributions to γ (T ,p O 2 ) for temperatures below 900 K is less than 45 meV/ Å2 .The configurational disorder approximated using the mixing entropy is less than 2 meV/ Å2 for temperatures less than 900 K. Thus it is justifiable to write Eq. (1) as where E O/Cu and E Cu(100) are the DFT energies of the copperoxide and Cu(100) slabs.The chemical potential of Cu, μ Cu , is approximated by the total DFT energy per unit cell of copper in bulk phase, E bulk Cu .The dependence of the surface free energy on the environmental variables follows from the dependence of the oxygen chemical potential on temperature and partial oxygen pressure.In the ideal gas approximation, where k B is the Boltzmann constant, E O 2 is the total energy of a free O 2 molecule, and μ O 2 (T ,p 0 ) is the chemical potential of O 2 .In our calculations, we set the standard pressure p 0 = 1 atm and use the oxygen chemical potential values μ O 2 (T ,p 0 ) calculated by Reuter et al. 42 based on thermochemical tables. 55he surface free energy can be written succinctly as where we have shifted the chemical potential of oxygen by the total energy of its most stable phase (O 2 gas), i.e., μ O = μ O − 1 2 E O 2 , and we have defined E B , the average binding energy of oxygen on the oxide surface to be Equation ( 4) is similar to what was presented in Ref. 32 except that in our case we normalize γ (T ,p O 2 ) by the surface area A surf and not 2A surf , as atoms are adsorbed only on one side of the Cu(100) slab in our calculations.

A. Overview
We herein investigate more than sixty surface structures to develop a better understanding of the early stages of Cu(100) oxidation using the ab initio thermodynamics approach.Some of these structures were published before. 31,32,34,35Although the thermodynamic stability of the Cu(100) surface was studied recently, 32 here we examine a wider variety of surface compositions, including ones beyond the MRR, to describe possible precursor phases for bulk Cu 2 O under strong oxidizing conditions.
We divide our list of structures into three categories: (1) structures that are seen experimentally (see Sec. III B), (2) structures predicted from DFT calculations related to the MRR but with on-surface or subsurface oxygen defects as well as copper and oxygen vacancies (see Sec. III C), and (3) structures predicted from DFT calculations describing a boundary between two missing-row nanodomains (see Sec. III D).We will integrate all of the results using the ab initio thermodynamic approach in Sec.IV and show the most stable phases under different thermodynamic equilibrium conditions.Motivated by previous experimental and computational results, we will also discuss structures that, while not being thermodynamically favored, may exist in the oxidation process due to kinetic hindrance.In Figs.1(a)-1(d), we show phases that have been observed experimentally during the oxidation of Cu(100): clean, 0.25 ML, c(2 × 2), 19 and MRR. 16Structures with intermediate oxygen coverages of 0.125 ML and 0.375 ML are shown in Figs.1(e) and 1(f).
Experimentally, the MRR, Fig. 1(d), is the most frequently observed intermediate state before Cu 2 O island formation.At high temperature and low-oxygen partial pressure conditions, however, the MRR is not observed experimentally and the Cu(100) surface remains unreconstructed. 27We therefore also consider oxygen coverages on the unreconstructed surface higher than 0.5 ML to have a full account of stable phase at different oxidation environments [see Figs.1(g) and 1(h)]. 37dditionally, during the transition from the MRR to bulk Cu 2 O island formation, Iddir et al. report that there is coppervacancy exchange above 473 K, which leads to a MRR with a 1/4 disordered vacancy [c(2 × 2)/0.25 DV] structure (referred to as "spaced vacancy structure" in Ref. 27).This phase is shown in Fig. 1(i).
To facilitate the thermodynamic analysis later, the average binding energy of oxygen, E B , the number of oxygen atoms in the unit cell, and the unit cell area (A surf ) are summarized in Table I.The surface free energy as a function of the oxygen chemical potential, Eq. ( 4), can be readily generated using this information.

C. Variation on MRR with oxygen adatom, penetration, and defect
The results presented in Sec.III B are not sufficient to elucidate the mechanisms of Cu 2 O formation.Recently, Lahtonen et al. found using XPS, XAES, and STM measurements that subsurface oxygen leads to the growth of disordered Cu 2 O. 30 In these experiments, as they increased the oxygen partial pressure from 8 × 10 −10 to 4 × 10 −5 atm, well ordered MRR regions were observed and the amount of oxygen on the surface increased from 0.4 to 0.6 ML.By increasing the oxygen impingement rate on the surface, the MRR surface became reactive to further oxidation, resulting in oxide structures that the surface-oxide energy for each of these structures with an additional oxygen atom.They found that the surfaceoxide energy decreases when subsurface oxygen is present. 33f particular interest is 2octa1234(fcc), as it shows a Cu 2 O-like structure [see Fig. 2(g)].
Another subset of studied structures include oxygen vacancies that can form on oxygen-covered Cu(100) surfaces, as seen experimentally. 27Additionally, copper diffusion has a low-energy barrier 56 and thus a copper vacancy induced by diffusion is a low-energy event.To test the effects of copper and/or oxygen vacancies on the stability of the missing-row reconstructed surface, we also investigate the MRR with additional copper and/or oxygen vacancies [see Figs.2(h)-2(m)].These structures are described as MRR with copper vacancy (h) not adjacent (Cu C ) and (i) adjacent (Cu S ) to the missing row; (j) MRR with an oxygen vacancy; MRR with two oxygen vacancies located (k) horizontally V H 2O , (l) vertically V V 2O , and (m) diagonally V D 2O .The information related to the copper oxide surfaces shown in Fig. 2

is summarized in Table
2) surface unit cell was used in all of these calculations.

D. Merged missing-row nanodomains
The last set of structures describe the boundary phases formed from merged missing-row nanodomains and are shown in Fig. 3. 35 In this case, there are six possible boundaries, depending on whether the missing rows of the two domains are perpendicular or parallel to each other and on the lattice mismatch between them.In the case of parallel orientation, there are two possible boundaries denoted BMRR 1 and BMRR 2 .In the case of perpendicular orientation, there are four possible boundaries, denoted BMRR ⊥ i with i = 1, . . .,4, which differ by the distance between the domain boundaries and the nearest [100] missing row.This distance can be 1/2, 1, 3/2, or 2 unit cell lengths.Subsurface oxygen atoms in all of the BMRR ⊥ i type boundaries as well as BMRR 2 showed a preference for octahedral site adsorption.On the other hand, the subsurface oxygen in BMRR 1 preferred tetrahedral site adsorption.Results related to these structures are summarized in Table III.

IV. DISCUSSION
The formation energy of Cu 2 O predicted from our DFT calculations is 1.26 eV, which is in good agreement with the value of 1.24 eV reported by Duan et al. 32 The DFT/PBE-predicted binding energy of O 2 molecules is 5.80 eV using ultrasoft pseudopotentials while the experimental binding energy of O 2 is 5.16 eV. 57Oxygen overbinding is a known problem with standard density functionals, 58 which is troublesome as it will influence the thermodynamic analysis and the stability of the phases.Different solutions have been proposed to remedy this issue, such as using the experimental binding energy of O 2 or the experimental cohesive energies of oxides to re-define the total energy of O 2 . 59We choose to use the total energy of O 2 as computed using the employed density functional in our study.One advantage of this approach is that there will be some cancellation of errors when computing energy differences.
As shown in Eq. ( 4), the surface free energy γ depends only on the oxygen chemical potential.The upper limit for μ O is determined by the heat of formation of Cu 2 O, −1.26 eV, above which a phase transformation to Cu 2 O takes place.In our results, however, we will show the stable phases beyond this limit.Our motivation is that the full oxidation to Cu 2 O might be hindered due to kinetic reasons, and thus, less thermodynamically stable phases might be observed experimentally.In fact, the existence of energy barriers in the oxidation of metals have been reported before. 47Furthermore, there are intrinsic errors in the DFT computed energies (e.g., the O 2 binding energy).These errors have the potential to modify the location of the boundaries between different surface compositions in the phase diagram.For example, using the experimentally corrected total energy of O 2 instead of the DFT computed value, the heat of formation of Cu 2 O will be −1.61eV instead of −1.26 eV.We note, however, that small modifications to the phase diagram do not affect our identification of potential precursors to Cu 2 O.
The surface free energies plotted in Fig. 4 summarize the regions of stability of the different surface compositions and the phases that are likely to be observed experimentally  clean Cu(100) surface transitions to the 0.25 ML phase and then to the MRR.This finding agrees with previous results obtained by Duan et al. 32 The phase transitions occur at μ O (T ,p O 2 ) values of −2.12 eV (clean → 0.25 ML), and −1.58 eV (0.25 ML → MRR).Our values of μ O (T ,p O 2 ) are in good agreement with previous results which found that the clean → 0.25 ML transition happens at −1.89 eV and that the 0.25 ML → MRR transition happens at −1.72 eV. 32The small differences between the two set of results are expected considering the differences in the computational approaches.For example, Duan et al. used a double-plus-polarization quality numerical basis set, while we used a plane wave basis set.Additionally, Duan et al. used a different flavor of the generalized gradient approximation functional than used in this study.
As shown in Table I, the average binding energy of oxygen, E B , decreases as the oxygen coverage increases from 0 ML to 0.5 ML.The 0.125 ML and 0.375 ML oxygen-covered surfaces, however, are not found to be the most stable structure at any finite temperature and pressure conditions.Furthermore, while the c(2 × 2), MRR, c(2 × 2)/0.25 DV, and 1 ML p(2 × 2) phases have the same number of oxygen atoms [i.e., the same slope in Fig. 4], the MRR is the most stable as it has the highest value of E B .However, both c(2 × 2) and c(2 × 2)/0.25 DV phases have slightly lower E B values by 0.03 and 0.06 eV/ Å2 compared to MRR.Thus, these two phases could be observed experimentally in addition to the MRR.In fact, Fujita et al. observed that the MRR is initiated from the c(2 × 2) reconstruction in their experiments. 19On the other hand, the 1 ML p(2 × 2) phase is not likely to be observed at any temperature or pressure, as it is more than 1.1 eV/ Å2 less stable than the MRR.
For higher values of the oxygen chemical potential, the surface energy diagram, Fig. 4, shows that the BMRR ⊥ 2 is the most stable phase although its region of stability is narrow, between oxygen chemical potentials of −1.055 and −1.037 eV.Interestingly, BMRR 1 , which is the only boundary structure that shows the tetrahedral site as the preferred subsurface oxygen adsorption site, is not thermodynamically stable.To further describe the regions of stability of the competitive boundary phases of Sec.III D, we plot in Fig. 5 γ = γ − γ MRR , the surface energies of these phases measured with respect to the MRR phase, as a function of the oxygen chemical potential.As can be seen from the figure, the BMRR 1 phase is the least thermodynamically stable among the boundary phases and, in particular, under strong oxidizing conditions.Thus the thermodynamic analysis suggests that the BMRR 1 phase is the least likely to transform into Cu 2 O phase under equilibrium conditions.
The phase boundaries, however, are formed from the merging of randomly nucleated two-dimensional missing-row domains, which suggests that the presence of the different boundary structures is random.Although some boundary structures are less stable, transformation from one boundary structure to another is kinetically hindered because of the required massive surface structure change involved.This hindrance suggests that once a boundary is formed from merging missing row nanodomains, it will persist even if it is not thermodynamically stable.The thermodynamic analysis provides a good additional argument for why Cu 2 O prefers to nucleate along BMRR 1 sites.If BMRR 1 is unstable, it has a    We show the surface phase diagram in Fig. 6.In this phase diagram, the oxygen environment changes from oxygen poor (lower-right corner) into oxygen rich (upper-left corner).In addition to bulk Cu 2 O, there are three distinct stable phases (clean, 0.25 ML, and MRR), consistent with previous experimental and theoretical results. 27,32For comparison, the 0.25 ML → MRR phase boundaries from experiment (Iddir et al. 27 ) and ab initio thermodynamic calculations (Duan et al. 32 ) are shown in Fig. 6.Our result agrees quantitatively with previous results. 32Additionally, the experimentally observed phase boundary (Iddir et al.), is in reasonable agreement with the two theoretical predictions, by considering differences between experimental measurements and simulations (i.e., sample size, accessible range of length scale, etc.).Also, while there is a difference of the oxygen partial pressure range for the 0.25 M → MRR phase boundary between ours and experiment, the temperature range matches well.
The simplified structures described in Secs.III C and III D for subsurface oxygen that are most stable are probably not the same as what Lahtonen et al. observed. 30It is interesting, however, to note that the temperature (373 K) and pressure range (8 × 10 −10 atm and 4 × 10 −5 atm) reported by Lahtonen et al. for oxygen penetration into the subsurface are in the region of the phase diagram where the BMRR 1 and the MRR structures with extra on-or subsurface oxygen atoms are the structures with the lowest surface free energy (see Fig. 6).

V. SUMMARY
We applied an ab initio atomistic thermodynamics framework to DFT-predicted surface structures to investigate phase transitions during the early stages of Cu(100) oxidation.This framework allows for the consideration of finite temperature  30 previously reported on subsurface oxidation as a result of penetration of the MRR phase.and pressure conditions.We confirmed that the most stable structures during the early stages of Cu(100) oxidation are the clean, 0.25 ML, and MRR.Our predictions are consistent with previous experimental 27 and theoretical results. 32We also tested the thermodynamic stabilities of two additional sets of structures that are potential precursors to the Cu 2 O nucleated phase.The first set includes unreconstructed surfaces with higher oxygen coverage and additional copper and/or oxygen vacancies on the missing-row reconstructed surface.The second set includes boundaries between two missing-row nanodomains that was investigated recently. 35While these structures are not thermodynamically favored in the range of 100-1000 K and 10 −15 -10 5 atm, we believe, based on computational and experimental evidence, that they may be precursors to Cu 2 O nucleation as a result of kinetic hindrance.

FIG. 2
FIG. 2. (Color online) Top and side views of relaxed structures that may exist between the missing-row reconstruction (MRR) and bulk Cu 2 O nucleation.The hol and octa sites and their notations are defined in Ref. 33.All these relaxed structures are visualized using the p(2 √ 2 × 2 √ 2) unit cell.under thermodynamic equilibrium.Examining Fig. 4, we see that as the relative oxygen chemical potential, μ O (T ,p O 2 ) is increased (i.e., by decreasing T or increasing p O 2 ), the FIG. 3. (Color online) Top view of relaxed structures of the six boundary phases formed from two merged missing-row nanodomains.
larger tendency to transform to the more stable structure Cu 2 O due to oxygen subsurface adsorption.This transformation would be much more kinetically favorable than transforming to the BMRR ⊥ 2 structure.For higher values of the oxygen chemical potential beyond the narrow region of stability of BMRR ⊥ 2 , some of the missing-row structures with extra on-or subsurface oxygen atoms become favorable.These include the 2hol7, 2octa23, 2octa123, and 2octa1234(fcc) phases shown in Fig. 2. The transition to these phases takes place at oxygen chemical potentials of −1.03 eV (BMRR ⊥ 2 → 2hol7), −0.96 eV (2hol7 → 2octa23), −0.70 eV (2octa23 → 2octa123), and −0.61 eV [2octa123 → 2octa1234(fcc)].These findings indicate that as an oxygen-poor environment changes into an oxygen rich environment, some missing-row structures with subsurface oxygen phases have competitive regions of stability similar to the boundary structures between the two merged nanodomains.

FIG. 5 .
FIG. 5. (Color online) The difference in the surface free energies between the boundary phases and the MRR value as a function of the oxygen chemical potential μ O .

FIG. 6 .
FIG. 6. (Color online) Phase diagram showing the stable phases as a function of temperature and oxygen partial pressure.The phase boundaries between 0.25 ML and MRR from experiment (Iddir et al.) 27 and a previous ab initio thermodynamic calculations (Duan et al.) 32 are shown for comparison.Additionally, we show the experimental range of oxygen partial pressures (solid brown line) at T = 373 K where Lahtonen et al.30 previously reported on subsurface oxidation as a result of penetration of the MRR phase.

TABLE I .
Average binding energy of oxygen E B of potential structures during the early stages of Cu(100) oxidation.The relaxed structures are shown in Fig. 1.The unit cell surface area A surf is 26.65 Å2 for the p(2 × 2) unit cell and 53.29 Å2 for the p(2

TABLE II .
Average binding energy of oxygen E B of potential intermediate states between MRR and bulk Cu 2 O nucleation.The relaxed structures are shown in Fig. 2. The unit cell surface area A surf for the p(2 √ 2 × 2 √ 2) unit cell, used in all calculations, is 53.29 Å2 .

TABLE III .
Average binding energy of oxygen E B for the boundary phases between two missing-row nanodomains.The relaxed structures are shown in Fig. 3.The unit cell surface area A surf is 213.16Å2 .