Two-Dimensional Forms of Robust CO$_2$ Reduction Photocatalysts

Novel photoelectrocatalysts that use sunlight to power the CO$_2$ reduction reaction will be crucial for carbon-neutral power and energy-efficient industrial processes. Scalable photoelectrocatalysts must satisfy a stringent set of criteria, such as stability under operating conditions, product selectivity, and efficient light absorption. Two-dimensional materials can offer high specific surface area, tunability, and potential for heterostructuring, providing a fresh landscape of candidate catalysts. From a set of promising bulk CO$_2$ reduction photoelectrocatalysts, we screen for candidate mono-layers of these materials, then study their catalytic feasibility and suitability. For stable monolayer candidates, we verify the presence of visible-light band gaps, check that band edges can support CO$_2$ reduction, determine exciton binding energies, and compute surface reactivity. We find for SiAs, ZnTe, and ZnSe monolayers, visible light absorption is possible, and reaction selectivity biases towards CO production. We thus identify SiAs, ZnTe, and ZnSe monolayers as targets for further investigation, expanding the chemical space for CO$_2$ photoreduction.


INTRODUCTION
Efficient, stable, scalable photoelectrocatalysts (PECs) which convert sunlight and CO 2 into useful products provide a desirable path towards achieving society's urgent carbon-neutral energy goals 1,2 . Three example applications of CO 2 reduction products include (i) short-term storage of solar energy using methane 3 , which forms a basis for decentralized solar electricity generation, (ii) generating syngas mixtures of CO and H 2 as feedstocks for the Fischer-Tropsch process 4 , or (iii) decreasing the carbon footprint of current industrial processes through efficient production of widely used feedstocks like formic acid. Efficient electrochemical reduction of CO 2 requires catalysts which can survive a strongly reducing environment and provide product selectivity at low overpotentials with respect to the complete reaction pathway 5 . Competing pathways for alternate reactions-for CO 2 reduction, most commonly hydrogen evolution 6,7 -further complicates this picture. Finally, photoelectrocatalysts must clear all the same hurdles while still efficiently capturing light and providing photoexcited electrons at the appropriate CO 2 reduction energy. The search for CO 2 reduction PECs is thus an active and challenging area of research.
In this work, we present a comprehensive study of twodimensional forms of recently suggested candidate CO 2 reduction catalysts. The chosen chemical systems were recently identified for their desirable properties in the bulk phase, such as stability under reducing conditions, suitable band structure, and appropriate band edges 8 . Specifically, we target novel chemistries for CO 2 photoelectroreduction: the compounds ZnTe, ZnSe, GaTe, GaSe, AlSb, SiAs, YbTe, and AlAs form the basis of our investigation. For these compounds, we explore the 2D structural landscape for new thermodynamically and dynamically stable monolayer structures, and evaluate their optoelectronic and reactive suitability as CO 2 reduction PECs.
Photoelectrocatalysis is distinguished from photocatalysis and electrocatalysis by the mechanism of electron excitation and transfer. In photocatalysis, a molecule's potential energy surface is modified by adsorption onto the catalyst, allowing a photon to directly interact with the molecule and induce the desired reaction. Electrocatalysis is characterized by the application of an electric potential U on an adsorbing electrode, tuning the electron chemical potential and causing a desired reaction on the surface to become spontaneous and kinetically facile 9 . The focus of this study is on cathodic photoelectrocatalysis, in which an electron within the catalyst is excited by light which then transfers to an adsorbate and facilitates a reaction 10,11 . In other words, the energy source powering the reaction comes from the incident light instead of an applied bias voltage.
The suite of desired properties for a photoelectrocatalyst-the ability to absorb visible light, allowing reaction intermediates to bind, and generating photoexcited electrons of sufficient energy to trigger a reaction of interest-fundamentally arise from the structure of a material [12][13][14][15][16] . Thus, unexplored structures may yield undiscovered functionality. Two-dimensional (2D) systems, characterized by a layered crystalline structure, offer a less well-studied and highly tunable avenue for future PECs 17,18 . The surprising effects of two-dimensional material interactions, such as ZnS/PbO and ZnSe/PbO heterostructures exhibiting enhanced solar absorption 19 , expands the space of possible 2D-material based reactors even further beyond monolayers to the massive combinatoric space of heterostructures and 2D material coatings 17,20,21 . addresses the feasibility of candidate 2D structures of a given material, and the suitability of these structures as PECs. We study feasibility by estimating the thermodynamic and dynamic stability of the structures, and suitability by assessing band gaps, band edges, exciton binding energies and adsorption energies of reaction intermediates.
Previous experimental and theoretical works have investigated the structural, electronic, and chemical characteristics of a few monolayer phases of GaSe 22 , GaTe 23 , ZnSe 19,24-32 , ZnTe 29,30,32 , and SiAs [33][34][35] . 2D layers of ZnSe in various structural forms have attracted attention for their photoabsorption properties 24,36 , as have low-dimensional forms of ZnTe bound to nanowires or in nanoparticle form [37][38][39] . However, to our knowledge, none of the 2D systems we consider have been studied for CO 2 photoelectrocatalyst applications.
Feasibility screening Thermodynamic screening. To "seed" our search for stable monolayers, we drew from the subset of 258 monolayer structures published by Mounet et al. 40 as well as the 2D prototype structures from the Computational 2D Materials Database by Haastrup et al. 41 . We also included the native monolayer form of layered bulk SiAs. For our target materials, we "mapped" elements into monolayer structures of binary compounds (e.g. mapping ZnTe into the MoS 2 or h-BN structures).
A recent data-driven study found the energy difference between the formation energy of the monolayer as compared to the equivalent bulk phase to be one of the most important predictors for MAX and MXene monolayer stability 42 . Here, we define the change in formation energy, ΔE F , as where E F,Monolayer is the formation energy of the 2D material and E F,Bulk is the formation energy of its most stable 3D bulk counterpart. For e.g. ZnTe, we used the structure mp-2176 in the Materials Project database 43 to define the bulk formation energy, and normalize each formation energy per atom in the corresponding unit cell. Fig. 2 shows the formation energies of the various 2D structures that were considered for the compositions AlAs, AlSb, YbTe, ZnSe, ZnTe, GaSe, GaTe, and SiAs. All formation energies were computed using density-functional theory with the SCAN+rVV10 functional (see the "Methods" section for more details).
We use a small ΔE F ≤ 200 meV/atom stability cutoff for our candidate monolayer structures, which Singh et al. 18 and Haastrup et al. 41 identified as a useful heuristic stability criteria. The structures lying in the green shaded region of Fig. 2 satisfy this criteria, and Fig. 3 depicts the stable structural prototypes. The ΔE F ≤ 200 meV/atom criteria resulted in eight different structural prototypes for five of the compounds, and excludes three compounds: YbTe, AlSb, and AlAs. We found only one candidate structure for ZnSe, but three for ZnTe and SiAs, four for GaTe, and five for GaSe. The naming conventions for all eight structure types are shown in Fig. 3. Later in the study, we will focus on ZnTe and ZnSe in structural prototypes that resemble the CuI and the CuBr structures 44 . For notational simplicity, we refer to them as the 'hexagonal' and 'tetragonal' structures.
Four structural prototypes (the GaSe, GaS, InSe, and native SiAs structures in Fig. 3) feature chemical environments in which the anions are coordinated by three cations, and each cation by three anions and another cation. The bulk structures of GaSe, SiAs, and GaTe all feature exactly this coordination, explaining their compatibility with most of the monolayer structures with the same coordination. SiAs in the InSe structure, which features the same coordination pattern, comes close to the cutoff at 216 meV/ atom above the bulk.
ZnTe and ZnSe in the bulk crystallize in the zincblende structure because they are 'octet compounds' 45 which attain chemical equilibrium by filling an eight-electron set of s and p valence orbitals. Additionally, the hexagonal and tetragonal prototypes can be understood respectively as a cleaved (110) plane and (100) plane of the bulk zincblende structure 31 . Unlike the GaSe and GaS structures, the tetragonal and hexagonal structure types feature fourfold opposite species coordination, explaining their ability to support ZnTe. Note that monolayer ZnSe meets the thermodynamic cutoff only in one configuration, the hexagonal structure. Previous studies 25,27,31,32 have found that ZnSe in the tetragonal structure is dynamically stable, though it did not meet our thermodynamic stability criterion and hence we did not include it in the next steps of our study; we found it to have ΔE F = 250 meV/ atom above the bulk using SCAN+rVV10. It was found (See Table 1 of Zhou et al. 31 ) using the PBE 46 exchange-correlation Fig. 1 Overview of workflow aimed at the following assessments for candidate structures. a Thermodynamic stability, as measured by comparing the formation energy of a monolayer against a cutoff. b Dynamic stability, as indicated by non-imaginary phonon frequencies in vacuum. c Solar photoreducibility, indicated by visible light band gaps and appropriate band edges. d Reactivity, as measured by the binding energy of adsorbates (COOH, CO), which are indicative of a CO 2 reduction pathway to form carbon monoxide. Fig. 2 Difference in formation energy per atom between a candidate 2D structure and the ground-state bulk, as computed using the SCAN+rVV10 functional. On the horizontal are bulk compounds we attempted to find 2D phases of. The horizontal offset of points within a column is for visual clarity. On the vertical are differences in energy. Structures which are below our threshold are marked with symbols annotated in Fig. 3, and those above it, with an X.
functional that the tetragonal structure form had a higher energy per atom than the hexagonal form, which is consistent with our findings, but Li et al. found using the HSE06 functional 47,48 a lower energy for the tetragonal structure 25 .
Low-dimensional GaSe 23,49,50 , GaTe 23 , and ZnSe 24 have all been experimentally studied with structures similar to those which we found. We are not aware of any experimental evidence of twodimensional SiAs or ZnTe, though both were featured in theoretical studies 29,30,33,34 . However, previous computational efforts considered different structural prototypes for ZnTe than in this work.
Dynamical stability screening. For all the structures which passed the thermodynamic screening, we sought to determine their dynamic stability. The full phonon band structures evidenced dynamic stability in vacuum, with most structures yielding nonimaginary frequencies. In many cases, we noted small, imaginary acoustic phonon modes close to the zone center; however, in the majority of cases, these disappeared upon increasing the size of the simulation cell. For ZnTe, ZnSe, and one GaTe structure, despite very large supercells (up to 8 × 8 × 1), we retained small instabilities in the elastic limit (See Supplementary Tables 1 and 2 for calculation details, and Supplementary Figs. 1-16 for full band structures). For these structures, perturbing along the negative displacement vector and then relaxing yielded in all cases a return to the starting structure, suggesting that the unstable elastic modes are artifacts of an insufficient supercell size and the known difficulty of fitting the quadratic Γ z-direction acoustic phonon 40 . GaSe in the GaS + T-Type structure type evidenced broad imaginary phonon modes off the Brillouin zone center, indicating a dynamically unstable structure.
Suitability screening HSE band gaps and band edges. We computed the electronic structure properties for all structures which passed the stability screening in order to determine band gaps and edges, which are relevant to catalytic application. The HSE06 hybrid exchangecorrelation functional 47,48 is known to exhibit an improved treatment of semiconductor bandgaps: in Fig. 4, we show that the band gaps of the 2D materials computed using HSE06 lie mostly within the visible light spectrum, and the band edges are appropriate for CO 2 reduction (see Supplementary Table 3 for values used to compute the work function, and Supplementary Figs. 17-46 for full band structure and density of states). We found the GaS + T-Type structure of GaSe was metallic.
To the best of our knowledge, the electronic structure of most of the 2D structural prototypes in this study have not been explored computationally before. Tong et al. 32 using PBE found ZnTe in the tetragonal form to have a band gap of 0.88 eV, lower than our HSE06 gap of about 1.6 eV. This is to be expected, as the PBE functional is known to underestimate band gaps. We reproduce Bai's monolayer SiAs indirect Y-Γ HSE gap of 2.39 eV 34 , though Cheng 33 reports a direct gap of 2.353 eV. Li et al. 25 found a band gap of 3.4 eV for the same ZnSe structure using HSE06, about 500 meV higher than what we find.
Quasiparticle gaps and exciton binding energies. We calculated the quasiparticle energies within the GW approximation for the self energy which captures the many-body effects in the electronic structure 51,52 of materials. To calculate the optical response of materials, we solve the Bethe-Salpeter equation (BSE) 53-55 , where we explicitly take the electron-hole interaction into account. The solution of the BSE can be used to compute the imaginary part of the dielectric function, which enables us to study the excitonic effects in the absorption spectrum. Low exciton binding energy is desirable for photovoltaic applications, as it enables easier electron-hole separation resulting in a higher device efficiency. However, in two-dimensional materials such as TMDCs, oftentimes large exciton energies (~1 eV) 56 have been observed, arising from large effective charge carrier mass, strong Coulomb interactions, and weak dielectric screening 57 among other factors.
In Table 1 Fig. 1. Parentheses below each panel indicate the materials which met the thermodynamic criteria for the given structure, and a symbol that is used to refer to that structure in Figs. 2 and 5. Naming conventions for structures are from Ashton et al. 44 ; we use the term "compressed SiAs structure'' as it was the result of relaxation from the SiAs structure for ZnTe. Compounds pictured by column are: first column, SiAs (blue/green). Second column, ZnTe (silver/gold). Third column, GaSe (Ga large and green/Se small and light green). Fourth column, on top, SiAs, on bottom, ZnTe.
absorption spectra). However, for a few of them we find low EBEs, such as 0.43 eV for ZnTe in the compressed SiAs-structure, 0.55 eV for SiAs in the native SiAs-structure, and 0.58 eV for GaTe in the GaS-structure.
The GW-BSE optical band gaps predicted for hexagonal ZnTe are in the visible light range (2.32 eV); those for SiAs in its native structure are slightly lower (1.33 eV) and for hexagonal ZnSe, slightly higher (3.43 eV). Notably, transition metal chalcogenides can experience dramatic changes in opto-electronic behavior in the single-layer limit 58 . For example, the photocurrent density of monolayer ZnSe increases dramatically 24 compared to the bulk. However, among the materials studied here, all but one of the 2D structures are predicted to be semiconductors, similar to their bulk counterparts. The HSE band gaps tended to increase for for SiAs, ZnSe, and GaSe, and were slightly lower for ZnTe 8 . Additionally, we find that ZnSe and ZnTe are, like their bulk parent structures, direct band-gap visible light semiconductors 8 , which suggests promising potential for efficient light absorption. Monolayer GaSe, as in its bulk form, exhibits a visible-light indirect band gap 8 . In the bulk, GaSe and SiAs have been identified as a material of interest for photovoltaic applications 49,59 . On the other hand, GaSe and GaTe were recently reported to exhibit a sharp decrease in photoluminesence during the transition from the many to fewlayer limit 23 .
Adsorption and reactivity. Finally, we studied adsorbate binding on the surfaces, which allows us to gauge the surface reactivity of the candidate photoelectrocatalyst monolayers. CO 2 reduction can proceed along several complex reaction pathways 5,60-62 . We examined the reaction pathway which predicts reactivity towards methane on copper 5 ; though, in finding that CO does not bind to the surfaces, we focus only on the steps in Eqs. (1) and (2) below, where * represents an arbitrary surface.
From adsorbate binding energies, we can compute the theoretical overpotential, which estimates the bias voltage that must be applied to the electrode in order for the reaction pathway to occur entirely downhill in free energy; this is equal to the greatest change in free energy. Note that under photoelectrocatalyst operating conditions, energy which facilitates the reaction would be provided by the photoexcited electron instead of the reaction being made spontaneous by an applied bias voltage. In this study we examine reactivity 'in the dark' to gain insights on reaction selectivity and binding propensity; the full mechanisms of photoexcited electron charge transfer from the surface are beyond the scope of this study. The studied surfaces were the pristine basal planes of the candidate materials, with no defects nor edges. Figure 5 presents the resulting binding energies and   Tables 5-7 for calculation parameters). We focused on ZnTe, ZnSe, and SiAs for this phase of our study, because of the experimental evidence that monolayer GaTe and GaSe absorb light poorly in the monolayer limit 23 . Examination of the bonding behavior of COOH and CO on the surfaces showed interesting results. For four structural prototypes (SiAs: GaSe, and GaS, ZnSe: Hexagonal, and ZnTe: Hexagonal) we found that COOH binds via inducing a surface reconstruction of a Si or Zn atom, as the COOH radical 'pulls' the cation out-of-plane. This behavior can be seen in the inset plots of Fig. 5. We rationalize this for Zn as the Zn atom maintaining fourfold coordination by bonding with the adsorbate over the Se/Te atom below it; we also find this induces the Se/Te atom directly below the binding Zn to shift down and out-of-plane on the bottom. Interestingly, for the native SiAs structure, the most favorable bonding environment is COOH adsorbed to the As atom highest out-of-plane, serving as a lower energy configuration than bonding to any Si atom.
For every prototype we examined, we found no evidence of CO bonding to the surface; during relaxations of CO instantiated on various sites, CO uniformly shifted to far away from the surface. The 'binding energies' of these configurations were approximately near zero. We examined the hydrocarbon-forming pathway from CO 2 to CHO 5 , but the lack of CO binding effectively terminates it at the second step. Selectivity can arise from several mechanisms 63,64 , but our results show promise for CO production selectivity, and select materials could be combined with a cocatalyst to facilitate other catalytic processes of interest.
For tetragonal ZnTe, we found that COOH did not bind to the surface, and the reconstruction seen in the hexagonal structure in which the Zn atom 'pops' out-of-plane was not observed. This suggests reactivity may be improved by the presence of vacancy or adatom defects which expose an undercoordinated Zn or Si atom. For instance, planar vacancy defects in transition metal chalcogenides have been found to reduce the binding energies of *COOH and *OH and reduce the limiting potentials for *COOH and *CO production (in some cases, increasing the density of vacancies changed binding energy by more than -0.3 eV) 65 .
While some 2D materials like TaS 2 and NbS 2 have been experimentally found to exhibit high basal plane reactivity 66 , the pristine structures we examined exhibit relatively unreactive overpotentials in the range of 1.61-2.02 eV, with SiAs in the GaSe structure around 2.85 eV. Due to their unreactivty, catalytic efficacy from these surfaces would have to originate from deviations from the pristine structure which could help reactivity, e.g. surface defects. The well known two-dimensional catalyst MoS 2 has an inert basal plane which can greatly increase in reactivity through defect engineering 67,68 . However, these conclusions would need to be carefully qualified. One common pitfall for CO 2 reduction catalysts is that because the CO 2 reduction reaction is driven by reducing potentials, hydrogen evolution presents a constant competitor.
The lack of CO binding indicates a possibility for selectivity towards CO production. This conclusion would also demand close study: weakly bound CO alone is not always solely responsible for favorable CO production. For instance, in Au, careful mechanistic studies revealed interfacial field effects as facilitating CO production at cathodic potentials 63,64 . The prediction of weak CO binding suggests that SiAs, ZnTe, and ZnSe belong to the class of materials characterized by the weak-binding leg of the volcano relationships derived in e.g. Kuhl et al 69 .
The central assumption of the computational hydrogen electrode (CHE) model 70 is that the electron chemical potential at which a given reaction step becomes downhill corresponds to a threshold of fast kinetics for that step, i.e. that the kinetic barrier between the two states of the step is consistent between materials. CHE is widely applied to metals, and a comparison of these kinetic barriers with semiconductors may be tentative, given that the electron transfer barriers that are normally very fast relative to the chemical kinetics in the proton-transfer portion of the proton-coupled electron transfer (PCET) may not be so in the case of a semiconductor interface. While it is true that semiconductors will not typically exhibit excited conduction-band electrons to contribute to a reaction, we assume that the experimental operation of these photoelectrocatalysts would occur under illumination. In this case, photoexcited electrons should be available for transfer to the reactant. Regardless, PCET reactions often represent an upper bound on the chemical activity.
It is an open question as to whether the intermediate kinetics on a semiconductor surface tend to resemble that of a metal surface. However, our primary interpretation of the results presented here are that these materials-in their pristine stateare not reactive, and so this distinction would not substantially affect our conclusions.

DISCUSSION
In conclusion, our study examined the structural, thermodynamic, dynamic, electronic, and reactive properties of 2D forms of CO 2 reduction PECs. The primary contributions of this work are (i) the uncommon chemistries of these materials in the CO 2 reduction literature, (ii) our thorough exploration of possible monolayer phases, and (iii) prediction of weak CO binding and therefore possible CO selectivity of SiAs, ZnSe, and ZnTe monolayers.
In native SiAs, hexagonal ZnTe, and hexagonal ZnSe, the comparatively lower overpotentials in tandem with ZnTe and ZnSe's direct band gaps present targets for further theoretical investigation, experimental verification, and property optimization. Possible means of engineering the catalytic activity of these structures includes tuning their reactivity and optical properties via defects, such as vacancies, adatoms, or dopants 71 . Heterostructuring could simultaneously allow for tuned opto-electronic properties or induced reactivity (as has been done for a tetragonal form of ZnSe with PbO) 19 . In particular, the predicted weak CO binding for SiAs, ZnSe, and ZnTe could facilitate selectivity, which would demand the presence of a co-catalyst to facilitate further reactions. Further excited-state studies could probe possible photoexcitation and reaction pathways 72 . For ZnTe and ZnSe in particular, the projected density of states indicate that the high valence bands tend to exhibit anionic (Se, Te) character and the low conduction bands have the character of the cation, Zn, similar to behavior seen in oxide systems 73 . Closer examination of the interactions between individual adsorbates and these states will be the subject of further study.
In conclusion, novel 2D chemistries and structures for CO 2 photoreduction are suggested and explored computationally. The structures we studied in this paper preserve their bulk counterparts' stability and semiconductor properties. Further enhancement of surface reactivity remains a challenge, and further work is needed to understand how chemical changes, defect engineering and surface treatments can be used to influence and tune the performance.

General calculation details
We performed all first-principles calculations in the Vienna Ab-Initio Simulation Package with the Projector Augmented-Wave method [74][75][76][77] . Phonon calculations were assisted by Phonopy 78 with high cutoff energies at and above 700 eV, and supercells ranging in size from 5 × 5 × 1 to 8 × 8 × 1. Automated workflows for relaxation, band gap estimation, and adsorption 79 were performed using Atomate 80 , using standard Materials Project parameters 43,81 with small modifications specified in the SI. Structure matching, mapping, and general calculation IO operations were performed with pymatgen 82 and FireWorks 83 . Work function analysis was performed using pymatgen's Surface Analyzer package 84,85 . More details on the calculation parameters, reference energies, and analysis are available in the SI.
For the energy calculation depicted in Fig. 2, we first relax structures with the PBE functional 46 , then PBE+DFT-D3 86 , and again with SCAN +rVV10 87 , a Meta-GGA functional with a van der Waals correction. In a small set of test cases, we found that the PBE+DFT-D3 functional tended to optimize structures closer to the final SCAN+rVV10 structures at a fraction of the computational cost. While this procedure introduced slightly more planning overhead to the production-scale calculations, it helped to reduce the time spent running SCAN+rVV10 calculations. We compare the energies with bulk structures using SCAN+rVV10 because it has shown good performance for estimating exfoliation energies of layered structures 88 . Due to its improved prediction of band gaps 89 , we relaxed the structures and computed the electronic structure using the HSE06 hybrid functional 47,48 . Exciton calculations were performed in VASP.

Mapping structures
In order to map the prototype structures from our two 2D material databases, (those of Mounet et al. 40 and Haastrup et al. 41 ) to the materials in our study, we performed the following procedure.
We used the 258-structure subset of the Mounet database which each contain six or fewer atoms per unit cell. We also considered all of the structural prototypes released by the Computational Two-Dimensional Materials Database with AB stoichiometry.
Iterating through each binary compound in the Singh photocatalyst database 8 , for every structure prototype which had two elements (e.g. h-BN, MoS 2 ), we mapped the elements from the material into the prototype, including the opposite mapping (i.e. ZnTe × MoS 2 → ZnTe 2 , Zn 2 Te). Since some prototypes like h-BN have AB sublattice symmetry, we also used the pymatgen structure matcher to prune repeats from the resulting set of structures this generated. For the Mounet database, this generated 36 structures; for the Haastrup database, 24. However, the structure matcher failed to catch some which were sufficiently different by angle, so there were repeats within and between the sets.

Stability assessment
The workflow for computing the formation energy of each structure is as follows.
Because we generate structures by substituting the atoms in structures with new elements, the resultant initial candidates deviate from their equilibrium structure and require relaxation. We performed the relaxation in three steps (Each to within force convergence of ≤0.01 eV/Å) for all of the generated two-dimensional structures, as well as their ground-state bulk phases as defined in the Materials Project database.
1. First, we use the PBE Functional with the MPRelax VASP input set from pymatgen, using the POT_GGA_PAW_PBE pseudopotentials, which use the Projector Augmented Wave method 77 . We used a kpoint density of 64 points per cubic angstrom (and manually set the number of k-points along the z axis to 1). 2. Second, we use the DFT-D3 functional of Grimme et al 86 . We do this in order to save time when using a more computationally expensive van der Waals functional for our final relaxation. Here we used a kpoint density of 100 per cubic angstrom. 3. Finally, we perform a relaxation using the SCAN+rVV10 functional 87 , which has been shown to have good performance for layered materials. This functional combines MetaGGA accuracy using SCAN 90 with a non-local van der Waals correction 91 . We used a kpoint density of 400 per cubic angstrom, and the PBE_54 VASP pseudopotentials. Once the energies of both the monolayer and bulk phases were known using the same set of pseudopotentials and VASP settings, we were able to compare the formation energy per atom between the two structures on equal footing. This calculation was performed by referencing the formation energy per atom in the bulk cell by comparing to the bulk forms of each element.
For example, for ZnTe, we computed the formation energy per atom of bulk Zn and bulk Te to obtain μ Zn and μ Te . Then, we find the formation energy of bulk ZnTe (which has 2 atoms each of Zn and Te) by subtracting E bulk À 2μ Zn À 2μTe n bulk where n bulk = 4 and the formation energy of a monolayer with i Zn and j Te atoms as E 2D À iμ Zn À jμ Te n 2D S.B. Torrisi et al.
where n 2D = i + j. Then, we can compare the normalized difference ΔE f = E 2D − E bulk . This made it easier to compare the formation energy of structures with different numbers of atoms per unit cell than the bulk.

Phonon calculations
Phonon spectra were computed using Phonopy 78 via the frozen-phonon method in VASP. One sample POSCAR, band.conf, and INCAR file is included in the SI, to demonstrate how they were generated. Calculation details were tweaked for individual materials, but all involved high ENCUT values of 700 + eV, ADDGRID set to True, and LREAL set to False (variable names and conventions appropriate to VASP 5.4.4). A supercell of size 6 × 6 or 7 × 7 was used for the different materials, in some cases increasing the displacement of the atoms to a radius of 0.05 angstrom, up from the default 0.01. This was necessary partially due to the notorious difficulty of converging the quadratic ZA acoustic phonon which exists for most hexagonal 2D materials, which often produced negative frequencies about the Γ point that vanished with increasing calculation precision and supercell size (See Supplementary Table 1). We found evidence that the negative frequencies were numerical artifacts for all but one structure (see caption of Supplementary Table 2). Because the acoustic phonons in the z-direction are quadratic, the perturbations involved in a frozen phonon calculation sample from very small forces, and thus they require minimal numerical noise in order to be sampled accurately. Using high cutoff energies of 700 + eV helped to ensure that the forces were well converged and that small perturbations of the atoms would recover the correct force response. Phonon band structures are presented in Supplementary Figs. 1-16.

HSE band gaps and band edges
Band gaps and band edges were computed using a workflow in the Atomate package for computing HSE band structures. First, the structure is relaxed using HSE (since small variations in lattice constant can have significant effects on the resultant band gap values 92 ). Then, the band gap is computed along an automatically generated k-point path using Pymatgen 93 . The Fermi energy is also extracted from these calculations. The two-dimensional surfaces were all oriented in the x-y plane, so the planar average was computed for each xy plane in the cell. The value of the electric potential in the vacuum was extracted (as the maximum planar average potential along the z-axis), and was used to find the work function ϕ such that ϕ = ϵ F + V vac . We then used the work function as the location of the valence band maximum. The values used to compute this quantity are all in Supplementary Table 3.
We present for each structure a computed band structure and a projected density of states on a per-element and a per-orbital basis. We compute the projected density of states in VASP using the default Wigner radius RWIGS specified in the pseudopotentials for each element in Supplementary Figs. 31-60.

GW-BSE
We used the Vienna Ab-initio Simulation Package (VASP) to perform GW 0 (by self-consistently iterating only G, but not W) and BSE calculations 94,95 . There are several other GW schemes available, such as partial selfconsistent GW and self consistent GW [96][97][98][99] , but due to their computational expense, we restrict ourselves to the GW 0 level of calculation. We include semicore states by using PAW pseudopotentials specifically designed for GW calculations to capture the strong exchange-interaction resulting from those states. The dielectric function was expanded in plane waves with energy up to 250 eV. We included 80 frequency points to perform the fully frequency dependent GW calculations. We have used 600 bands in our GW calculations, which is sufficient to converge the QP energies to <0.1 eV for all the materials reported in this study. In Supplementary Table 4 we report the k-grids used for the BSE calculations of different materials.
We included 10 valence and 10 conduction bands in the BSE calculations, which is sufficient to ensure the convergence of the absorption spectrum for energies up to~6 eV for all materials. In Supplementary Figs. 17-30, we show the absorption spectra of all the materials obtained from BSE calculations. We have shown the QPG and the OPG for each of them by using blue and red vertical lines respectively.

Adsorption analysis
Reference energies. We used the RPBE 100 exchange correlation functional for all adsorption calculations. Parameters such as plane-wave cutoff energy came from the MVLSlab set of VASP input parameters as found in the pymatgen python package. We used a standard dipole correction along the z-axis as implemented in VASP via the IDIPOL flag for calculations involving molecules adsorbed on surfaces.
We performed a DFT calculation for the isolated molecules H 2 O, CO, and H 2 , each in a 10 × 10 × 10 angstrom cell using only the gamma K-point. The total energies of these isolated reference 'molecules in a box' E DFT were used, after free energy and other corrections, as a reference to compute the chemical potentials of hydrogen, carbon, and oxygen atoms as E ref .
The obtained values are the final values in Supplementary Table 5. The reference schemes are made explicit below.
Energy corrections. Entropic and thermodynamic effects appear as terms in the Gibbs free energy, influencing the thermodynamic favorability of certain states. We therefore must compute free energy corrections based on the degrees of freedom of molecules and adsorbates. Because the reaction pathway we study ends at CO (with CO not binding to the surface) we computed the vibrational modes for five molecules in total. They are: isolated H 2 , isolated H 2 O, isolated CO, isolated CO 2 , and COOH adsorbed onto five different surfaces. The surfaces are SiAs in the GaSe, GaS, and native structure, ZnTe in the hexagonal structure, and ZnSe in the hexagonal structure. COOH and CO both failed to bind to ZnTe in the tetragonal structure. We used all 3N = 12 vibrational degrees of freedom for COOH adsorbed on the surface. For the nonlinear molecule H 2 O we used 3N-6 = 3 degrees of freedom. For the linear molecules CO 2 and H 2 , we used 3N-5 = 2 degrees of freedom. Following Jones 101 and Peterson 5 , we use the harmonic approximation for the adsorbed COOH, neglecting the contributions from configurational and rotational entropy for the adsorbed COOH. We treat the gas-phase molecules with full configurational, vibrational, and rotational entropy. The vibrational frequencies used for our free energy correction are detailed in Supplementary Table 5, and the correction values used are in Supplementary Table 6.
The free energy of gas molecules require both temperature and pressure. We use 3534 Pa for H 2 O as the experimental vapor pressure of water at room temperature. For CO, CO 2 , and H 2 , we use standard state of 1 atm (101,325 Pa). Equation 5 of Jones 101 is given as (with minor notational adjustment to our case): Here, we use E Ref as the relaxed DFT-computed energy of a given molecule, E ZPE as the zero-point-energy (the energy contribution from the occupation of the lowest vibrational states), ΔU 0,T as the energy associated with the specific heat of various degrees of freedom, and S as the entropy associated with a given temperature. We use T = 300 K for all cases.
To find vibrational modes, we use used the IBRION = 5 flag in VASP which computes the Hessian matrix associated with small displacements. We neglect the adsorbate-induced change in vibrational state of the surface, and so used selective dynamics to only compute the molecule's vibrational modes. We additionally used NFREE = 4, which computes the dynamical matrix using a four-point finite-difference stencil, and POTIM = 0.015, which is the displacement distance associated with the finite difference calculation in angstrom. VASP computes the eigenfrequencies from the Hessian matrix, which we then analyzed using the Ideal-GasThermo and HarmonicThermo functions of the Thermochemistry package in the atomic simulation environment (ASE) Python library 102 .
Computational details. The Atomate adsorption workflow 79 was used to generate structures with the adsorbate molecules placed on top of highsymmetry sites detected using Delaunay triangulation on the structures of interest. A dipole correction was applied using VASP's built-in features.
We performed adsorption calculations across five compounds and three unique structure prototypes. For SiAs, we used the GaSe, GaS, and native SiAs structure types. For GaSe we used the GaSe structure type, and for GaTe, the GaS structure type. For ZnSe, we used the hexagonal type, and for ZnTe, we used the hexagonal and tetragonal structure type. The adsorbates were placed on top of a 3 × 3 supercell of atoms on high-S.B. Torrisi et al. symmetry sites. There were six high-symmetry sites identified for the hexagonal (GaSe/GaS/Hexagonal) cells with chemical formula AB: atop A, atop B, between A and B, bridge A and the center, bridge B and the center, and above the center of the hexagon. For the one tetragonal structure, the sites were three: atop A, atop B, and above vacuum.
The lowest energy configuration of a given molecule was chosen to use for the binding energy associated with that molecule. The total energies per atom of pristine unit cells were computed for each slab and used to compute the energy of the pristine slabs.
Free energy diagrams. We referenced carbon with CO, hydrogen with H 2 , and oxygen with H 2 O.
For instructive purposes, we show the full justification below for the reaction ðCO 2 Þ g þ Ã þ ðH þ þ e À Þ ! COOHÃ; (7) for which the change in free energy can be understood as (approximating the energetic contribution of a proton as the contribution from a hydrogen atom): ΔGðCOOHÃÞ ¼ GðCOOHÃÞ À EðÃÞ À ðEðH þ Þ þ Eðe À ÞÞ À G CO 2 ð Þ g (8) ¼ GðCOOHÃÞ À EðÃÞ À GðH 2 Þ 2 þ eU À ΔGðCO 2 Þ g À μðCÞ À 2μðOÞ (9) ¼ GðCOOHÃÞ À EðÃÞ À μðHÞ þ eU À ΔGðCO 2 Þ ðgÞ À μðCÞ À 2μðOÞ (10) We use μ(C), μ(H), and μ(O) as the reference energies of individual atoms. Where G(COOH*) comes from the DFT energy of COOH adsorbed onto the slab *, plus the corrections due to the vibrational degrees of freedom of the adsorbate, and a stabilization of an adsorbate with an *OH group of −0.25 eV following Peterson 5 . For G(H 2 ), we include the energy from H 2 in DFT as well as the Gibbs corrections detailed in Supplementary Table 7. We neglect free-energy corrections for the surface atoms and simply use the energy found using DFT.
Because we use CO and H 2 O as our reference for carbon and hydrogen respectively, we obtain the energy of CO 2(g) referenced against them as ΔGðCO 2 Þ ¼ GðCO 2 Þ À GðCOÞ þ ðGðH 2 Þ À GðH 2 OÞÞ; (11) where G includes the DFT energy, the contributions from trans-rovibrational degrees of freedom, and all other aforementioned corrections in Supplementary Table 7. Note that the corrections used lead to a positive energy associated with the formation of CO 2 gas. Our expression for the binding energy referenced against the slab and adsorbed molecule is E B ðCOOHÞ ¼ GðCOOHÃÞ À EðÃÞ À μðCOOHÞ (12) so re-arrangement and substituting E(COOH*) in yields ΔG ¼ E B ðCOOHÞ þ EðÃÞ À EðÃÞ þ μðCOOHÞ À μðCO 2 Þ À μðHÞ (13) ¼ E B ðCOOHÞ þ μðCÞ À μðCÞ þ 2ðμðOÞ À μðOÞÞ þ μðHÞ À μðHÞ þ eU À ΔGðCO 2 Þ ðgÞ (15) ¼ E B ðCOOHÞ þ eU À ΔGðCO 2 Þ ðgÞ : Above, e is the magnitude of the elementary charge (which can be treated as 1), U is the bias voltage applied on the electrode, and * represents the surface, either for the energy of pristine surfaces or to indicate a molecule is adsorbed. Therefore, we get the useful simplification that the differences in free energy between reaction steps can be expressed as the difference in the binding energies for each individual adsorbate.

Reconstruction
In order to verify that the deformation of the surface is not facile (which could suggest instability in the presence of adsorbates), we computed the energy associated with the distortions produced by binding of *COOH to three different structure types in 3 × 3 supercells using the same VASP parameters as for the adsorbate calculation. We also compared the anioncation bond lengths to show that the local chemical environments of the distorted atoms do not substantially change.
The distorted structures and associated energies were hexagonal ZnSe (+0.791 eV), hexagonal ZnTe (+1.261 eV), and GaSe-type SiAs (+3.209 eV). For ZnSe, the pristine Zn-Se bond length is~2.48 Å, and the translated atom has mean nearest-neighbor Se bond length of 2.50 Å. For ZnTe, the pristine Zn-Te bond length is approximately 2.69 Å, and the translated atom has mean nearest-neighbor Te bond length of 2.76 Å. Finally, for SiAs, the pristine Si-As bond length is~2.41 Å, which was unchanged for the translated atom.
Input files and graphics