<?xml-model href='http://www.tei-c.org/release/xml/tei/custom/schema/relaxng/tei_all.rng' schematypens='http://relaxng.org/ns/structure/1.0'?><TEI xmlns="http://www.tei-c.org/ns/1.0">
	<teiHeader>
		<fileDesc>
			<titleStmt><title level='a'>Insight into Subsurface Adsorption Derived from a Lattice-Gas Model and Monte Carlo Simulations</title></titleStmt>
			<publicationStmt>
				<publisher></publisher>
				<date>03/24/2022</date>
			</publicationStmt>
			<sourceDesc>
				<bibl> 
					<idno type="par_id">10320492</idno>
					<idno type="doi">10.1021/acs.jpcc.2c00342</idno>
					<title level='j'>The Journal of Physical Chemistry C</title>
<idno>1932-7447</idno>
<biblScope unit="volume">126</biblScope>
<biblScope unit="issue">11</biblScope>					

					<author>Carson J. Mize</author><author>Lonnie D. Crosby</author><author>Sara B. Isbill</author><author>Sharani Roy</author>
				</bibl>
			</sourceDesc>
		</fileDesc>
		<profileDesc>
			<abstract><ab><![CDATA[Theoretical gas-surface models that describe adsorption over a broad range of adsorbate concentrations can provide qualitative insight into chemical phenomena, such as subsurface adsorption, surface reconstruction, and industrial heterogeneous catalysis. However, most atomistic, quantum-mechanical models of gas-surface adsorption are limited to low adsorbate coverage due to the large computational cost of models built using many surface atoms and adsorbates. To investigate adsorption in the subsurface of a crystalline solid with increasing coverage, we present a lattice-gas adsorption model that includes surface and subsurface sites of the solid and is fully parametrized using density functional theory. We apply the model to study the competition between the surface and subsurface adsorption of atomic oxygen on the Ag(111) surface. Oxygen population distributions calculated using the model in combination with Monte Carlo simulations show the onset of subsurface adsorption above a total coverage of 0.375 monolayer and a greater accumulation of oxygen in the second than in the first subsurface at total coverages between 0.5 and 1 monolayer. Our simulations also show that oxygen atoms do not percolate into the bulk region of silver for total coverages of up to 1 monolayer, indicating that oxygen adsorbed in the subsurface is distinct from oxygen absorbed in the bulk in this coverage range. Computations of core-electron binding energies and projected density of states for the equilibrium oxygen distribution at 0.5 monolayer reveal qualitative differences in the oxygen-silver bonding at the surface and subsurface, suggesting that oxygen adsorbed in the two regions could play distinct roles in the surface chemistry.]]></ab></abstract>
		</profileDesc>
	</teiHeader>
	<text><body xmlns="http://www.tei-c.org/ns/1.0" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xmlns:xlink="http://www.w3.org/1999/xlink">
<div xmlns="http://www.tei-c.org/ns/1.0"><head n="1.">INTRODUCTION</head><p>The adsorption of an atom or molecule on a solid surface is an elementary step of many chemical processes, ranging from corrosion, gas storage, and heterogeneous catalysis to chemical sensing and the formation of self-assembled monolayers. Whereas advances in periodic density functional theory (DFT) have enabled the quantum-mechanical computation of adsorption properties at low adsorbate concentrations with high accuracy and efficiency, large adsorbate concentrations remain a challenge to compute due to their complex surface chemistry and steep computational cost. The cost can be kept low by keeping the number of adsorbates fixed and decreasing the lateral dimensions of the supercell, but the coverages constructed in this matter enforce highly periodic adlayer structures. Chemically meaningful, flexible, yet efficient theoretical models capable of describing adsorption over a wide range of adsorbate concentrations can provide mechanistic insight into processes such as surface reconstruction, surface oxidation, and industrial heterogeneous catalysis.</p><p>Lattice-gas models and cluster-expansion methods are powerful theoretical approaches to studying coverage-dependent adsorption and the effects of interadsorbate interactions as coverage is increased. These approaches describe the total adsorption energy as the sum of the noninteracting energies of the adsorbates at different lattice sites plus the n-body interaction energies between the adsorbates up to an appropriate and affordable cutoff for n. <ref type="bibr">[1]</ref><ref type="bibr">[2]</ref><ref type="bibr">[3]</ref><ref type="bibr">[4]</ref><ref type="bibr">[5]</ref> Some recent examples in which the site-adsorption energies and interaction energies were calculated using DFT include coverage-dependent adsorption studies of hydrogen on Fe(100) <ref type="bibr">6</ref> and atomic oxygen on Pt(111) <ref type="bibr">7</ref> and lateral adsorbate interactions in the NO-CO reaction system on Rh(100) and Rh(111). <ref type="bibr">8</ref> Clusterexpansion methods are also used in conjunction with Monte Carlo and kinetic Monte Carlo simulations, such as in the studies of ethylene epoxidation on Ag(100), Ag(110), and Ag(111) <ref type="bibr">9</ref> and adlayer phases in H/Cu(100), H/Ag(100), and O/Cu(100). <ref type="bibr">10</ref> Whereas lattice-gas models have been created to study adsorption in diverse gas-solid systems, to our knowledge, they have exclusively addressed adsorption on surface sites and have not explored subsurface adsorption, that is, adsorption just beneath the surface of a crystalline solid. <ref type="bibr">11</ref> The subsurface adsorption of light atoms, such as hydrogen, carbon, and oxygen, has been investigated on several transition metals, including aluminum, nickel, copper, rhodium, palladium, silver, and platinum, to understand their role in catalytic reactions, surface reconstruction, and surface oxidation. However, the conditions of subsurface adsorption, the chemical properties of subsurface adsorbates, and the microscopic mechanisms of their participation in surface processes are not yet completely understood. <ref type="bibr">[12]</ref><ref type="bibr">[13]</ref><ref type="bibr">[14]</ref><ref type="bibr">[15]</ref><ref type="bibr">[16]</ref><ref type="bibr">[17]</ref><ref type="bibr">[18]</ref><ref type="bibr">[19]</ref><ref type="bibr">[20]</ref><ref type="bibr">[21]</ref><ref type="bibr">[22]</ref><ref type="bibr">[23]</ref><ref type="bibr">[24]</ref><ref type="bibr">[25]</ref> Lattice-gas models that include surface and subsurface sites can help to narrow the gap in fundamental knowledge by capturing the interplay between adatoms bound above and below the surface of a crystalline solid, ultimately advancing our conceptual understanding of the nature and prevalence of subsurface adsorption in surface chemistry.</p><p>To this end, we have developed a lattice-gas adsorption model of atomic oxygen on Ag(111) that includes surface and subsurface sites and is fully parametrized using DFT. The application to O/Ag(111) was motivated by evidence from several experimental and ab initio computational studies that strongly indicate the formation of subsurface oxygen in Ag(111) and its participation in surface reconstruction, oxide formation, and oxidation catalysis. <ref type="bibr">18,</ref><ref type="bibr">[26]</ref><ref type="bibr">[27]</ref><ref type="bibr">[28]</ref><ref type="bibr">[29]</ref><ref type="bibr">[30]</ref><ref type="bibr">[31]</ref><ref type="bibr">[32]</ref><ref type="bibr">[33]</ref> To our knowledge, this is the first lattice-gas adsorption model that describes both surface and subsurface adsorption. In its present form, the model describes oxygen adsorption at the most strongly binding sites but can be extended to include all types of surface and subsurface sites on Ag(111). The effects of coverage are described using pairwise interactions between oxygen atoms coadsorbed on different sites of the metal. Using this model in combination with Monte Carlo simulations, we computed equilibrium population distributions of oxygen on the surface and subsurface of Ag(111) as a function of the total coverage. On the basis of the computed population distributions, we constructed a DFT model of a total coverage of 0.5 monolayer (ML) to study the differences in the electronic properties of atomic oxygen at the surface and subsurface. We note that because the current form of the model contains oxygenoxygen pair interactions on a fixed silver lattice, it is not able to explicitly describe any type of surface reconstruction or oxide formation. <ref type="bibr">34,</ref><ref type="bibr">35</ref> However, the formation of subsurface oxygen is influenced by interactions with oxygen atoms that are either integrated in a reconstructed surface or are chemisorbed on unreconstructed domains of the surface. Understanding the effects of such interactions on the formation, thermal stability, and chemical nature of subsurface oxygen is the central goal of this work. Overall, our study demonstrates that the lattice-gas adsorption model provides a simple and transferable theoretical framework for exploring the competition between surface and subsurface adsorption in gas-surface systems.</p></div>
<div xmlns="http://www.tei-c.org/ns/1.0"><head n="2.">METHODS</head><p>2.1. Lattice-Gas Adsorption Model. In our lattice-gas model, the total adsorption energy (E total ) of N oxygen atoms on the Ag(111) surface and subsurface is defined as</p><p>where E i (E j ) is the adsorption energy of the ith (jth) oxygen atom in the absence of any neighboring adsorbates and E ij is the total adsorption energy of the ith and jth oxygen atoms when both are bound to Ag(111). Therefore, &#916;E ij is the pairinteraction energy between the ith and jth adsorbed oxygen atoms. E i is negative (positive) for stable (unstable) adsorption and &#916;E ij is negative (positive) for attractive or cooperative (repulsive or uncooperative) interactions. The values of E i and &#916;E ij were calculated using the DFT method described in Section 2.2 with the cutoff distance for &#916;E ij set to a 2 , where a is the lattice constant.</p><p>The (111) surface of a face-centered cubic (fcc) crystal contains four types of high-symmetry sites on the surface, top, bridge, fcc hollow, and hexagonal close-packed (hcp) hollow, and three types of high-symmetry sites in the subsurface, tetrahedral, inverted-tetrahedral, and octahedral. Images and descriptions of the sites are presented in Figure <ref type="figure">S1</ref> and Section 1 of the Supporting Information (SI). As we penetrate deeper into the solid, the inverted-tetrahedral site becomes the same as the tetrahedral site, such that there are two unique types of sites inside the bulk of silver metal, tetrahedral and octahedral.</p><p>In the present study, we included fcc hollow, hcp hollow, and octahedral sites in the lattice-gas model because our DFT results (Section 3.1) showed that these types of sites adsorb atomic oxygen most strongly on the surface, subsurface, and bulk of Ag(111). On the basis of the DFT-computed adsorption energies reported in Table <ref type="table">1</ref>, Table <ref type="table">S1</ref> shows that the chosen types of sites have much greater adsorption probabilities relative to other types of sites in the respective regions. Furthermore, we assessed the strength of three-body interactions (&#916;E ijk = E ijk -E i -E j -E k ) by calculating several cases of &#916;E ijk of three coadsorbed oxygen atoms on Ag(111) using DFT. On the basis of the results (Table <ref type="table">S3</ref>), which show that three-body interactions contribute up to a maximum of &#8764;2.3% of the total adsorption energy, we concluded that manybody interactions between three or more coadsorbed oxygen atoms make a minor contribution to the total adsorption energy, and we did not include them in our lattice-gas model of O/Ag(111).</p><p>2.2. Density Functional Theory Calculations. We performed DFT calculations using the revised Perdew, Burke, and Ernzerhof with Pade&#769;approximation (RPBE) <ref type="bibr">[36]</ref><ref type="bibr">[37]</ref><ref type="bibr">[38]</ref><ref type="bibr">[39]</ref> exchange-correlation functional, as implemented within the Vienna ab initio Simulation Package (VASP). <ref type="bibr">[40]</ref><ref type="bibr">[41]</ref><ref type="bibr">[42]</ref><ref type="bibr">[43]</ref>  The Journal of Physical Chemistry C cutoff of 400 eV was used for the plane-wave basis set, and interactions between nuclei and core electrons were described using the projector augmented wave method. <ref type="bibr">44,</ref><ref type="bibr">45</ref> The Methfessel-Paxton <ref type="bibr">46</ref> smearing method was applied with a smearing width of 0.2 eV, giving an entropy change of &lt;0.01 meV/atom. The RPBE functional, combined with a 15 &#215; 15 &#215; 15 Monkhorst-Pack 47 k-point grid, gave a bulk lattice constant (a) of 4.213 &#197; for Ag, which is 3.0% greater than the corresponding experimental value of 4.09 &#197;. <ref type="bibr">48</ref> The Ag(111) surface was modeled using a p(3 &#215; 4) supercell with six layers of Ag atoms and a vacuum layer twice as thick as the surface slab along the z direction, as shown in Figure <ref type="figure">1</ref>. The dimensions of the supercell were 8.93 &#197; &#215; 10.32 &#197; &#215; 43.78 &#197;. A 5 &#215; 5 &#215; 1 Monkhorst-Pack 47 k-point grid was used to sample the first Brillouin zone. The bottom three layers of the supercell were fixed at the bulk crystal geometry, whereas the top three layers and all adsorbates were allowed to relax in geometry optimizations until the forces on the mobile nuclei were &lt;0.05 eV/&#197;. The adsorption energy, E i , of an oxygen atom at a surface or subsurface site on Ag(111) is defined as</p><p>where E O/Ag(111) is the energy of the relaxed O/Ag(111) supercell, E Ag(111) is the energy of the relaxed Ag(111) supercell, and E O is the energy of an isolated oxygen atom placed inside the empty supercell. Adsorption energies of atomic oxygen were calculated at the high-symmetry sites on the surface, first subsurface (between the first and second silver layers), and second subsurface (between the second and third silver layers) of Ag(111). A negative adsorption energy indicates stable adsorption. If the supercell contains N adsorbed oxygen atoms, then E i is defined as an average over N atoms</p><p>The fractional oxygen coverage (&#952;) is defined as N N Ag,surface , where N Ag,surface is the number of Ag atoms in the topmost layer of the surface. Because the supercell has 12 Ag atoms in each layer, the lowest investigated coverage was &#8776; 0.08 1  12   ML. The electronic energies were corrected for the dipole formed due to charge transfer from the surface to the adsorbates. Charges on Ag and O atoms were calculated using Bader's method of charge density partitioning. <ref type="bibr">[49]</ref><ref type="bibr">[50]</ref><ref type="bibr">[51]</ref> The relative shifts in binding energies of the 1s core electrons of adsorbed oxygen atoms were calculated by exciting the O 1s electron into the valence band using the method implemented in VASP. <ref type="bibr">52,</ref><ref type="bibr">53</ref> The shift values were averaged over the oxygen atoms to study the effects of adsorbate coverage on the core-electron binding energies.</p><p>Absorption energies of oxygen atoms inside the bulk of Ag were calculated using the same p(3 &#215; 4) supercell but excluding the vacuum layer. A 5 &#215; 5 &#215; 5 Monkhorst-Pack kpoint grid was used to sample the first Brillouin zone. Geometry optimizations performed to calculate absorption energies allowed all Ag and O atoms in the supercell to relax. Absorption energies were calculated using eq 3, in which E O/Ag(111) , E Ag(111) , and E O were computed using the bulk supercell.</p><p>2.3. Calculation of Model Parameters Using Density Functional Theory. As shown in eq 1, our lattice-gas model contains two types of terms: the adsorption energy (E i ) of an oxygen atom independent of any coadsorbed oxygen atoms and a pair interaction energy (&#916;E ij ) between the oxygen and any coadsorbed oxygen within a cutoff distance of a 2 . The E i values were determined from DFT-computed adsorption energies, as shown in eq 3, at the lowest coverage afforded by our supercell (&#952; = 1 12 ML). Pair interaction energies were calculated by computing the total adsorption energy (E ij ) of two coadsorbed oxygen atoms using DFT</p><p>and subtracting the E i values of the two atoms from E ij , in accordance with eq 2. The nearest-neighbor distance between two oxygen atoms occupying like sites on the Ag(111) surface is a</p></div>
<div xmlns="http://www.tei-c.org/ns/1.0"><head>2</head><p>, and the cutoff distance of a 2 includes up to thirdneighbor interactions for oxygen pairs occupying like or unlike sites in the surface or subsurface. E i and E ij parameters for absorption and pair interactions inside the solid were calculated using the bulk supercell. The E i parameters are presented in Table <ref type="table">1</ref>, and the values of the &#916;E ij parameters are shown in Table <ref type="table">S2</ref>.</p><p>Figure <ref type="figure">S2</ref> shows a parity plot comparing adsorption energies computed using DFT to the corresponding adsorption energies calculated using the lattice-gas model in the coverage range of 1 4 ML to <ref type="bibr">3 4</ref> ML. The linear-fit equation of the parity plot of y = 0.990x -0.172 with R 2 = 0.998 shows that the model well reproduces the DFT-computed adsorption energies. The inset of the plot shows the parity plot for the subset of oxygen coverages from 1 2 ML to <ref type="bibr">3 4</ref> ML, demonstrating that the approximation of pairwise interactions remains sufficiently accurate at the higher coverages investigated in The Journal of Physical Chemistry C this study. Table <ref type="table">S4</ref> shows the DFT-computed and corresponding model-derived adsorption energies used to create the plot in Figure <ref type="figure">S2</ref>.</p><p>2.4. Monte Carlo Simulations. We performed Monte Carlo (MC) simulations using the lattice-gas model to calculate equilibrium distributions of atomic oxygen on the surface, subsurface, and bulk regions of Ag(111) as functions of the total oxygen coverage and surface temperature. The Metropolis MC algorithm <ref type="bibr">54</ref> was applied within the canonical ensemble to sample oxygen atoms over fcc hollow, hcp hollow, and octahedral sites in the solid, keeping the positions of the Ag atoms fixed. A hexagonal Ag(111) supercell was used, as shown in Figure <ref type="figure">2</ref>, consisting of 10 layers stacked along the z axis of the supercell, each containing 1024 Ag atoms. The hexagonal periodic boundary conditions relate opposite faces of the hexagon in the xy Cartesian plane by symmetry. Octahedral sites in the fourth subsurface (between the fourth and fifth surface layers) and deeper regions of the MC supercell were treated as octahedral sites in the bulk crystal. As described in Sections 2.2 and 2.3, E i and &#916;E ij parameters for adsorption sites on the surface down to the second subsurface region and absorption sites in the bulk of the solid (i.e., fourth subsurface and below) were explicitly calculated using DFT. Adsorption and pair interaction energies in the third subsurface (between the third and fourth surface layers) region were calculated as the averages between corresponding values in the second subsurface and the bulk solid. The choice to transition parameter values from the second surface to the bulk smoothly over two surface layers was motivated by two findings: first, that the parameter values in the second subsurface and the bulk were close to but different from each other, and, second, that the energetics of the inverted tetrahedral site were still different from that of the tetrahedral site in the second subsurface.</p><p>An initial system configuration was obtained by randomly assigning a defined number of oxygen atoms to adsorption sites in the surface and first three subsurface regions. Each move in the MC simulation randomly chose an adsorbed oxygen atom and moved it to a random, available adsorption site anywhere in the 10-layer supercell. Ensemble averages were obtained by using 10 independently initialized Markov chains, each consisting of 1 000 001 configurations that are utilized for ensemble averages.</p></div>
<div xmlns="http://www.tei-c.org/ns/1.0"><head n="3.">RESULTS AND DISCUSSION</head><p>3.1. Adsorption at Low Coverage. Table <ref type="table">1</ref> shows the DFT-computed adsorption properties of atomic oxygen at high-symmetry sites on the surface and subsurface of Ag(111). These properties include the adsorption energy, Bader charge on oxygen, and perpendicular oxygen-surface distance at &#952; = 1 12 ML on top, bridge, hcp hollow, fcc hollow, tetrahedral, inverted tetrahedral, and octahedral sites. Whereas past ab initio studies have computed the adsorption energies of atomic oxygen in the first subsurface of Ag(111), <ref type="bibr">17,</ref><ref type="bibr">27,</ref><ref type="bibr">28,</ref><ref type="bibr">[55]</ref><ref type="bibr">[56]</ref><ref type="bibr">[57]</ref><ref type="bibr">[58]</ref><ref type="bibr">[59]</ref> our study is the first to compute adsorption in the second subsurface and bulk using DFT and incorporate the results in a lattice-gas model. The depth of oxygen atoms shown in Table <ref type="table">1</ref> suggests that oxygen in the second subsurface could be detected in surface-sensitive experiments, such as X-ray photoelectron spectroscopy (XPS). <ref type="bibr">32</ref> We also show the absorption energies of atomic oxygen in octahedral and tetrahedral sites within the bulk of the Ag crystal.</p><p>O binds most strongly to fcc hollow sites on the surface and to octahedral sites in the subsurface. This site preference agrees with previous studies reported in the literature. <ref type="bibr">27,</ref><ref type="bibr">60</ref> The inverted tetrahedral site in the first subsurface is not included in the table, because at this coverage, the oxygen atom rises to the surface and nestles in the hcp hollow site. Overall, our results show that at a low coverage, O binds more strongly to the surface than to the subsurface. Bader charges show that 0.7 electron to 1.1 electron is transferred from the metal surface to the adsorbed oxygen. This charge transfer has a strong effect on the adsorption behavior at high coverages, as described in Section 3.2.</p><p>3.2. Effects of Coverage on Adsorption Energy. Figure <ref type="figure">3</ref> displays the adsorption energy of atomic oxygen on Ag(111) as a function of the total coverage (&#952; total ) at fcc hollow, hcp hollow, octahedral (first subsurface), and octahedral (second subsurface) sites, calculated using DFT (points) and the lattice-gas adsorption model (lines). A multipanel version of  the plot is shown in Figure <ref type="figure">S3</ref>. For each coverage shown in the figure, the oxygen atoms are placed as far apart as possible from each other in the supercell without approaching oxygen atoms in replica supercells. As a result, oxygen pair interactions occur across second-neighbor and third-neighbor distances at &#952; &lt; total 1</p><p>3 ML. However, at &#952; &gt; total 1</p><p>3 ML, first-neighbor pair interactions begin to accumulate and subsequently dominate the effects of oxygen coverage on adsorption energy. The effects of first-neighbor pair interactions on the adsorption energy are discussed later in this section. The model was applied to the same arrangement of oxygen atoms as that in the DFT calculation for each surface and subsurface coverage shown in Figure <ref type="figure">3</ref>. We emphasize that each coverage shown in the plot consists of oxygen atoms occupying a single type of site; therefore, these adsorption energies do not reflect the effects of the adsorption at one type of site on the adsorption at a different type of site. They also do not reflect the effects of the surface adsorption on the subsurface adsorption or vice versa.</p><p>The results show that whereas the adsorption weakens with coverage at fcc hollow and hcp hollow sites, it strengthens with coverage at octahedral sites in the second subsurface. Octahedral sites in the first subsurface show an intermediate trend in which adsorption weakens with coverage for &#952; &lt; total 1 3</p><p>ML and subsequently strengthens with &#952; for &#952; &gt; total 1</p><p>3 ML. Because of opposite trends between the surface and subsurface, atomic oxygen becomes more strongly bound to the subsurface than to the surface for &#952; &gt; total 1</p><p>2 ML. The adsorption energy of -3.38 eV from DFT for a coverage of <ref type="bibr">1 3</ref> ML on fcc hollow sites is in good agreement with the value of -3.44 eV extracted by Campbell from temperature-programmed desorption (TPD) measurements at 0.375 ML. <ref type="bibr">61</ref> Our lattice-gas model (lines) well reproduces the DFT-computed adsorption energies (points) at each of the sites and provides qualitative insight into how differences in interadsorbate interactions lead to a contrast in adsorption behavior between the surface and subsurface.</p><p>We note that the gradient of adsorption energy with respect to coverage changes abruptly at &#952; = total 1</p><p>3 ML for each type of site shown in Figure <ref type="figure">3</ref>. This change occurs because coverages of greater than <ref type="bibr">1 3</ref> ML at a single type of site include at least one first-neighbor interaction due to the symmetry of the Ag(111) lattice. The pairwise interaction energies from our model show that whereas the first-neighbor interaction is strongly repulsive between oxygen atoms at fcc hollow sites (+0.532 eV) and hcp hollow sites (+0.493 eV), it is weakly attractive between oxygen atoms at octahedral sites in the first subsurface (-0.076 eV) and moderately attractive between oxygen atoms at octahedral sites in the second subsurface (-0.181 eV). Figure <ref type="figure">S4</ref> illustrates the increase in adsorption strength due to the increase in the number of first-neighbor interactions between oxygen atoms at octahedral sites in the second subsurface. First-neighbor interadsorbate interactions are very unfavorable on the surface due to negative charges on the oxygen atoms (Table <ref type="table">1</ref>). Second-neighbor interactions, which contribute the most to adsorption energies at low coverage, are repulsive at fcc hollow sites (+0.035 eV), hcp hollow sites (+0.043 eV), and octahedral sites in the first subsurface (+0.092 eV) but attractive at octahedral sites in the second subsurface (-0.034 eV), thus explaining the variations in adsorption energy at coverage of less than 1 3 ML. Adsorption energies at the hcp hollow site vary similarly with coverage as at the fcc hollow site; however, adsorption at the bridge and top sites is not stable at higher coverages. In both the first and second subsurface regions, adsorption energies at the tetrahedral and inverse tetrahedral sites vary similarly with coverage as at the octahedral site. Furthermore, we found that the adsorption of atomic oxygen to the inverted tetrahedral site in the first subsurface becomes stable at higher coverage.</p><p>3.3. Equilibrium Oxygen Distributions from Monte Carlo Simulations. Whereas DFT calculations become expensive with increasing coverage, especially when atomic oxygen occupies more than one type of site, our lattice-gas model can efficiently calculate the adsorption energy of any coverage distributed over different types of sites included in the model. As noted in Section 2.1, the present form of our model includes the most strongly binding sites on the surface, subsurface, and bulk, which are fcc hollow, hcp hollow, and octahedral. Using this model with canonical MC simulations, we calculated the equilibrium population distributions of atomic oxygen on Ag(111) as a function of coverage at room temperature (300 K), shown in Figure <ref type="figure">4</ref>. A key difference The Journal of Physical Chemistry C between Figures <ref type="figure">3</ref> and<ref type="figure">4</ref> is that an oxygen coverage in Figure <ref type="figure">3</ref> was created by binding oxygen atoms to a single type of site, whereas the MC simulations in Figure <ref type="figure">4</ref> allowed oxygen to bind to any combination of site types included in the model. Therefore, in contrast with the adsorption energies shown in Figure <ref type="figure">3</ref>, the results of the MC simulations include the effects of adsorption at one type of site or at one type of region (surface, subsurface, or bulk) on the adsorption at a different type of site or region.</p><p>Figure <ref type="figure">4a</ref> shows that at thermal equilibrium at room temperature, oxygen exclusively populates surface sites up to a saturation limit of 0.375 ML, and further oxygen populates the subsurface. In the range of 0.375 to 0.5 ML, the excess oxygen uniformly occupies octahedral sites in the first and second subsurfaces until the first subsurface saturates at &#8764;0.05 ML; then, oxygen accumulates in the second subsurface between 0.5 to 1 ML. The saturation of surface coverage results from strong interoxygen repulsion, especially between first neighbors on the surface. Figure <ref type="figure">4b</ref> shows that surface oxygen is uniformly distributed between fcc hollow and hcp hollow sites at all studied coverages. Our computed surface saturation coverage of 0.375 ML agrees well with surface-oxygen concentrations of common oxygen-induced reconstructions of Ag(111), such as p(4 &#215; 4) (&#952; = 0.375 ML), p &#215; (4 5 3) (&#952; = 0.375 ML), c &#215; (3 5 3)(&#952; = 0.4 ML), and c(4 &#215; 8) (&#952; = 0.5 ML). This similarity suggests that interoxygen interactions have a strong effect on surface reconstructions of silver.</p><p>The equilibrium oxygen distribution at 300 K (Figure <ref type="figure">5</ref>) shows that whereas oxygen atoms form hexagonal patterns on the surface to minimize repulsive first-neighbor interactions at fcc hollow and hcp hollow sites, they cluster to form islands in the second subsurface to maximize attractive first-neighbor interactions. The low saturation coverage of &#8764;0.05 ML in the first subsurface results from repulsive first-, second-, and thirdneighbor interactions between the oxygen at fcc hollow or hcp hollow sites on the surface and the oxygen at octahedral sites in the first subsurface. This result appears to contradict the results in Figure <ref type="figure">3</ref>, which show that oxygen adsorbs more strongly to the first subsurface than to the second subsurface up to 1 3 ML.</p><p>However, interadsorbate interactions between the surface and first subsurface are absent in Figure <ref type="figure">3</ref> because oxygen atoms occupy a single type of site for each coverage shown in the plot. In Figures <ref type="figure">4</ref> and<ref type="figure">5</ref>, the combination of multiple repulsive interoxygen interactions between the surface and first subsurface, along with a moderately attractive first-neighbor interaction within the first subsurface (&#916;E ij = -0.076 eV), makes it less favorable for oxygen to populate the first subsurface. In contrast, oxygen atoms occupying octahedral sites in the second subsurface have more attractive firstneighbor interactions with each other (&#916;E ij = -0.181 eV) and either weakly attractive or weakly repulsive interactions with oxygen atoms in the first subsurface and the third subsurface. Consequently, atomic oxygen accumulates more in the second subsurface than in the first subsurface in the MC simulations.</p><p>A notable qualitative result of the simulations is that oxygen does not percolate down to the bulk (i.e., fourth subsurface and below in the MC supercell) for total coverages up to 1 ML, despite interoxygen interactions being similar across adjacent subsurface regions from the second subsurface down to the bulk. This stems from a less attractive interaction between oxygen atoms occupying first-neighbor octahedral sites in the bulk (&#916;E ij = -0.142 eV) compared with oxygen atoms occupying first-neighbor octahedral sites in the second subsurface (&#916;E ij = -0.181 eV). These results suggest that oxygen adsorbed in the subsurface is distinct from oxygen dissolved in the bulk of silver for coverages up to 1 ML. In the range of 1 to 2 ML, oxygen atoms start to occupy the third subsurface (first-neighbor &#916;E ij = -0.157 eV), as shown in Figure <ref type="figure">S6</ref>. Finally, we note that the population distributions of a fixed oxygen coverage do not vary in the temperature range of 300 to 600 K because adsorption-energy differences between the surface, first subsurface, and second subsurface are considerably greater than the thermal energy k B T (0.026 to 0.052 eV), that is, the product of the Boltzmann constant (k B ) To summarize the results of our lattice-gas model implemented via Monte Carlo simulations, oxygen atoms adsorb to the surface of Ag(111) up to a saturation coverage of 0.375 ML on the surface, following which any excess oxygen up to a total coverage of 1 ML accumulates in the subsurface, in particular, in the second subsurface. Importantly, the equilibrium oxygen distributions from MC simulations offer compelling computational evidence that the increase in total oxygen coverage on Ag(111) observed by the Killelea group, <ref type="bibr">22,</ref><ref type="bibr">32</ref> despite a fixed maximum surface coverage of 0.375 ML, occurs due to oxygen adsorption in the subsurface. The population distributions also show that subsurface oxygen will emerge to the surface if surface oxygen is depleted, such as in a catalytic oxidation reaction, thereby demonstrating that the subsurface acts as a dynamic reservoir of oxygen.</p><p>3.4. Electronic and Bonding Properties of Surface Oxygen and Subsurface Oxygen. To probe the differences in electronic properties between adsorbates on the surface and adsorbates in the subsurface, we calculated the core-electron binding energies (CBEs) of oxygen atoms at fcc hollow sites on the surface, octahedral sites in the first subsurface, and octahedral sites in the second subsurface of Ag(111) as a function of coverage using DFT, as shown in Figure <ref type="figure">6</ref>. The CBE of an adsorbed oxygen atom increases when its 1s electrons are less shielded from the nucleus by the 2s and 2p electrons, indicating that the valence electrons are more strongly shared with neighboring nuclei. The coverages and oxygen arrangements in Figure <ref type="figure">6</ref> are the same as those in Figure <ref type="figure">3</ref>. The CBE values are plotted as shifts with respect to the CBE of &#952; = ML, appearing to approach invariance with coverage. In contrast, the CBE of oxygen in the first and second subsurface regions more or less steadily decreases with coverage over the studied range. Furthermore, the CBE of a total coverage of subsurface oxygen is greater than the CBE of the corresponding total coverage of surface oxygen, creating a gap of up to +0.8 eV at the lowest coverage, whereas the CBEs of oxygen atoms in the two subsurface regions are similar to each other at all coverages. As the total coverage increases and more oxygen atoms share the same number of Ag atoms, the individual O-Ag interactions weaken, resulting in greater shielding of the core electrons by the valence electrons and consequently lower CBE values.</p><p>Because each coverage studied in Figures <ref type="figure">3</ref> and<ref type="figure">6</ref> was constructed by populating a single type of binding site, the CBE shifts in Figure <ref type="figure">6</ref> do not include the effects of oxygen atoms coadsorbed on different types of sites within the same surface or subsurface region or across two regions. Therefore, we calculated the CBEs of oxygen atoms in the supercell shown in Figure <ref type="figure">7</ref> using DFT. The total coverage in the supercell was chosen to be 0.5 ML, divided into 1 3 ML on the surface, <ref type="bibr">1 12</ref> ML in the first subsurface, and <ref type="bibr">1 12</ref> ML in the second subsurface in accordance with the equilibrium population distributions for &#952; = 0.5 ML in Figure <ref type="figure">4</ref>. The surface oxygen atoms were equally distributed between fcc hollow and hcp hollow sites, as predicted by the model, and the subsurface oxygen atoms were adsorbed to octahedral sites. We note that the average surface adsorption energy at 1 3 ML using this geometry but without any subsurface oxygen was calculated to be -3.40 eV from both DFT and the model, in close agreement with the measured value of -3.44 eV at 0.375 ML by Campbell. <ref type="bibr">61</ref> The calculated shifts with respect to the CBE of &#952; = 1 12 ML at the fcc hollow site, shown in Figure <ref type="figure">8</ref>, reveal significant qualitative differences between the electronic properties of coadsorbed surface and subsurface oxygen. First, the CBE of surface oxygen is similar to the corresponding value in the absence of subsurface oxygen, indicating that subsurface oxygen has a small impact on the CBE of surface oxygen. Second, the CBE of subsurface oxygen undergoes a large change in the presence of surface oxygen, such that the CBE gap between subsurface and surface oxygen increases from +0.45 eV, when present independently, to +0.85 eV when present together. This change suggests that surface oxygen has a considerable effect on the CBE of subsurface oxygen,  ML in the second subsurface of Ag(111). The surface oxygen atoms are uniformly distributed on the fcc hollow and hcp hollow sites, and the subsurface oxygen atoms are bound to octahedral sites. The color scheme is the same as that used in Figure <ref type="figure">1</ref>.</p><p>The Journal of Physical Chemistry C resulting in significantly different electronic properties of atomic oxygen adsorbed to the two regions. The calculated CBE gap of +0.85 eV is in very good agreement with the CBE gap between the low-coverage and high-coverage oxygen species in multiple XPS studies. <ref type="bibr">32,</ref><ref type="bibr">33,</ref><ref type="bibr">57</ref> In these studies, the high-coverage peak at greater CBE was attributed to subsurface oxygen, which is congruent with our results. Finally, the CBEs of oxygen in the two subsurface regions remain similar to each other in the presence or absence of surface oxygen. We note that the qualitative results are robust with a change in the arrangement of the oxygen atoms in each region. For example, when the oxygen atom in the second subsurface was separated from the oxygen atom in the first subsurface by a thirdneighbor instead of the first-neighbor interaction shown in Figure <ref type="figure">7</ref>, the CBE shift of oxygen on the surface, first subsurface, and second subsurface changed to +0.40, +1.45, and +1.33 eV, respectively.</p><p>In the final step of our study, we calculated the projected density of states (PDOS) of the adsorbed oxygen atoms and their neighboring silver atoms from the optimized geometry shown in Figure <ref type="figure">7</ref> to investigate the origin of the CBE difference between surface and subsurface oxygen. Figure <ref type="figure">9a</ref> shows the PDOS of the valence orbitals of an oxygen adsorbed to an fcc hollow site (dotted lines) and the 4d orbitals of a silver atom below it in the first surface layer (solid lines). The silver atom occupies one of the vertices of the three-fold hollow below the oxygen. Figure <ref type="figure">9b</ref> shows the valence orbitals of the oxygen and the 5s-5p orbitals of the same silver atom. The spin-down electrons have the same PDOS as the spin-up electrons and are therefore omitted for clarity. The energies of the states are reported with respect to the Fermi energy of the system. The results show that the energy width of the valence oxygen orbitals in the oxygen-silver bonding region ranges from approximately -5.65 to -2.5 eV and contains a negligible contribution from the 2s orbital. The antibonding region begins at approximately -2.5 eV and extends beyond 0 eV into the unoccupied states of the band. The 4d orbitals of Ag extend from -6.0 to -2.5 eV in the bonding region and make a negligible contribution to the antibonding region, whereas the 5s orbitals extend from -6.5 to -2.5 eV in the bonding region and also make a minor contribution to the antibonding region.</p><p>The PDOS plots in Figure <ref type="figure">9</ref> show a broad overlap between O 2p orbitals and Ag 4d orbitals and a weak overlap between O 2p orbitals and Ag 5s orbitals in the bonding region. The analogous PDOS plot in Figure <ref type="figure">10</ref> for an oxygen atom in the octahedral site of the first subsurface and a silver atom below it in the second layer shows that O 2p orbitals drop to lower energy in the bonding region, ranging from approximately -6.25 to -2.5 eV. The Ag 4d orbitals also drop lower to -6.25 eV, whereas the Ag 5s orbitals remain within an energy range of -6.5 to -2.5 eV, similar to that of the Ag atom in the first surface layer. The O 2p orbitals show a relatively narrow overlap with the left edge of the Ag 4d orbitals and a stronger overlap with the Ag 5s orbital compared with the O-Ag interaction on the surface.</p><p>The O 2p orbitals and Ag 5s orbitals in the second subsurface (Figure <ref type="figure">11</ref>) drop to an even lower energy in the bonding region, starting from -6.5 to -6.8 eV, respectively, whereas the Ag 4d orbitals remain in the same energy range as the Ag 4d orbitals in the first subsurface (Figure <ref type="figure">10</ref>). These changes lead to weaker overlap between O 2p and Ag 4d but stronger overlap between O 2p and Ag 5s compared with the first subsurface.</p><p>Together, the PDOS plots show that the O and Ag orbitals are energetically more stable in the subsurface regions compared with the surface region. Furthermore, the energy alignment between the Ag and O orbitals changes between the surface and subsurface, because of which the O 2p orbitals primarily overlap with Ag 4d orbitals on the surface but with  The Journal of Physical Chemistry C both Ag 4d and 5s orbitals in the subsurface. The overlap between O 2p and Ag 5s is greater in the second subsurface than in the first subsurface. The qualitative differences in O-Ag bonding between the surface and the subsurface result in a significantly greater CBE of subsurface oxygen compared with that of surface oxygen. Importantly, the increasing overlap between O 2p and Ag 5s in the subsurface regions screens the oxygen from neighboring oxygen atoms and decreases the interoxygen electrostatic repulsion, resulting in attractive or cooperative first-neighbor O-O interactions in the subsurface regions. We note that the CBE and PDOS of other surface oxygen atoms from Figure <ref type="figure">7</ref> are qualitatively similar to the plots shown in Figures <ref type="figure">8</ref> and<ref type="figure">9</ref>, and the same trend in orbital stabilization and alignment was found when we considered silver atoms above instead of below the subsurface oxygen atoms. The combination of CBE and PDOS results highlights the qualitative differences in the electronic properties and chemical bonding of atomic oxygen between the surface and subsurface of Ag(111).</p></div>
<div xmlns="http://www.tei-c.org/ns/1.0"><head n="4.">CONCLUSIONS</head><p>In summary, we have developed a lattice-gas model that includes both surface and subsurface adsorption sites in a crystalline solid and describes the effects of adsorbate coverage using pairwise interactions between coadsorbates. We have parametrized the model using DFT to study the competition between the surface and subsurface adsorption of atomic oxygen on the Ag(111) surface as a function of increasing coverage. Canonical Monte Carlo simulations using the model show that at thermal equilibrium at room temperature, oxygen occupies the subsurface of Ag(111) at total coverages of &gt;0.375 ML and prefers to occupy the second subsurface rather than the first subsurface at coverages in the range of 0.5 to 1 ML. The stronger overlap between the 2p orbitals of oxygen and the 5s orbital of silver in the subsurface screens the oxygen atoms from each other more effectively, resulting in less interoxygen repulsion and the consequent accumulation of oxygen in the subsurface. This result is manifested in the model via the transformation in first-neighbor O-O interaction energies from strongly repulsive on the surface, to weakly attractive in the first subsurface, to moderately attractive in the second subsurface. The model shows that the subsurface of the metal acts as a reservoir of oxygen, and the qualitative differences in oxygen-silver bonding between the surface and the subsurface suggest that the two species could play distinct roles in the surface chemistry. In future work, we will apply the model to perform grand canonical Monte Carlo simulations to explore the effects of the surface temperature and the oxygen desorption on the total oxygen coverage and the interplay between the surface and subsurface adsorption for all high-symmetry adsorption sites. Further, we will add explicit oxygen-silver and silver-silver pairwise interactions to the model to study the effects of subsurface oxygen on surface reconstruction.</p></div>
<div xmlns="http://www.tei-c.org/ns/1.0"><head>&#9632; ASSOCIATED CONTENT</head></div>
<div xmlns="http://www.tei-c.org/ns/1.0"><head>* s&#305; Supporting Information</head><p>The Supporting Information is available free of charge at <ref type="url">https://pubs.acs.org/doi/10.1021/acs.jpcc.2c00342</ref>.</p><p>Images of surface sites, adsorption probabilities, parameters of the adsorption model, assessment of three-body interactions, parity plots and data tables  The Journal of Physical Chemistry C</p><p>comparing the model to DFT, multipanel representation of Figure <ref type="figure">3</ref>, illustration of cooperative interactions between oxygen atoms adsorbed in the second subsurface, and plots from MC simulations (PDF) Lattice vectors, final geometries, and energies from DFT geometry optimizations performed at coverages of &#952; = total 1</p><p>12 ML for atomic oxygen on the high-symmetry surface sites, subsurface sites, and bulk sites of Ag(111) and final geometries and energies from DFT geometry optimizations performed to obtain pair-interaction parameters of the lattice-gas model (XYZ)</p></div>
<div xmlns="http://www.tei-c.org/ns/1.0"><head>&#9632; AUTHOR INFORMATION Corresponding Author</head><p>Sharani Roy -Department of Chemistry, University of Tennessee, Knoxville, Tennessee 37996, United States; orcid.org/0000-0002-1020-482X; Email: sharani.roy@ utk.edu</p></div><note xmlns="http://www.tei-c.org/ns/1.0" place="foot" xml:id="foot_0"><p>https://doi.org/10.1021/acs.jpcc.2c00342 J. Phys. Chem. C 2022, 126, 5343-5353</p></note>
		</body>
		</text>
</TEI>
