Introduction

Efficient, stable, scalable photoelectrocatalysts (PECs) which convert sunlight and CO2 into useful products provide a desirable path towards achieving society’s urgent carbon-neutral energy goals1,2. Three example applications of CO2 reduction products include (i) short-term storage of solar energy using methane3, which forms a basis for decentralized solar electricity generation, (ii) generating syngas mixtures of CO and H2 as feedstocks for the Fischer–Tropsch process4, or (iii) decreasing the carbon footprint of current industrial processes through efficient production of widely used feedstocks like formic acid. Efficient electrochemical reduction of CO2 requires catalysts which can survive a strongly reducing environment and provide product selectivity at low overpotentials with respect to the complete reaction pathway5. Competing pathways for alternate reactions—for CO2 reduction, most commonly hydrogen evolution6,7—further complicates this picture. Finally, photoelectrocatalysts must clear all the same hurdles while still efficiently capturing light and providing photo-excited electrons at the appropriate CO2 reduction energy. The search for CO2 reduction PECs is thus an active and challenging area of research.

In this work, we present a comprehensive study of two-dimensional forms of recently suggested candidate CO2 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 edges8. Specifically, we target novel chemistries for CO2 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 CO2 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 facile9. 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 reaction10,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 material12,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 PECs17,18. The surprising effects of two-dimensional material interactions, such as ZnS/PbO and ZnSe/PbO heterostructures exhibiting enhanced solar absorption19, expands the space of possible 2D-material based reactors even further beyond monolayers to the massive combinatoric space of heterostructures and 2D material coatings17,20,21.

Results

Workflow

To probe the diversity of 2D structure types for each compound and investigate the feasibility of CO2 reduction, we devise a computational workflow (see Fig. 1). Conceptually, the workflow 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.

Fig. 1: Overview of workflow aimed at the following assessments for candidate structures.
figure 1

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 CO2 reduction pathway to form carbon monoxide.

Full size image

Previous experimental and theoretical works have investigated the structural, electronic, and chemical characteristics of a few monolayer phases of GaSe22, GaTe23, ZnSe19,24,25,26,27,28,29,30,31,32, ZnTe29,30,32, and SiAs33,34,35. 2D layers of ZnSe in various structural forms have attracted attention for their photoabsorption properties24,36, as have low-dimensional forms of ZnTe bound to nanowires or in nanoparticle form37,38,39. However, to our knowledge, none of the 2D systems we consider have been studied for CO2 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 MoS2 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 stability42. Here, we define the change in formation energy, ΔEF, as

$$\Delta {E}_{{\rm{F}}}={E}_{{\rm{F}},{\rm{Monolayer}}}-{E}_{{\rm{F}},{\rm{Bulk}}}$$

where EF,Monolayer is the formation energy of the 2D material and EF,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 database43 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).

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.
figure 2

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.

Full size image

We use a small ΔEF ≤ 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 ΔEF ≤ 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 structures44. For notational simplicity, we refer to them as the ‘hexagonal’ and ‘tetragonal’ structures.

Fig. 3: Visuals of structures which cleared the thermodynamic screening process depicted in Fig. 1.
figure 3

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.

Full size image

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 structure31. 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 studies25,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 ΔEF = 250 meV/atom above the bulk using SCAN+rVV10. It was found (See Table 1 of Zhou et al.31) using the PBE46 exchange-correlation 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 functional47,48 a lower energy for the tetragonal structure25.

Low-dimensional GaSe23,49,50, GaTe23, and ZnSe24 have all been experimentally studied with structures similar to those which we found. We are not aware of any experimental evidence of two-dimensional SiAs or ZnTe, though both were featured in theoretical studies29,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 non-imaginary 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 phonon40. 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 exchange-correlation functional47,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 CO2 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.

Fig. 4: The HSE0647,48 band gaps tend to lie within the visible light spectrum.
figure 4

a Points indicate band gap values associated with different structures. Boxed points indicate direct band gaps. The y-axis denotes band gap values. b The band edges fulfill the requisite energy for CO2 reduction. Blue (red) denotes valence (conduction) band maxima (minima). On both subplots, the x-axis groups band gaps and edges by material. Reaction energies on right from Fig. 3 in Singh et al8.

Full size image

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 eV34, though Cheng33 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 structure51,52 of materials. To calculate the optical response of materials, we solve the Bethe–Salpeter equation (BSE)53,54,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 screening57 among other factors.

In Table 1 we show the quasiparticle gap (QPG), optical gap (OPG), and exciton binding energy (EBE) obtained from GW-BSE calculations. The QPG values we report here are the minimum direct bandgap for these materials, as we have not included any optical transitions with finite momentum transfer in our BSE calculations. The optical gap values reported in the table are the lowest excitation energies obtained by diagonalizing the BSE Hamiltonian. The EBE is then computed as the difference between the QPG and OPG. Our calculation shows for most of the materials the EBEs are large, similar to TMDCs (See Supplementary Table 4 for K-point grids used, and Supplementary Figs. 47–60 for 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.

Table 1 Exciton calculation results.
Full size table

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 limit58. For example, the photocurrent density of monolayer ZnSe increases dramatically24 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 ZnTe8. Additionally, we find that ZnSe and ZnTe are, like their bulk parent structures, direct band-gap visible light semiconductors8, which suggests promising potential for efficient light absorption. Monolayer GaSe, as in its bulk form, exhibits a visible-light indirect band gap8. In the bulk, GaSe and SiAs have been identified as a material of interest for photovoltaic applications49,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 few-layer limit23.

Adsorption and reactivity

Finally, we studied adsorbate binding on the surfaces, which allows us to gauge the surface reactivity of the candidate photoelectrocatalyst monolayers. CO2 reduction can proceed along several complex reaction pathways5,60,61,62. We examined the reaction pathway which predicts reactivity towards methane on copper5; 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.

$${{\rm{C}}{\rm{O}}}_{2(g)}+(\text{H}^{+}+\text{e}^{-})+* \to ^*\!\!{\rm{COOH}}$$
(1)
$$^*{\rm{COOH}}+(\text{H}^{+}+\text{e}^{-})\to {{\rm{C}}{\rm{O}}}_{(g)}+{{\rm{H}}}_{2}{{\rm{O}}}_{(g)}+*$$
(2)

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 the free energy reaction pathways for ZnTe, ZnSe, and SiAs (see “Methods” section for full details, and Supplementary 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 limit23.

Fig. 5: Free-energy diagrams for COOH bonding on SiAs, ZnTe, and ZnSe present relatively unreactive surfaces.
figure 5

Inset are images from the side (upper left) and at an angle (upper right) of COOH bonding on surfaces; the reconstruction is clearly visible when contrasting the side view with those from Fig. 3. Structure types are from Fig. 3. SiAs in the native structure presents the lowest barrier. For complete computational details, see the “Methods” section.

Full size image

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 CO2 to CHO5, but the lack of CO binding effectively terminates it at the second step. Selectivity can arise from several mechanisms63,64, but our results show promise for CO production selectivity, and select materials could be combined with a co-catalyst 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 TaS2 and NbS2 have been experimentally found to exhibit high basal plane reactivity66, 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 MoS2 has an inert basal plane which can greatly increase in reactivity through defect engineering67,68. However, these conclusions would need to be carefully qualified. One common pitfall for CO2 reduction catalysts is that because the CO2 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 potentials63,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 al69.

The central assumption of the computational hydrogen electrode (CHE) model70 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 state—are 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 CO2 reduction PECs. The primary contributions of this work are (i) the uncommon chemistries of these materials in the CO2 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 dopants71. 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 pathways72. 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 systems73. 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 CO2 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.

Methods

General calculation details

We performed all first-principles calculations in the Vienna Ab-Initio Simulation Package with the Projector Augmented-Wave method74,75,76,77. Phonon calculations were assisted by Phonopy78 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 adsorption79 were performed using Atomate80, using standard Materials Project parameters43,81 with small modifications specified in the SI. Structure matching, mapping, and general calculation IO operations were performed with pymatgen82 and FireWorks83. Work function analysis was performed using pymatgen’s Surface Analyzer package84,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 functional46, then PBE+DFT-D386, and again with SCAN+rVV1087, 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 structures88. Due to its improved prediction of band gaps89, we relaxed the structures and computed the electronic structure using the HSE06 hybrid functional47,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 database8, for every structure prototype which had two elements (e.g. h-BN, MoS2), we mapped the elements from the material into the prototype, including the opposite mapping (i.e. ZnTe × MoS2 → ZnTe2, Zn2Te). 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. 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 method77. We used a k-point density of 64 points per cubic angstrom (and manually set the number of k-points along the z axis to 1).

  2. 2.

    Second, we use the DFT-D3 functional of Grimme et al86. 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 k-point density of 100 per cubic angstrom.

  3. 3.

    Finally, we perform a relaxation using the SCAN+rVV10 functional87, which has been shown to have good performance for layered materials. This functional combines MetaGGA accuracy using SCAN90 with a non-local van der Waals correction91. We used a k-point 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

$$\frac{{E}_{{\rm{bulk}}}-2{\mu }_{{\mathrm{Zn}}}-2\mu {\mathrm{Te}}}{{n}_{{\rm{bulk}}}}$$

where nbulk = 4 and the formation energy of a monolayer with i Zn and j Te atoms as

$$\frac{{E}_{2{\rm{D}}}-i{\mu }_{{\rm{Zn}}}-j{\mu }_{{\rm{Te}}}}{{n}_{2{\rm{D}}}}$$

where n2D = i + j. Then, we can compare the normalized difference ΔEf = E2D − Ebulk. 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 Phonopy78 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 values92). Then, the band gap is computed along an automatically generated k-point path using Pymatgen93. 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 + Vvac.

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 GW0 (by self-consistently iterating only G, but not W) and BSE calculations94,95. There are several other GW schemes available, such as partial self-consistent GW and self consistent GW96,97,98,99, but due to their computational expense, we restrict ourselves to the GW0 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 RPBE100 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 H2O, CO, and H2, 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’ EDFT were used, after free energy and other corrections, as a reference to compute the chemical potentials of hydrogen, carbon, and oxygen atoms as Eref. The obtained values are the final values in Supplementary Table 5. The reference schemes are made explicit below.

$$\mu ({\rm{H}})=\frac{{E}_{{\rm{ref}}}({{\rm{H}}}_{2})}{2}$$
(3)
$$\mu ({\rm{C}})={E}_{{\rm{ref}}}({\rm{C}}{\rm{O}})-{E}_{{\rm{ref}}}({{\rm{H}}}_{2}{\rm{O}})+{E}_{{\rm{ref}}}({{\rm{H}}}_{2})$$
(4)
$$\mu ({\rm{O}})={E}_{{\rm{ref}}}({{\rm{H}}}_{2}{\rm{O}})-{E}_{{\rm{ref}}}({{\rm{H}}}_{2})$$
(5)

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 H2, isolated H2O, isolated CO, isolated CO2, 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 H2O we used 3N–6 = 3 degrees of freedom. For the linear molecules CO2 and H2, we used 3N–5 = 2 degrees of freedom. Following Jones101 and Peterson5, 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 H2O as the experimental vapor pressure of water at room temperature. For CO, CO2, and H2, we use standard state of 1 atm (101,325 Pa). Equation 5 of Jones101 is given as (with minor notational adjustment to our case):

$$G={E}_{{\rm{Ref}}}+{E}_{{\rm{ZPE}}}+\Delta {U}^{0,T}-TS.$$
(6)

Here, we use ERef as the relaxed DFT-computed energy of a given molecule, EZPE as the zero-point-energy (the energy contribution from the occupation of the lowest vibrational states), ΔU0,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 IdealGasThermo and HarmonicThermo functions of the Thermochemistry package in the atomic simulation environment (ASE) Python library102.

Computational details

The Atomate adsorption workflow79 was used to generate structures with the adsorbate molecules placed on top of high-symmetry 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-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 H2, and oxygen with H2O.

For instructive purposes, we show the full justification below for the reaction

$${({\rm{C}}{{\rm{O}}}_{2})}_{g}+* +(\text{H}^{+}+\text{e}^{-})\to {\rm{C}}{\rm{O}}{\rm{O}}{\rm{H}}* ,$$
(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):

$$\Delta G({\rm{C}}{\rm{O}}{\rm{O}}{\rm{H}}* )=G({\rm{C}}{\rm{O}}{\rm{O}}{\rm{H}}* )-E(*)-(E({\text{H}}^{+})+E({\text{e}}^{-}))-G{\left({\rm{C}}{{\rm{O}}}_{2}\right)}_{g}$$
(8)
$$=G({\rm{C}}{\rm{O}}{\rm{O}}{\rm{H}}* )-E(* )-\frac{G({{\rm{H}}}_{2})}{2}+eU-\Delta G{({\rm{C}}{{\rm{O}}}_{2})}_{g}-\mu ({\rm{C}})-2\mu ({\rm{O}})$$
(9)
$$=G({\rm{C}}{\rm{O}}{\rm{O}}{\rm{H}}* )-E(* )-\mu ({\rm{H}})+eU-\Delta G{({\rm{C}}{{\rm{O}}}_{2})}_{(g)}-\mu ({\rm{C}})-2\mu ({\rm{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 Peterson5. For G(H2), we include the energy from H2 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 H2O as our reference for carbon and hydrogen respectively, we obtain the energy of CO2(g) referenced against them as

$$\Delta G({\rm{C}}{{\rm{O}}}_{2})=G({\rm{C}}{{\rm{O}}}_{2})-G({\rm{C}}{\rm{O}})+(G({{\rm{H}}}_{2})-G({{\rm{H}}}_{2}{\rm{O}})),$$
(11)

where G includes the DFT energy, the contributions from trans-ro-vibrational 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 CO2 gas.

Our expression for the binding energy referenced against the slab and adsorbed molecule is

$${E}_{{\rm{B}}}({\rm{C}}{\rm{O}}{\rm{O}}{\rm{H}})=G({\rm{C}}{\rm{O}}{\rm{O}}{\rm{H}}* )-E(* )-\mu ({\rm{C}}{\rm{O}}{\rm{O}}{\rm{H}})$$
(12)

so re-arrangement and substituting E(COOH*) in yields

$$\Delta G={E}_{{\rm{B}}}({\rm{C}}{\rm{O}}{\rm{O}}{\rm{H}})+E(* )-E(* )+\mu ({\rm{C}}{\rm{O}}{\rm{O}}{\rm{H}})-\mu (\rm{C}{{\rm{O}}}_{2})-\mu ({\rm{H}})$$
(13)
$$-\Delta G{({\rm{C}}{{\rm{O}}}_{2})}_{(g)}+eU$$
(14)
$$={E}_{{\rm{B}}}({\rm{C}}{\rm{O}}{\rm{O}}{\rm{H}})+\mu ({\rm{C}})-\mu ({\rm{C}})+2(\mu ({\rm{O}})-\mu ({\rm{O}}))+\mu ({\rm{H}})-\mu ({\rm{H}})+eU-\Delta G{({\rm{C}}{{\rm{O}}}_{2})}_{(g)}$$
(15)
$$={E}_{{\rm{B}}}({\rm{C}}{\rm{O}}{\rm{O}}{\rm{H}})+eU-\Delta G{({\rm{C}}{{\rm{O}}}_{2})}_{(g)}.$$
(16)

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 anion-cation 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

Crystal structure visualization used the VESTA software103. VASP Input files, including POSCAR and INCAR files, are supplied under the Supplementary Methods section in the Supplementary Material.