Introduction

As the first elementary step of surface chemistry, adsorption, i.e., the binding of gaseous species to a solid surface, lays the foundations for complex phenomena, such as the natural corrosion of metals and the catalytic transformations of atoms and molecules on surfaces. Whereas most studies of adsorption focus on the nature and strength of the interaction between the bound atom or molecule (adsorbate) and the surface, an increase in gas concentration on the surface can increase interadsorbate interactions, which could significantly affect the outcomes of surface processes. [1][2][3] Such interactions could integrate adsorbed atoms into the lattice causing surface reconstructions or push them beneath the surface into the subsurface of the crystal. Subsurface adsorption is difficult to observe directly, however, numerous experiments indicate that light atoms can accumulate beneath different crystal surfaces, with important effects on surface reactions. For example, subsurface hydrogen has been shown to participate in catalytic hydrogenation reactions of olefins on nickel and palladium surfaces, and the role of subsurface oxygen has been investigated in oxide formation and catalytic oxidation reactions on platinum, palladium, copper, rhodium, and silver surfaces. [4][5][6][7][8][9][10][11] Theoretical models that accurately describe both gas-solid and coadsorbate interactions can provide valuable qualitative insight into subsurface adsorption and, in general, into the evolution of adsorption mechanisms with increase in adsorbate coverage on the surface. However, quantum-chemical computations of multiple adsorbates on a surface can become prohibitively expensive. A popular approach to balance the accuracy of concentration-dependent adsorption with the underlying computational expense is the ''lattice gas'' model in which gas particles bind to sites on a periodic grid or lattice. [12][13][14][15][16][17] Such models have been applied to study lateral interactions in adlayers on flat and stepped surfaces [18][19][20][21] as well as alloy formation. 22 Since conventional lattice-gas models treat the solid as a single plane, they exclusively focus on adsorbate interactions on the surface. Recently, we developed a three-dimensional form of the lattice-gas model, which includes binding sites on the surface, subsurface, and bulk of a crystalline solid. 23 This form enables the study of subsurface adsorption and coadsorption effects within and across different layers of the solid. We parameterized the model for the adsorption of atomic oxygen on the Ag(111) surface using density functional theory (DFT) and studied thermal oxygen distributions in silver using canonical Monte Carlo simulations. Further, we determined that surface and subsurface oxygen have notable chemical differences by computing their core-electron binding energies and projected density of states. The choice of this system is strongly motivated by past experimental and theoretical results suggesting that subsurface oxygen plays important roles in inducing surface reconstructions of Ag(111), forming silver oxide, and promoting partial oxidation reactions on silver catalysts. [24][25][26][27][28][29][30][31][32] Here, we advance our previous lattice-gas model, which contained the four most strongly binding sites, to an all-site lattice-gas model for atomic oxygen on Ag(111) that contains all types of high-symmetry sites in the solid. Similar to the previous model, we parameterize the new model using DFT without any empirical parameters or function fitting. However, we apply the new model to perform Monte Carlo simulations in two ensembles, first, in the canonical ensemble to compute oxygen distributions as a function of coverage, and second, in the grand-canonical ensemble to calculate oxygen distributions as a function of O 2 pressure. In the latter category, we investigate subsurface adsorption under low pressures typically employed in surface-science laboratories and under high pressures used in industry, covering a pressure range of ten orders of magnitude. While this model describes chemisorbed atomic oxygen on Ag(111) and cannot capture ordered surface oxides, we find that our simulated surface coverages of oxygen are similar to the oxygen coverages in the well-studied surface reconstructions of Ag(111). 29,[33][34][35] This assures us that the oxygen-oxygen interactions between the surface and subsurface that are critical in stabilizing or destabilizing subsurface deposition in Ag(111) are represented realistically by our model.

Finally, using the results from grand-canonical Monte Carlo simulations, we construct a phase diagram of O/Ag(111) to determine the pressure-temperature window within which subsurface oxygen thermodynamically co-exists with surface oxygen in Ag(111). Whereas ab initio thermodynamics approaches have been used to calculate phase diagrams of surface oxides and bulk oxides of O/Ag(111), [36][37][38] to our knowledge, our phase diagram is the first to show conditions under which subsurface oxygen is thermodynamically stable.

Method

2.1 Lattice-gas model with surface, subsurface, and bulk sites Our lattice-gas model is formulated as a truncated cluster expansion of the form

where E total is the total energy of N adsorbates, E i is the site energy of the i-th adsorbate, and DE ij is the 2-body or pairinteraction energy between the i-th and j-th adsorbates. Since O 2 binds very weakly compared to O on Ag(111) 39 and is too large to intercalate below the surface, we only consider the dissociative binding (chemisorption) of O 2 as O atoms to the Ag(111) surface. A negative E i indicates stable adsorption and a negative DE ij indicates a cooperative or ''attractive'' pair interaction that strengthens the binding of both adsorbates. We include every type of high-symmetry binding site in the adsorption model. These are the atop (A), bridge (B), hexagonal close-packed (HCP, H) hollow, and face-centered cubic (FCC, F) hollow sites on the surface; tetrahedral (T1, T2, T3), invertedtetrahedral (I1, I2, I3), octahedral sites (O1, O2, O3) in the two subsurface and near-bulk regions, respectively; and tetrahedral (T) and octahedral sites (O) in the bulk of the crystal (Fig. S1, ESI †).

In the interior of the crystal, the inverted-tetrahedral site is geometrically and energetically identical to a tetrahedral site, therefore, we only consider tetrahedral and octahedral sites for oxygen absorption in the bulk of silver. Overall, our Ag(111) model contains five regions to which O atoms can bind, the surface ('surf'), the first subsurface (between first and second layers, 'sub1'), the second subsurface (between second and third layers, 'sub2'), the near-bulk (between third and fourth layers), and the bulk (below the fourth layer). We exclude three-body interactions from this all-site model as we did from our 4-site lattice-gas model of O/Ag(111). 23 The error in E total from neglecting higher-body interactions is discussed in Section 2.2. Ultimately, considering all types of binding locations on the surface, in the subsurface, and inside the bulk crystal, we computed 15 types of site energies and 307 types of pair-interaction energies for O/Ag(111) using DFT. Since the Ag atoms remain fixed in their lattice positions, our model can only describe chemisorbed and dissolved oxygen, not surface or bulk oxides that have different lattice structures than Ag(111) or bulk Ag. However, our results and their connections to experiments discussed in Section 3 strongly suggest that our model provides meaningful insight into the formation of subsurface oxygen and its interactions with surface oxygen. We add that the effects of lattice relaxation are included in E i and DE ij through the geometry optimizations described in Section 2.2.

Model parameterization using periodic DFT

The values of E i and DE ij for the lattice-gas model of O/Ag(111) were computed using spin-polarized DFT within the Vienna ab initio Simulation Package (VASP, version 5.4). [40][41][42][43] Exchange and correlation between electrons were described by the revised Perdew, Burke, and Ernzerhof with Pade Approximation (RPBE) functional. [44][45][46][47] The bulk lattice constant (a) of Ag was computed using the RPBE functional in combination with a 15 Â 15 Â 15 Monkhorst-Pack 48 k-point grid to be 4.213 Å, in good agreement with the experimental lattice constant of 4.09 Å. 49 The Ag(111) surface was modeled using an orthorhombic p(3 Â 4) periodic supercell, consisting of six layers of Ag atoms and a vacuum layer twice as thick as the surface slab (Fig. 1(a)). The dimensions of the supercell were 8.93 Å Â 10.32 Å Â 43.78 Å, within which the thickness of the vacuum layer was 29.19 Å. The lowest fractional oxygen coverage, y, that could be calculated in the supercell was 1 12 % 0:08 monolayer (ML), where 1 ML contains one oxygen atom per silver atom in the top layer. Bulk Ag was modeled using the analogous p(3 Â 4) supercell without any vacuum layer. Electronic energies were computed using a plane-wave cutoff of 400 eV in combination with a 5 Â 5 Â 1 Monkhorst-Pack 48 k-point grid for the Ag(111) surface and a 5 Â 5 Â 5 Monkhorst-Pack k-point grid for bulk Ag. The projector augmented wave (PAW) method as implemented in Version 54 of the VASP PAW-PBE pseudopotentials was used to describe interactions between core electrons and nuclei, 50,51 with energy cutoffs of 250 eV and 400 eV for Ag and O, respectively. In addition, a Methfessel-Paxton 52 smearing of width 0.2 eV was applied to minimize electronic free energies. To test the sensitivity of the model parameters to the choice of pseudopotential, 53 we calculated E i for O adsorbed on the FCC site and DE ij for a pair of O atoms adsorbed on first-neighbor FCC sites using PAW-PBE pseudopotentials with varying energy cutoffs. The results (Table S1, ESI †) show that the use of higher cutoff energies change the model parameters by up to 40 meV.

To compute minimum-energy adsorption geometries, we constrained the bottom three layers of the supercell (Fig. 1(b)) at the bulk crystal geometry and allowed the top three layers and all adsorbates to relax until the forces on the mobile nuclei were less than 0.05 eV per Å. Electronic energies were dipole-corrected for charge transfer from the surface to the adsorbates. To calculate absorption energies in the bulk, absorption structures of O/Ag were optimized by allowing all Ag and O nuclei in the bulk supercell to relax down to forces less than 0.05 eV per Å. Vibrational frequencies of optimized geometries were calculated from the Hessian matrix computed using the finite-difference method. Atomic charges were calculated using Bader's method of charge density partitioning. [54][55][56] The site adsorption energy, E i , was defined as

where E O/Ag(111) was the energy of the relaxed O/Ag(111) supercell, E Ag(111) was the energy of the relaxed Ag(111) supercell, and E O was the energy of an oxygen atom inside an empty supercell.

Similarly, E ij of two coadsorbed oxygen atoms was defined as,

Using E i and E ij , the pair-interaction energy between two coadsorbed oxygen atoms (Fig. 1(c)), DE ij , was defined as

In this way, we calculated the E i and DE ij values for oxygen atoms binding to the surface, first subsurface, second subsurface, and bulk of Ag( 111 The region between the third and fourth layers of the solid was treated as the ''near bulk'' region and its parameters were calculated via linear interpolation between the parameters of the second subsurface and those of the bulk crystal. For example, the E i of oxygen in the near-bulk octahedral (O3) site (À2.74 eV) was obtained by interpolating between E i in the second-subsurface octahedral (O2) site (À2.71 eV) and E i in the bulk octahedral (O) site (À2.77 eV). Similarly, the DE ij between two oxygen atoms occupying first-neighbor O3-O3 sites (À0.16 eV) was obtained by interpolating between the corresponding DE ij of O2-O2 sites (À0.18 eV) and DE ij of O-O sites (À0.14 eV). All sites below the fourth layers were described using bulk-absorption parameters. The parameters of the lattice-gas model are provided in the ESI. †

To assess the quality of the lattice-gas model, we created a parity plot comparing 219 DFT-computed adsorption energies over a wide range of coverages and involving different binding sites to the corresponding model-derived adsorption energies. This plot (Fig. S3, ESI †) has a linear-fit equation of y = 0.073 + 0.994x with R 2 = 0.989 and root-mean-squared error (RMSE) of 0.646 eV. Further, we compared the average adsorption energy per oxygen atom between DFT and the lattice-gas model for the same 219 geometries because our canonical and grand-canonical Monte Carlo simulations, as described in Section 2.3, sample geometries by moving one oxygen atom in a single step. This second parity plot (Fig. S4, ESI †) has a linear-fit equation of y = À0.029 + 0.978x with R 2 = 0.872 and RMSE of 0.132 eV over an adsorption-energy range of À1.88 eV to À4.80 eV per oxygen atom. Importantly, the RMSE calculations of average adsorption energies at specific total coverages (Table S2, ESI †) show that the RMSE per oxygen remains stable with increasing coverage. Overall, this analysis suggests that our lattice-gas model represents DFT adsorption energies with reasonable and consistent accuracy over a wide range of oxygen coverages on Ag(111).

As in our earlier 4-site model, 23 we excluded three-body interactions between adsorbed oxygen atoms in the present, all-site lattice-gas model. While developing the 4-site model, we examined the contributions of three-body interactions by computing the total adsorption energies of 20 geometries, each containing 3 O atoms on Ag(111) using DFT. The geometries were chosen to include various combinations of the most strongly binding sites on the surface and in the subsurface of Ag(111). Then we calculated the 3-body interaction energy in each geometry as the difference between the adsorption energy from DFT and the corresponding energy from our DFT-parameterized pairwise-interaction model. This set of results, reported in Table S3 (ESI †) of our previous work, 23 showed that 3-body interaction energies constitute on average 1% of the total adsorption energy. However, as the oxygen coverage increases, the number of 3-body interactions increases more rapidly than the number of 2-body interactions and could still result in significant cumulative 3-body contributions at higher coverages. Fortunately, further examination of our calculated 3-body interactions revealed that they are fairly evenly distributed between cooperative and anticooperative. Therefore, we posit that the considerable cancellation of 3-body effects in O/Ag(111) enables our pairwise lattice-gas models to capture the DFT results reliably over the studied coverage range. This cancellation is likely why the RMSE per oxygen atom shown in Table S2 (ESI †) is fairly insensitive to coverage.

We note that vibrational free energies of adsorbed oxygen were not included in the parameters of the lattice-gas model. Consequently, the site energies, E i , and pair-interaction energies, DE ij are electronic energies, not free energies. To examine the effects of neglecting the vibrational contributions, we calculated the freeenergy corrections to E i and DE ij , respectively, for O atoms adsorbed on/in the energetically preferred sites at 298.15 K. The results (Table S3, ESI †) show that free-energy corrections are relatively small, $ 1 100 th of the respective values of E i and DE ij and, importantly, similar between the different sites. These results suggest that using adsorption free energies in the Monte Carlo simulations would produce similar thermal distributions of oxygen atoms on Ag(111) as using adsorption (electronic) energies. We clarify that our 4-site and all-site models apply the DFT-computed parameters without any model fitting to DFT calculations at different coverages.

Monte Carlo simulations

Monte Carlo (MC) simulations were implemented using the lattice-gas model to determine equilibrium oxygen distributions in Ag(111) at chosen coverages, pressures, and temperatures. The simulations were performed on a hexagonal periodic O/Ag(111) supercell consisting of ten layers, each containing 1024 Ag atoms (Fig. 1(d)). The hexagonal periodic boundary conditions relate opposite faces of the hexagon in the xy cartesian plane by symmetry. In accordance with the DFT parameterization of the lattice-gas model, the binding sites below the fourth layer of Ag atoms in the MC supercell were treated as absorption sites in the bulk of the solid. The Metropolis MC algorithm 57 was applied within the canonical and grandcanonical ensembles to sample configurations of oxygen atoms over the available adsorption sites in the solid, keeping Ag atoms fixed in their crystal positions. The purpose of the canonical (N,V,T) MC simulations was to determine how a fixed total number (N) of oxygen atoms distributes between the surface, subsurface, and bulk regions of Ag(111) at a specified temperature. Grand-canonical (m,V,T) MC simulations add and remove oxygen atoms from Ag(111) at a fixed chemical potential (m) corresponding to the partial pressure of O 2 (g). The equations used to calculate m O 2 (g) for a particular pressure and temperature are provided in the ESI. † The purpose of the grand-canonical MC simulations was to determine both the total oxygen coverages and corresponding oxygen distributions on Ag(111) at specified temperatures and applied partial O 2 pressures. Additionally, the grand-canonical MC simulations utilized replica-exchange (RE) 58 in order to improve configurational space sampling. RE was implemented by allowing the configurations between grandcanonical MC simulations performed at different but neighboring temperatures exchange every 1000 MC steps. We initialized the system configuration by randomly assigning a defined number of oxygen atoms to sites on the surface and between the first four silver layers. Each move in the MC simulation randomly chose a bound oxygen atom and moved it to a random, available site anywhere in the supercell. For simulations in the grand-canonical ensemble, an initial system configuration was randomly sampled from a converged oxygen distribution obtained from canonical ensemble simulations at a total coverage of 0.75 ML. The initial total coverage of 0.75 ML was chosen to balance the sampling of surface and subsurface oxygen adsorption with the likely bulk oxidation of the silver substrate at higher coverages (Fig. S5, ESI †). In addition to the random MC moves described above, oxygen atoms were also randomly added from the surrounding gas reservoir to empty sites in the supercell or removed from the supercell into the gas. In both canonical and grand-canonical MC simulations, ensemble averages were obtained by using ten independently initialized Markov chains, each consisting of 3 000 001 configurations for which only the last 2 million configurations were utilized for these ensemble averages. The average geometries of oxygen adsorption from the simulations were calculated using the adsorption probability of single and pairs of adsorption sites, which take into account correlations between any two coadsorbed oxygen atoms. A negative perpendicular distance between the surface and O indicates that O is below the surface. We do not report perpendicular O-surface distances or Bader charges for the near-bulk sites as those adsorption energies were calculated via interpolation as described in Section 2.2 instead of directly from DFT calculations. The Bader charges show significant charge transfer of BÀ1 from metal to adsorbate and highlight the anionic nature of adsorbed oxygen on silver. This result, combined with projected density of states from our previous work 23 showing no unpaired electron density on the O atom up to B0.6 ML, indicates that the adsorbed oxygen is qualitatively O 2À but is numerically underestimated by the Bader charge analysis as O À . We found that oxygen adsorbs most strongly to FCC and HCP hollow sites on the Ag(111) surface and to octahedral sites in the subsurface and bulk. During DFT geometry optimizations at y ¼ 1 12 ML ML, oxygen in the inverted tetrahedral site (Table 1, I1) rose to the HCP hollow site directly above it, resulting in stronger binding and a slightly positive O-surface distance. We verified that this oxygen remains in the I1 site at higher coverages. Our calculated adsorption energies agree well with previous DFT studies at low coverages, 27,[59][60][61] which were limited to ab initio investigations of adsorption on the surface and in the first subsurface. Fig. 2 shows a pixel map of first-neighbor, second-neighbor, and third-neighbor pair interactions (DE ij ) between coadsorbed oxygen atoms calculated using DFT. First-neighbor interactions are highly anti-cooperative (repulsive) on the surface but become progressively cooperative (attractive) in the subsurface and bulk. For example, first-neighbor DE ij values between like sites on the surface range from 0.27-0.53 eV, whereas the analogous first-neighbor DE ij values in the first subsurface and second subsurface range from À0.08 to -0.06 eV and À0.22 to -0.18 eV, respectively. These trends show that coadsorbate interactions on the surface are unfavorable due to the strong anionic character of surface oxygens (Table 1), but become favorable in the subsurface due to greater charge screening by silver atoms above and below the oxygen atoms. Importantly, most first-neighbor O surf -O sub1 interactions are anti-cooperative, whereas many first-neighbor O sub1 -O sub2 interactions are cooperative. Such differences in DE ij represent significant variations in O-Ag and O-O interactions between the surface and subsurface regions, which in turn have profound effects on the oxygen distributions on Ag(111).

Results and discussion

Oxygen adsorption energies and pair interactions

The DFT-computed E i and DE ij values discussed in this section were inserted into eqn (1) to complete the lattice-gas model for O/Ag(111). We note that the model adsorption energy per oxygen of À3.40 eV for a total coverage of 1 3 ML uniformly occupying FCC and HCP hollow sites agrees well with the surface adsorption energy of À3.44 eV at 0.375 ML measured by Campbell using temperature-programmed desorption (TPD). 39 Sections 3.2-3.4 discuss the equilibrium oxygen distributions and phase diagram computed by integrating the model with MC simulations.

Oxygen distributions with increasing total coverage

Fig. 3 shows atomic oxygen distributions on the surface, subsurface, and bulk of Ag(111) at different total oxygen coverages 1.

(y total ) computed using canonical MC simulations at 300 K. Tables S4 and S5 (ESI †) detail the regional coverages, site coverages, and standard deviations corresponding to Fig. 3a. Within our coverage range of 0.125 ML r y total r 1 ML, oxygen atoms primarily occupy the FCC and HCP hollow sites on the surface and octahedral sites in the subsurface regions. Notably, the surface coverage (y surf ) reaches a maximum of 0.37-0.39 ML above which excess oxygen deposits in the subsurface of Ag(111). The surface coverage saturates due to repulsions between charged oxygen atoms, which are especially strong between oxygen atoms coadsorbed on first-neighbor surface sites. We note that nearest-neighbor occupations cannot be avoided above a total coverage of 1 3 ML due to the symmetry of the FCC lattice. Fig. 3(a) shows that for y total 4 0.375 ML, oxygen deposits into the first subsurface but only up to B0.14 ML due to O surf -O sub1 repulsions between surface oxygen and first-subsurface oxygen. Thus, excess oxygen primarily collects in the second subsurface, where interactions with other O sub2 are cooperative, and interactions with oxygen atoms in layers above and below are either weakly cooperative or weakly anti-cooperative. There is negligible oxygen deposition in the near-bulk and bulk regions within the studied coverage range. Our distributions show that subsurface oxygen is thermodynamically stable only in the presence of surface oxygen at y total greater than the surface saturation coverage. The computed regional distributions change negligibly between 300-600 K because the energy differences between various bound oxygen species are greater than the available thermal energy in this temperature range (Table S6, ESI †).

Our computed saturation coverage agrees well with the surface concentrations of several oxygen-induced reconstructions or surface oxides of Ag(111), such as p(4 Â 4) containing

containing y O = 0.4 ML, and c(4 Â 8) containing y O = 0.5 ML. [33][34][35]62 Our calculated y surf decreases slightly for y total Z 0.875 ML in order to lower repulsion between O surf and O sub1 . It is encouraging that our maximum y surf captures the experimental surface oxygen concentrations without including surface reconstructions in the model. This agreement suggests that the O surf -O surf and O surf -O sub1/sub2 interactions that control surface saturation are represented realistically by our model. Further, we note that the subsurface populations shown in Fig. 3 compare well to results from TPD experiments on O/Ag(111) by the Killelea group, which show that y surf saturates at 0.375 ML but y total increases up to 1.12 ML. 63 X-ray photoelectron spectroscopy (XPS) measurements in the same work indicate that the excess oxygen of up to 0.75 ML deposits within 5.7 Å below the surface. This result aligns well with our subsurface coverage of B0.65 ML at y total = 1 ML, distributed in the first and second subsurface regions, with maximum depth of B4 Å (O-surface perpendicular distances in Table 1).

Fig. 3(b) and (c) show the average geometries at y total = 0.5 ML and 0.75 ML, respectively, obtained from the corresponding MC distributions in Fig. 3(a). The surface oxygen of y total = 0.5 ML is distributed uniformly between FCC and HCP hollow sites (B0.19 ML each) with only B0.01 ML collective occupation of the top and bridge sites, resulting in y surf = 0.39 ML. The first subsurface contains B0.04 ML of oxygen of which B0.03 ML occupy octahedral sites, whereas the second subsurface contains B0.07 ML of which B0.06 ML occupy octahedral sites. At a higher y total = 0.75 ML, the surface oxygen coverage and distribution are similar to the previous case with y surf of 0.38 ML mostly populating FCC (B0.16 ML) and HCP (B0.19 ML) sites. The first-subsurface oxygen population increases to B0.08 ML and consists of B0.06 ML in octahedral sites, whereas the second-subsurface oxygen increases to B0.29 ML and consists of B0.24 ML in octahedral sites. The site coverages from the distributions in Fig. 3(a) are listed in Table S3 (ESI †). Interestingly, Fig. 3(c) shows that the oxygen atoms in the second subsurface cluster together due to cooperative oxygen-oxygen interactions in this region, as discussed in Section 3.1 and depicted in Fig. 2. Based on this result, we propose that the formation of such oxygen islands in the second subsurface is a precursor to the formation of bulk silver oxide and treat our simulated distributions for y total 4 0.75 ML as proxies for the Ag 2 O oxide. Interestingly, our prediction is corroborated by Lundgren and coworkers whose XPS results on O/Ag(111) suggest that Ag 2 O starts to form for y total Z 0.75 ML. 64 We return to this point of bulk oxidation in Section 3.3. Finally, apart from relatively minor quantitative differences, our simulated distributions agree with the corresponding distributions derived from the 4-site model, which only included the most favorable binding sites, i.e., FCC, HCP, and octahedral sites. This agreement demonstrates that these sites collect most of the oxygen atoms on or in Ag(111) even when other types of sites are available. However, an important change due to advancing to the all-site model is that the spread of oxygen concentrations to other types of sites modulates the cumulative pair-interaction energy between the surface and first subsurface and increases the maximum oxygen coverage in the first subsurface from B0.05 ML to B0.14 ML.

Although the MC simulations demonstrate that subsurface oxygen is thermodynamically allowed for intermediate to high coverages, they do not shed light on the kinetics of oxygen diffusion into the subsurface of Ag(111). We performed nudged-elastic-band computations within VASP to determine the energy barriers to percolation of the oxygen atom from the surface to the second subsurface of Ag(111) at the lowest coverage, y total = 0.08 ML. Specifically, we calculated the vertical diffusion of an oxygen atom from the FCC site to the octahedral site (O1) directly below in the first subsurface, followed by lateral diffusion from the octahedral site to the tetrahedral site (T1), and finally vertical diffusion from the tetrahedral site to the octahedral site (O2) directly below in the second subsurface. We chose this path because it involves the energetically preferred surface and subsurface sites for O on Ag(111) (Table 1). The calculated diffusion path (Fig. S6(a) (ESI †)) shows that barriers into and within the first subsurface (FCC -O1 and O1 -T1) are very close to the adsorption-energy differences between the initial and final sites, and diffusion into the second subsurface (T1 -O2) has a relatively low barrier of B0.3 eV. In addition, we calculated the adsorption energy of the diffusing oxygen at the four sites in the presence of a fixed y surf of 0.42 ML (y total = 0.5 ML) using DFT. The results (Fig. S6(b) (ESI †)) demonstrate how the presence of a fixed y surf similar to the surface saturation coverage transforms the ''ascending'' staircase of adsorption energies at y total = 0.08 ML into a ''descending'' staircase at y total = 0.5 ML, suggesting that the barrier to oxygen diffusion from the surface into the first subsurface of Ag(111) decreases significantly with increasing surface coverage.

Oxygen distributions with increasing oxygen pressure

While canonical MC simulations can provide valuable insight into how a fixed total oxygen coverage thermally populates the surface and subsurface of Ag(111) and distributes between different types of binding sites, it cannot capture the exchange of oxygen atoms between Ag(111) and the surrounding O 2 gas. To this end, we performed grand-canonical MC in which Ag( 111) is held at equilibrium with a O 2 gas reservoir, enabling the exchange of atomic oxygen between the gas and the solid. If the chemical potential of an oxygen atom is lower (more favorable) in O 2 (g) than on Ag(111) then it is more likely to transfer to the gas, whereas if it is lower on Ag(111) then it is more likely to bind to the solid. Fig. 4(a) shows the oxygen regional distributions on Ag(111) from 300-450 K at a low O 2 gas pressure of 5 Â 10 À7 torr (6.67 Â 10 À10 bar), chosen to represent the vacuum to ultrahigh vacuum conditions employed in a surface-science laboratory. Fig. 4(b) shows the progressive decrease in oxygen coverage in the temperature subset from 425-450 K at the same O 2 (g) pressure of 5 Â 10 À7 torr. Tables S7 and S8 (ESI †) provide the regional coverages and respective standard deviations corresponding to Fig. 4.

Oxygen distributions in Fig. 4 show that surface adsorption on Ag(111) becomes thermodynamically favorable below a temperature of 450 K. Above 450 K, oxygen desorbs from the surface due to a lower chemical potential of 1 2 O 2 (g) relative to the chemical potential of O surf . These results are consistent with TPD results from Campbell, 39 which show up to only B0.03 ML of oxygen atoms deposited on Ag(111) at 490 K under UHV O 2 pressures. Later, Schlo ¨gl and coworkers observed using scanning tunneling microscopy (STM) that the Ag(111) surface remained essentially clean down to 440 K under UHV O 2 pressures. 65 Further, we find that y total does not exceed the surface saturation coverage between 300-450 K because the chemical potential of subsurface oxygen remains greater than the chemical potential of 1 2 O 2 (g). As a result, oxygen in excess of the saturation surface coverage does not percolate beneath the surface but rather desorbs into the gas under these temperatures and pressure.

In sharp contrast, at higher partial O 2 pressures, such as p O 2 = 0.01 bar, 0.1 bar, and 1 bar, oxygen deposits on Ag(111) between 400-550 K and saturates the surface over the entire temperature range (Fig. 5 and Tables S9-S11 (ESI †)). The results at 0.01 bar compare favorably to TPD and low-energy electron diffraction (LEED) results from Campbell, which show the p(4 Â 4) surface reconstruction (y O = 0.375 ML) at 490 K and p O 2 = 5 Torr E 0.007 bar. 39 Using high-pressure STM, Schlo ¨gl and coworkers observed p(4 Â 4) and p 4 Â 5 ffiffi ffi 3 p À Á (y O = 0.375 ML) surface reconstructions between 450-480 K for p O 2 Z 1 mbar. 65 Further, our simulations show the presence of subsurface oxygen within 400-550 K at all three studied pressures. Therefore, we find that the low-pressure (p O 2 = 6.67 Â 10 À10 bar) and high-pressure (p O 2 = 0.01-1 bar) regimes produce qualitatively different adsorption outcomes due to the occurrence of subsurface adsorption at higher O 2 pressures.

Phase diagram of subsurface oxygen in Ag(111)

To determine the temperature ranges within which subsurface oxygen is stable at higher O 2 pressures, we fit the computed y total at the three pressures from Fig. 5 to exponential decay functions of temperature, as shown in Fig. S7 (ESI †). These fits and the corresponding equations can be used to calculate the temperature corresponding to a specific y total for p O 2 = 0.01 bar, 0.1 bar, and 1 bar. Based on our grand-canonical MC distributions from Fig. 5 (Tables S9-S11, ESI †), we used 0.46 ML as the y total with an initial, small concentration of subsurface oxygen, and based on Fig. 3 with the accompanying discussion in Section 3.2, we used 0.75 ML as the y total above which the bulk Ag 2 O oxide becomes thermodynamically favored over dissolved oxygen. Using the fits in combination with these coverage thresholds, we calculated that subsurface oxygen is thermodynamically stable between 394.4-426.0 K at p O 2 = 0.01 bar, 416.8-472.0 K at p O 2 = 0.1 bar, and 451.3-505.1 K at p O 2 = 1 bar. For temperatures above each range and up to 550 K, atomic oxygen is bound to the surface of Ag(111) either as chemisorbed oxygen or as a surface oxide, and for temperatures below each range, the thermodynamically preferred phase is Ag 2 O. We note that our computed equilibrium temperatures (T eqm ) for Ag 2 O formation agree well with T eqm for the reaction, , as shown in eqn (5). For subsurface adsorption (red points), a = 507.8 K and b = À0.0907, and for Ag 2 O formation (black points), a = 450.1 K and b = À0.0728. The blue dashed horizontal line demarcates the lower limit of the temperature range of 475-550 K commonly used in industrial partial oxidation reactions.

2Ag!sÞ þ 1 2 O 2 !gÞ ! Ag 2 O!sÞ, calculated to be 362.8 K at p O 2 = 0.01 bar, 408.2 K at p O 2 = 0.1 bar, and 466.7 K at p O 2 = 1 bar using thermodynamic properties. 66 Our T eqm values for Ag 2 O formation also compare well to those calculated in previous studies using ab initio thermodynamics. [36][37][38] Finally, in Fig. 6, we plot the calculated equilibrium temperatures as functions of log p

Conclusions

In summary, we have extended the theoretical framework of lattice-gas models from two to three dimensions to investigate how coadsorption effects with increasing surface coverage lead to subsurface adsorption in crystalline solids. Whereas we previously developed a 4-site ab initio lattice-gas model to explore the interplay between surface and subsurface adsorption of atomic oxygen on the Ag(111) surface, here we advance our treatment to an all-site ab initio model and integrate it with not only canonical but grand-canonical Monte Carlo simulations to examine oxygen adsorption over wide ranges of coverage, pressure, and temperature. We find that subsurface oxygen does not form in Ag(111) under high vacuum or lower O 2 pressures even at room temperature, but coexists with surface oxygen under the high pressures and temperature ranges used in industrial partial oxidation reactions. Therefore, our study reveals that oxygen-oxygen interactions on silver lead to qualitatively distinct adsorption behaviors across different temperature and pressure regimes. Using the Monte Carlo results, we construct the first phase diagram of O/Ag(111) that includes subsurface oxygen and show the temperatures and pressures under which subsurface oxygen is thermodynamically stable in Ag(111). We hope that our insight can be applied to harness the availability, properties, and reactivity of subsurface oxygen for various applications. Future work will develop lattice-gas models to explore how cooperative effects between surface and subsurface adsorbates at different coverages influence surface reconstructions of crystalline solids.