The Computational 2D Materials Database: high-throughput modeling and discovery of atomically thin crystals

The workflow is initiated with a crystal structure defined by its unit cell (Bravais lattice and atomic basis). The crystal lattice is typically that of an experimentally known prototype (the 'seed structure') decorated with atoms picked from a subset of the periodic table, see figure 2. We refer to materials by the chemical formula of their unit cell followed by the crystal structure. The latter is indicated by a representative material of that prototype, as described in section 3.1. For example, monolayer MoS₂ in the hexagonal H and T phases are denoted MoS₂-MoS₂ and MoS₂-CdI₂, respectively. Now, MoS₂ is in fact not stable in the T phase, but undergoes a $2\times 1$ distortion to the so-called T' phase. Because the T' phase is the thermodynamically stable phase of WTe₂, we denote MoS₂ in the distorted T phase by Mo₂S₄-WTe₂. In the following, we shall refer to the unit cell with which the workflow is initiated, i.e. the unit cell of the seed structure, as the primitive cell or the $1\times 1$ cell, even if this cell is not dynamically stable for the considered material (see section 2.4).

Zoom In Zoom Out Reset image size

The unit cell and internal coordinates of the atoms are relaxed in both a spin-paired (NM), ferromagnetic (FM), and anti-ferromagnetic (AFM) configuration. Calculations for the AFM configuration are performed only for unit cells containing at least two metal atoms. The symmetries of the initial seed structure are kept during relaxation. All relevant computational details are provided in table 1.

After relaxation, we check that the structure has remained a covalently connected 2D material and not disintegrated into 1D or 0D clusters. This is done by defining clusters of atoms using the covalent radius [39]  +  30% as a measure for covalent bonds between atoms. The dimensionality of a cluster is obtained from the scaling of the number of atoms in a cluster upon repetition of the unit cell following the method described by Ashton et al [16]. Only materials containing exactly one cluster of dimensionality 2 are given further consideration (an exception is made for the metal-organic perovskites (prototype PbA₂I₄) for which the metal atom inside the octahedron represents a 0D cluster embedded in a 2D cluster). To illustrate the effect of the covalent radius  +  30% threshold, figure 3 shows the distribution of the candidate structures in the database as a function of the covalent factor needed to fully connect the structure. Most materials have a critical covalent factor below 1.3 and fall in the green shaded region. There is, however, a tail of around 100 disconnected materials (red region); these materials are not included in the database (see first green box in figure 1).

Zoom In Zoom Out Reset image size

We also check that the material is not already contained in the database (second green box in figure 1). This is done by measuring the root mean square distance (RMSD) [40] relative to all other materials in the C2DB with the same reduced chemical formula. A threshold of 0.01 Å is used for this test.

In case of multiple metastable magnetic configurations (in practice, if both a FM and AFM ground state are found), these are regarded as different phases of the same material and will be treated separately throughout the rest of the workflow. To indicate the magnetic phase we add the extensions 'FM' or 'AFM' to the material name. The total energy of the spin-paired ground state is always stored, even when it is not the lowest. If the energy of the non-magnetic state is higher than the most stable magnetic state by less than 10 meV/atom, the workflow is also performed for the non-magnetic state. This is done in recognition of the finite accuracy of DFT for predicting the correct energetic ordering of different magnetic states.

We have compared the lattice constants of 29 monolayers with those reported in [41], which were obtained with the VASP code using PBE and very similar numerical settings and find a mean absolute deviation of 0.024 Å corresponding to 0.4%. The small yet finite deviations are ascribed to differences in the employed PAW potentials.

2.2.1. Symmetry

2.2.2. Prototypes

The materials are classified into crystal structure prototypes based on the symmetry of the crystals. For two materials to belong to the same prototype, we require that they have the same space group, the same stoichiometry, and comparable thicknesses. The last requirement is included to distinguish between materials with the same symmetry and stoichiometry but with different number of atomic layers, see for example monolayer BN and GaS in figure 4. Each prototype is labelled by a specific representative material. For prototypes which have been previously investigated, we comply with the established conventions. However, since the field of 2D materials is still young and because C2DB contains a large number of never-synthesised materials, some of the considered crystal structures fall outside the known prototypes. In these cases we have chosen the representative material to be the one with the lowest energy with respect to the convex hull. Some of the crystal structure prototypes presently contained in the C2DB are shown in figure 4.

Zoom In Zoom Out Reset image size

The heat of formation, $\Delta H$, is defined as the energy of the material with respect to the standard states of its constituent elements. For example, the heat of formation per atom of a binary compound, A_xB_y, is

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \Delta H = (E({{\rm A}}_x{{\rm B}}_y) - xE({{\rm A}}) - yE({{\rm B}}))/(x+y), \nonumber \end{align} \tag{ 1 }$$

where $E({{\rm A}}_x{{\rm B}}_y)$ is the total energy of the material A_xB_y, and $E({{\rm A}})$ and $E({{\rm B}})$ are the total energies of the elements A and B in their standard state. When assessing the stability of a material in the C2DB, it should be kept in mind that the accuracy of the PBE functional for the heat of formation is only around 0.2 eV/atom on average [43]. Other materials databases, e.g. OQMD, Materials Project, and AFLOW, employ fitted elementary reference energies (FERE) [44] and apply a Hubbard U term [45] for the rare earth and transition metal atoms (or a selected subset of them). While such correction schemes in general improve $\Delta H$ they also introduce some ambiguity, e.g. the dataset from which the FERE are determined or the exact form of the orbitals on which the U term is applied. Thus in order not to compromise the transparency and reproducibility of the data we use the pure PBE energies.

For a material to be thermodynamically stable it is necessary but not sufficient that $\Delta H<0$. Indeed, thermodynamic stability requires that $\Delta H$ be negative not only relative to its pure elemental phases but relative to all other competing phases, i.e. its energy must be below the convex hull [46]. We stress, however, that in general, but for 2D materials in particular, this definition cannot be directly applied as a criterion for stability and synthesisability. The most important reasons for this are (i) the intrinsic uncertainty on the DFT energies stemming from the approximate xc-functional (ii) substrate interactions or other external effects that can stabilise the monolayer (iii) kinetic barriers that separate the monolayer from other lower energy phases rendering the monolayer (meta)stable for all practical purposes.

We calculate the energy of the 2D material relative to the convex hull of competing bulk phases, $\Delta H_{{\rm hull}}$. The convex hull is currently constructed from the 2836 most stable binary bulk compounds which were obtained from the OQMD [24]. The energies of the bulk phases were recalculated with GPAW using the PBE xc-functional and the same numerical settings as applied for the 2D materials (but the structure was not re-optimised). Because the bulk reference structures from OQMD were optimised with the VASP code and with Hubbard U corrections for materials containing 3d elements, and because the PBE misses attractive vdW interaction, the bulk energies could be slightly overestimated relative to the monolayers. As a consequence, monolayers that also exist in a layered bulk phase could have $\Delta H_{{\rm hull}}<0$, even if the layered bulk phase is part of the convex hull and thus should be energetically more stable than the monolayer. Comparing our $\Delta H_{\mathrm{hull}}$ values for 35 compounds with the exfoliation energies calculated in [18] employing vdW compliant xc-functionals for both bulk and monolayer, we estimate the errors in the convex hull energies to be below 0.1 eV/atom.

As an example, the convex hull for Fe_xSe_1−x is shown in figure 5. The convex hull as defined by the bulk binaries is indicated by the blue lines. The labels for the 2D materials refer to the crystal prototype and magnetic order. Clearly, most 2D materials lie above the convex hull and are thus predicted to be thermodynamically unstable in freestanding form under standard conditions. However, as mentioned above, depending on the material, errors on the PBE formation energies can be sizable and thus the hull diagram should only be taken as guideline. Nevertheless, in the present example we find that FeSe (which is itself a prototype) with anti-ferromagnetic ordering lies slightly below the convex hull and is thus predicted to be thermodynamically stable. This prediction is consistent with the recent experimental observation that monolayer FeSe deposited on SrTiO₃ exhibits AFM order [47].

Zoom In Zoom Out Reset image size

Due to the applied symmetry constraints and/or the limited size of the unit cell, there is a risk that the structure obtained after relaxation does not represent a local minimum of the potential energy surface, but only a saddle point. We test for dynamical stability by calculating the Γ-point phonons of a $2 \times 2$ repeated cell (without re-optimising the structure) as well as the elastic constants (see section 2.5). These quantities represent second-order derivatives of the total energy with respect to atom displacements and unit cell lengths, respectively, and negative values for either quantity indicate a structural instability.

The Γ-point phonons of the $2 \times 2$ supercell are obtained using the finite displacement method [48]. We displace each atom in the primitive cell by  ±0.01 Å, and calculate the forces induced on all the atoms in the supercell. From the forces we construct the dynamical matrix, which is diagonalised to obtain the Γ-point phonons of the $2\times 2$ cell (or equivalently the Γ-point and zone boundary phonons of the primitive cell). The eigenvalues of the dynamical matrix correspond to the square of the mass-renormalised phonon frequencies, $\tilde{\omega}$. Negative eigenvalues are equivalent to imaginary frequencies and signal a saddle point.

Our procedure explicitly tests for stability against local distortions of periodicities up to $2 \times 2$ and thus provides a necessary, but not sufficient condition for dynamic stability. We stress, however, that even in cases where a material would spontaneously relax into a structure with periodicity larger than $2 \times 2$, the Γ-point dynamical matrix of the $2\times 2$ cell could exhibit negative eigenvalues. Our test is thus more stringent than it might seem at first glance. In principle, a rigorous test for dynamic stability would require the calculation of the full phonon band structure. Mathematically, the instabilities missed by our approach are those that result in imaginary phonons in the interior of the BZ but not at the zone boundary. Physically, such modes could be out of plane buckling or charge density wave-driven reconstructions with periodicities of several unit cells. In general, however, these types of instabilities are typically rather weak (as measured by the magnitude of the imaginary frequency) as compared to more local distortions such as the T to T' distortion considered below. Moreover, they could well be a special property of the isolated monolayer and become stabilised by the ubiquitous interactions of the 2D material with its environment, e.g. substrates. This is in fact supported by the full phonon calculations by Mounet et al for $\sim 250$ isolated monolayers predicted to be easily exfoliable from experimentally known layered bulk phases [18]. Indeed, most of the instabilities revealed by their calculations are of the type described above and would thus be missed by our test. However, these instabilities cannot be too critical as the monolayers are known to be stable in the vdW bonded layered bulk structure.

As an example, figure 6 compares the dynamical stability of a subset of transition metal dichalcogenides and -oxides in the T and T' phases (CdI₂ and WTe₂ prototypes). The two upper panels show the smallest eigenvalue of the Γ-point dynamical matrix of the $2\times 2$ cell. Only materials above the dashed line are considered dynamically stable (for this example we do not consider the sign of the elastic constants which could further reduce the set of dynamically stable materials). Since the unit cell of the T' phase contains that of the T phase it is likely that a material initially set up in the T' phase relaxes back to the T phase. To identify these cases, and thereby avoid the presence of duplicates in the database, the third panel shows the root mean square distance (RMSD) between the structures obtained after relaxations starting in the T- and T' phase, respectively. Structures below the dashed line are considered identical. The color of each symbol refers to the four different potential energy surfaces illustrated at the bottom of the figure.

Zoom In Zoom Out Reset image size

2.4.1. Stability criteria

To assess the stability of the materials in the C2DB, we turn to the set of experimentally synthesised/exfoliated monolayers. For these materials, the calculated energy above the convex hull and minimum eigenvalue of the dynamical matrix are shown on figure 7. It is clear that all but five known monolayers have a hull energy below 0.2 eV/atom, and three of these have only been synthesised on a metal substrate. Turning to the dynamical stability, all but one of the experimentally known monolayers have a minimum eigenvalue of the dynamical matrix above $-2~\mathrm{eV}~\mathring{\rm A}^{-2}$, and 70% have a minimum eigenvalue above $-1\times10^{-5}~\mathrm{eV}~\mathring{\rm A}^{-2}$.

Zoom In Zoom Out Reset image size

Guided by these considerations, we assign each material in the C2DB a stability level (low, medium or high) for both dynamical and thermodynamic stability, as illustrated in table 2. For ease of reference, we also define the overall stability level of a given material as the lower of the dynamical and thermodynamic stability levels. If a material has 'low' overall stability (marked by bold in the table), we consider it unstable and do not carry out the rest of the workflow. Materials with 'high' overall stability are considered likely to be stable and thus potentially synthesisable. Materials in the 'medium' stability category, while unlikely to be stable as freestanding monolayers, cannot be discarded and might be metastable and possible to synthesise under the right conditions. For example, free-standing silicene has a heat of formation of 0.66 eV/atom, but can be grown on a silver substrate. Likewise, the T' phase of MoS₂ (WTe₂ prototype) has an energy of 0.27 eV/atom higher than the thermodynamically stable H phase, but can be stabilised by electron doping.

Table 2. The materials in the C2DB distributed over the nine stability categories defined by the three levels (high, medium and low) of dynamical stability (columns) and thermodynamic stability (rows). The overall stability of the materials is defined as the lower of the two separate stability scales. Materials with low overall stability (bold) are considered unstable.

Thermodynamic stability (eV/atom)	Dynamic stability (eVÅ−2)			Total
	$\left| \tilde \omega_{\mathrm{min}}^2\right| > 2$ or Cii  <  0	$10^{-5} < \left|\tilde\omega_{\mathrm{min}}^{2}\right| < 2$, Cii  >  0	$\quad \left|\tilde\omega_{\mathrm{min}}^2\right| < 10^{-5}$, $\ \ C_{ii} > 0$
$\Delta H > 0.2$	6.0%	4.2%	1.7%	12.0%
$\Delta H < 0.2$	14.9%	10.9%	6.4%	32.2%
$\Delta H_\mathrm{hull} < 0.2$	11.4%	24.1%	20.3%	55.8%
Total	32.3%	39.2%	28.5%

The elastic constants of a material are defined by the generalised Hooke's law,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:hooke} \sigma_{ij} = C_{ijkl} \varepsilon_{kl} \nonumber \end{align} \tag{ 2 }$$

where $\sigma_{ij}$, C_ijkl and $\varepsilon_{kl}$ are the stress, stiffness and strain tensors, respectively, and where we have used the Einstein summation convention. In two dimensions, the stress and strain tensors have three independent components, namely planar stress/strain in the x and y directions, as well as shear stress/strain. The stiffness tensor is a symmetric linear map between these two tensors, and therefore has up to six independent components. Disregarding shear deformations, the relationship between planar strain and stress is

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \left[\begin{array}{@{}ccc@{}} \sigma_{xx} \nonumber \\ \sigma_{yy} \end{array} \right] = \left[\begin{array}{@{}ccc@{}} C_{11} & C_{12} \nonumber \\ C_{12} & C_{22} \end{array} \right] \left[\begin{array}{@{}ccc@{}} \epsilon_{xx} \nonumber \\ \epsilon_{yy} \end{array} \right]. \nonumber \end{align} \tag{ 3 }$$

For all materials in the C2DB, we calculate the planar elastic stiffness coefficients C₁₁, C₂₂, and C₁₂. These are calculated using a central difference approximation to the derivative of the stress tensor: the material is strained along one of the coordinate axes, x or y, and the stress tensor is calculated after the ions have relaxed. We use strains of $\pm$$1\%$ which we have found to be sufficiently large to eliminate effects of numerical noise and sufficiently small to stay within the linear response regime.

Table 3 shows the calculated planar stiffness coefficients of a set of 2D materials. As can be seen the values from the C2DB are in very good agreement with previously published PBE results. For the isotropic materials MoS₂, WSe₂ and WS₂, C₁₁ and C₂₂ should be identical, and we see a variation of up to 0.6%. This provides a test of how well converged the values are with respect to numerical settings.

The energy dependence on the direction of magnetisation, or magnetic anisotropy (MA), arises from spin–orbit coupling (SOC). According to the magnetic force theorem [96] this can be evaluated from the eigenvalue differences such that the correction to the energy becomes

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:anisotropy} \Delta E(\hat{\mathbf n}) = \sum_{{\bf k}n}f(\varepsilon_{{\bf k}n}^{\hat{\mathbf{ n}}})\varepsilon_{{\bf k}n}^{\hat {\mathbf n}}-\sum_{{\bf k}n}f(\varepsilon_{{\bf k}n}^0)\varepsilon_{{\bf k}n}^0, \nonumber \end{align} \tag{ 4 }$$

where $\varepsilon_{{{\bf k}}n}^{\hat{\mathbf n}}$ and $f(\varepsilon_{{\bf k}n}^{\hat {\mathbf n}})$ are the eigenenergies and occupation numbers, respectively, obtained by diagonalising the Kohn–Sham Hamiltonian including SOC in a basis of collinear spinors aligned along the direction $\hat{\mathbf n}$, while $\varepsilon_{{\bf k}n}^0$ and $f(\varepsilon_{{\bf k}n}^0)$ are the bare Kohn–Sham eigenenergies and occupation numbers without SOC.

For all magnetic materials we have calculated the energy difference between out-of-plane and in-plane magnetisation $E_{{\rm MA}}(i)=\Delta E(\hat{\mathbf z})-\Delta E(i)$, ($i=\hat{\mathbf x}, \hat{\mathbf y}$). Negative values of $E_{{\rm MA}}(i)$ thus indicate that there is an out-of-plane easy axis of magnetisation.

Calculations for the ground state have been performed with plane-wave cutoff and energetic convergence threshold set to 800 eV and 0.5 meV/atom respectively. For all calculations we have used a Γ-centered Monkhorst–Pack k-point with a density of 20/$\mathring{\rm A}^{-1}$. The SOC contribution is introduced via a non-self-consistent diagonalisation of the Kohn–Sham Hamiltonian evaluated in the projector-augmented wave formalism [97].

The projected density of states (PDOS) is a useful tool for identifying which atomic orbitals comprise a band. It is defined as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:pdos_paw} \rho^S_l(\varepsilon) = \sum_{a \in S} \sum_{{\bf k}n}\sum_{m} |\langle \phi_{l, m}^a|\psi_{{\bf k}n}\rangle|^2 \delta(\varepsilon - \varepsilon_{{\bf k}n}), \nonumber \end{align} \tag{ 5 }$$

where $\psi_{{\bf k}n}$ are the Kohn–Sham wave functions with eigenvalues $\varepsilon_{{\bf k }n }$ and $\phi_{l, m}^a$ are the spin-paired Kohn–Sham orbitals of atomic species S with angular momentum l ($s, p, d, f$). We sum over all atoms belonging to species S so every atomic species has one entry per angular momentum channel. In the PAW formalism this can be approximated as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:pdos_paw2} \rho^S_l(\varepsilon) = \sum_{a \in S} \sum_{{\bf k}n}\sum_{m} |\langle \tilde{p}_{l, m}^a|\tilde{\psi}_{{\bf k}n}\rangle|^2 \delta(\varepsilon - \varepsilon_{{\bf k}n}) \nonumber \end{align} \tag{ 6 }$$

where $\tilde{\psi}_{{\bf k}n}$ are the pseudo wave functions and $\tilde{p}_{l, m}^a$ are the PAW projectors associated with the atomic orbitals $\phi_{l, m}^a$. The PDOS is calculated from equation (6) using linear tetrahedron interpolation [98] (LTI) of energy eigenvalues obtained from a ground state calculation with a k-point sampling of $36/\mathring{\rm A}^{-1}$. In contrast to other techniques for calculating the PDOS using smearing, the PDOS yielded by the LTI method returns exactly zero at energies with no states. Examples of PDOS are shown in figure 9 (right) for respectively the ferromagnetic metal VO₂ and the semiconductor WS₂ in the H phase (MoS₂ prototype).

Electronic band structures are calculated along the high symmetry paths shown in figure 8 for the five different types of 2D Bravais lattices. The band energies are computed within DFT using three different xc-functionals, namely PBE, HSE06, and GLLBSC. These single-particle approaches are complemented by many-body ${\rm G_0W_0}$ calculations for materials with a finite gap and up to four atoms in the unit cell (currently around 250 materials). For all methods, SOC is included by non-selfconsistent diagonalisation in the full basis of Kohn–Sham eigenstates. Band energies always refer to the vacuum level defined as the asymptotic limit of the Hartree potential, see figure 12. Below we outline the employed methodology while section 3.2.1 provides an overview and comparison of the band energies obtained with the different methods.

Zoom In Zoom Out Reset image size

2.8.1. PBE band structure

The electron density is determined self-consistently on a uniform k-point grid of density 12.0/$\mathring{\rm A}^{-1}$. From this density, the PBE band structure is computed non-selfconsistently at 400 k-points distributed along the band path (see figure 8). Examples of PBE band structures are shown in figure 9 for the ferromagnetic metal VO₂ and the semiconductor WS₂ both in the MoS₂ prototype structure. The expectation value of the out-of-plane spin component, $\langle \chi_{n{\bf k}\sigma}|\hat S_z|\chi_{n {\bf k} \sigma}\rangle$, is evaluated for each spinorial wave function, $\chi_{n{\bf k}\sigma}=(\psi_{n {\bf k} \uparrow}, \psi_{n {\bf k} \downarrow})$, and is indicated by the color of the band. For materials with inversion symmetry, the SOC cannot induce band splitting, meaning that $\langle \chi_{n {\bf k} \sigma}|\hat S_z|\chi_{n {\bf k} \sigma}\rangle$ is ill-defined and no color coding is used. The band structure without SOC is indicated by a dashed grey line. We have compared our PBE  +  SOC band gaps of 29 different monolayers with those obtained with the VASP code in [41] and find a mean absolute deviation of 0.041 eV.

Zoom In Zoom Out Reset image size

2.8.2. HSE band structure

The band structure is calculated non-selfconsistently using the range-separated hybrid functional HSE06 [99] on top of a PBE calculation with k-point density 12.0/$\mathring{\rm A}^{-1}$ and 800 eV plane wave cutoff. We have checked for selected systems that the HSE band structure is well converged with these settings. The energies along the band path are obtained by spline interpolation from the uniform k-point grid. As an example, the HSE band structure of WS₂ is shown in the left panel of figure 10 (black line) together with the PBE band structure (grey dashed). The PBE band gap increases from 1.52 eV to 2.05 eV with the HSE06 functional in good agreement with earlier work reporting band gaps of 1.50 eV (PBE) and 1.90 eV (HSE) [100] and 1.55 eV (PBE) and 1.98 eV (HSE) [101], respectively. A more systematic comparison of our results with the HSE  +  SOC band gaps obtained with the VASP code in [41] for 29 monolayers yield a mean absolute deviation of 0.14 eV. We suspect this small but non-zero deviation is due to differences in the employed PAW potentials and the non-selfconsistent treatment of the HSE in our calculations.

Zoom In Zoom Out Reset image size

2.8.3. GLLBSC fundamental gap

2.8.4. G₀W₀ band structure

For materials with finite PBE band gap the quasiparticle (QP) band structure is calculated using the ${\rm G_0W_0}$ approximation on top of PBE following our earlier work [105, 106]. Currently, this resource demanding step is performed only for materials with up to four atoms in the unit cell. The number of plane waves and the number of unoccupied bands included in the calculation of the non-interacting density response function and the GW self-energy are always set equal. The individual QP energies are extrapolated to the infinite basis set limit from calculations at plane wave cutoffs of 170, 185 and 200 eV, following the standard $1/N_\mathrm{G}$ dependence [107, 108], see figure 11 (right). The screened Coulomb interaction is represented on a non-linear real frequency grid ranging from 0 eV to 230 eV and includes around 250 frequency points. The exchange contribution to the self-energy is calculated using a Wigner–Seitz truncation scheme [109] for a more efficient treatment of the long range part of the exchange potential. For the correlation part of the self-energy, a 2D truncation of the Coulomb interaction is used [110, 111]. We stress that the use of a truncated Coulomb interaction is essential to avoid unphysical screening from periodically repeated layers which otherwise leads to significant band gap reductions.

Zoom In Zoom Out Reset image size

Importantly, the use of a truncated Coulomb interaction leads to much slower k-point convergence because of the strong q-dependence of the 2D dielectric function around q  =  0. We alleviate this problem by using an analytical expression for the screened interaction when performing the BZ integral around q  =  0 [106]. This allows us to obtain well converged results with a relatively low k-point density of 5.0/$\mathring{\rm A}^{-1}$ (corresponding to $12\times12$ k-points for MoS₂). For example, with this setting the ${\rm G_0W_0}$ band gap of MoS₂ is converged to within 0.05 eV, see figure 11 (left). In comparison, standard BZ sampling with no special treatment of the q  =  0 limit, requires around $40\times40$ k-points to reach the same accuracy.

Figure 10 (right) shows the PBE and ${\rm G_0W_0}$ band structures of WS₂ (including SOC). The ${\rm G_0W_0}$ self-energy opens the PBE band gap by 1.00 eV and the HSE gap by 0.47 eV, in good agreement with previous studies [112]. We note in passing that our previously published ${\rm G_0W_0}$ band gaps for 51 monolayer TMDCs [36] are in good agreement with the results obtained using the workflow described here. The mean absolute error between the two data sets is around 0.1 eV and can be ascribed to the use of PBE rather than LDA as starting point and the use of the analytical expression for W around q  =  0.

A detailed comparison of our results with previously published ${\rm G_0W_0}$ data is not meaningful because of the rather large differences in the employed implementations/parameter settings. In particular, most reported calculations do not employ a truncated Coulomb interaction and thus suffer from spurious screening effects, which are then corrected for in different ways. Moreover, they differ in the amount of vacuum included the supercell, the employed k-point grids and basis sets, the in-plane lattice constants, and the DFT starting points. For example, published values for the QP band gap of monolayer MoS₂ vary from from 2.40 to 2.90 eV [113–120] (see [119] for a detailed overview). The rather large variation in published GW results for 2D materials is a result of the significant numerical complexity of such calculations and underlines the importance of establishing large and consistently produced benchmark data sets like the present.

For bulk materials, self-consistency in the Green's function part of the self-energy, i.e. the GW₀ method, has been shown to increase the ${\rm G_0W_0}$ band gaps and improve the agreement with experiments [121]. The trend of band gap opening is also observed for 2D materials [106, 120, 120, 122], however, no systematic improvement with respect to experiments has been established [122]. For both bulk and 2D materials, the fully self-consistent GW self-energy systematically overestimates the band gap [121, 122] due to the neglect of vertex corrections [122, 123]. In ${\rm G_0W_0}$ the neglect of vertex corrections is partially compensated by the smaller band gap of the non-interacting Kohn–Sham Green's function compared to the true interacting Green's function. In this case, the vertex corrections corrections will affect mainly the absolute position of the bands relative to vacuum while the effect on the band gap is relatively minor [122].

In table 4 we compare calculated band gaps from C2DB with experimental band gaps for three monolayer TMDCs and phosphorene. The experimental data has been corrected for substrate interactions [122, 124], but not for zero-point motion, which is expected to be small (<$0.1 \mathrm{~eV}$). The ${\rm G_0W_0}$ results are all within 0.2 eV of the experiments. A further (indirect) test of the ${\rm G_0W_0}$ band gaps against experimental values is provided by the comparison of our BSE spectra against experimental photoluminescence data in table 7, where we have used a ${\rm G_0W_0}$ scissors operator. Finally, we stress that the employed PAW potentials are not norm-conserving, and this can lead to errors for bands with highly localised states (mainly 4f and 3d orbitals), as shown in [108]. Inclusion of vertex corrections and use of norm conserving potentials will be the focus of future work on the C2DB.

Table 4. Comparison between calculated and experimental band gaps for four freestanding monolayers. The experimental values have been corrected for substrate screening. MAD refers to the mean absolute deviation between the predicted values and the measured values.

Material	Band gap (eV)
	PBE	HSE06	GLLBSC	${\rm G_0W_0}$	Experiment
MoS2	1.58	2.09	2.21	2.53	2.50 [125]
MoSe2	1.32	1.80	1.88	2.12	2.31 [126]
WS2	1.53	2.05	2.16	2.53	2.72 [127]
P (phosphorene)	0.90	1.51	1.75	2.03	2.20 [124]
MAD w.r.t. experiment	1.10	0.57	0.43	0.15	—

For materials with a finite band gap, the positions of the valence band maximum (VBM) and conduction band minimum (CBM) within the BZ are identified together with their energies relative to the vacuum level. The latter is defined as the asymptotic value of the electrostatic potential, see figure 12. The PBE electrostatic potential is used to define the vacuum level in the non-selfconsistent HSE and G₀W₀ calculations. For materials with an out-of-plane dipole moment, a dipole correction is applied during the selfconsistent DFT calculation, and the vacuum level is defined as the average of the asymptotic electrostatic potentials on the two sides of the structure. The PBE vacuum level shift is also stored in the database.

Zoom In Zoom Out Reset image size

The Fermi surface is calculated using the PBE xc-functional including SOC for all metallic compounds in the database. Based on a ground state calculation with a k-point density of at least 20/$\mathring{\rm A}^{-1}$, the eigenvalues are interpolated with quadratic splines and plotted within the first BZ. Figure 13 (left) shows an example of the Fermi surface for VO₂-MoS₂ with color code indicating the out-of-plane spin projection $\langle S_z \rangle$. The band structure refers to the ferromagnetic ground state of VO₂-MoS₂, which has a magnetic moment of 0.70 $\mu_{\rm B}$ per unit cell, characterised by alternating spin-polarised lobes with $\langle S_z \rangle = \pm 1$.

Zoom In Zoom Out Reset image size

For materials with a finite PBE gap, the effective electron and hole masses are calculated from the PBE eigenvalues; initially these are calculated on an ultrafine k-point mesh of density $45.0/\mathring{\rm A}^{-1}$ uniformly distributed inside a circle of radius 0.015 $\mathring{\rm A}^{-1}$ centered at the VBM and CBM, respectively. The radius is chosen to be safely above the noise level of the calculated eigenvalues but still within the harmonic regime; it corresponds to a spread of eigenvalues of about 1 meV within the circle for an effective mass of 1 m₀. For each band within an energy window of 100 meV above/below the CBM/VBM, the band curvature is obtained by fitting a third order polynomial. Even though the masses represent the second derivative of the band energies, we have found that the inclusion of 3rd order terms stabilises the fitting procedure and yields masses that are less sensitive to the details of the employed k-point grids. For each band the mass tensor is diagonalised to yield heavy and light masses in case of anisotropic band curvatures. The masses (in two directions) and the energetic splitting of all bands within 100 meV of the band extremum are calculated both with and without SOC and stored in the database. Other approaches exist for calculating effective masses, such as $\mathbf{k} \cdot \mathbf{p}$ perturbation theory (see e.g. [128] and references therein); the present scheme was chosen for its simplicity and ease of application to a wide range of different materials.

In addition to the effective masses at the VBM and CBM, the exciton reduced mass is calculated by applying the above procedure to the direct valence-conduction band transition energies, $\varepsilon_{\mathrm{v-c}}(\bf k)=\varepsilon_{\mathrm{c}}(k)-\varepsilon_{\mathrm{v}}(k)$. For direct band gap materials the exciton reduced mass is related to the electron and hole masses by $1/\mu_{\mathrm{ex}}=1/m^*_{\mathrm{e}}+1/m^*_{\mathrm{h}}$, but in the more typical case of indirect band gaps, this relation does not hold.

As an example, figure 14 shows a zoom of the band structure of SnS-GeSe around the VBM and CBM (upper and lower panels). The second order fits to the band energies (extracted from the fitted 3rd order polynomial) are shown by red dashed lines. It can be seen that both the conduction and valence bands are anisotropic leading to a heavy and light mass direction (left and right panels, respectively). The valence band is split by the SOC resulting in two bands separated by $\sim 10$ meV and with slightly different curvatures. The conduction band exhibits a non-trivial band splitting in one of the two directions. The peculiar band splitting is due to a Rashba effect arising from the combination of spin–orbit coupling and the finite perpendicular electric field created by the permanent dipole of the SnS structure where Sn and S atoms are displaced in the out of plane direction leading to a sizable vacuum level difference of 1.13 eV, see figure 12.

Zoom In Zoom Out Reset image size

Table 5 shows a comparison between selected effective masses from the C2DB and previously published data also obtained with the PBE xc-correlation functional and including SOC. Overall, the agreement is very satisfactory.

For metallic compounds, the work function is obtained as the difference between the Fermi energy and the asymptotic value of the electrostatic potential in the vacuum region, see figure 12. The work function is determined for both PBE and HSE band structures (both including SOC) on a uniform k-point grid of density 12.0/$\mathring{\rm A}^{-1}$. Since the SOC is evaluated non-selfconsistently, the Fermi energy is adjusted afterwards based on a charge neutrality condition.

For semiconductors, the deformation potentials quantify the shift in band edge energies (VBM or CBM) upon a linear deformation of the lattice. The uniaxial absolute deformation potential along axis i ($i=x, y$) is defined as [129, 130]

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:deformation} D_{ii}^\alpha=\frac{\Delta E_\alpha}{\varepsilon_{ii}},\quad \alpha=\mathrm{VBM,CBM} \nonumber \end{align} \tag{ 7 }$$

where $\Delta E_\alpha$ is the energy shift upon strain and $\varepsilon_{ii}$ are the strains in the i-directions.

The deformation potentials are important physical quantities as they provide an estimate of the strength of the (acoustic) electron-phonon interaction, see section 3.2.2. Moreover, they are obviously of interest in the context of strain-engineering of band gaps and they can be used to can be used to infer an error bar on the band gap or band edge positions due to a known or estimated error bar on the lattice constant.

The calculation of $D^\alpha_{ii}$ is based on a central difference approximation to the derivative. A strain of $\pm$$1\%$ is applied separately in the x and y directions and the ions are allowed to relax while keeping the unit cell fixed. Calculations are performed with the PBE xc-functional, a plane wave cutoff of 800 eV, and a k-point density of 12/$\mathring{\rm A}^{-1}$.

The change in band energy, $\Delta E_\alpha$, is measured relative to the vacuum level. In cases with nearly degenerate bands, care must be taken to track the correct bands as different bands might cross under strain. In this case, we use the expectation value $\langle \hat S_z \rangle$ to follow the correct band under strain. Figure 15 shows how the band structure of MoS₂ changes as a function of strain. Both the VBM and the CBM shift down (relative to the vacuum level) when tensile strain is applied in the x direction, but the conduction band shows a much larger shift, leading to an effective band gap closing under tensile strain.

Zoom In Zoom Out Reset image size

Table 6 shows a comparison between the deformation potentials in the C2DB, and literature values obtained using similar methods. There is generally good agreement, and part of the discrepancy can be ascribed to the fact that, in contrast to [131], our numbers include spin–orbit coupling.

The dielectric response of a 2D material is described by its 2D polarisability, $\alpha^\mathrm{2D}$ (see section 2.15 for a general introduction of this quantity). For metals, it can be separated into contributions from intraband and interband transitions, i.e. $\alpha^\mathrm{2D} = \alpha^\mathrm{2D, intra} + \alpha^\mathrm{2D, inter}$. We have found that local field effects (LFEs) are negligible for the intraband component, which consequently can be treated separately and evaluated as an integral over the Fermi surface. Specifically, this leads to the Drude expression for the polarisability in the long wave length limit $\alpha^\mathrm{2D, intra}(\omega) = -\omega_\mathrm{P, 2D}^2 / (2\pi\omega^2)$ where $\omega_\mathrm{P, 2D}$ is the 2D plasma frequency, which in atomic units is given by

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \omega^2_\mathrm{P,2D} = \frac{4\pi}{A} \sum_{sn\mathbf{k}} | \hat{\mathbf{q}} \cdot \mathbf{v}_{sn\mathbf{k}}|^2 \delta(\varepsilon_{sn\mathbf{k}} - \varepsilon_\mathrm{F}), \nonumber \end{align} \tag{ 8 }$$

where $\mathbf{v}_{sn\mathbf{k}}= \langle sn\mathbf{k} | {\hat{\mathbf p}} / m_0 | sn\mathbf{k} \rangle$ is a velocity matrix element (with m₀ the electron mass), $\hat{\mathbf{q}} = \mathbf{q} / q$ is the polarisation direction, $s, n, \mathbf{k}$ denote spin, band and momentum indices, and A is the supercell area. The 2D plasma frequency is related to the conventional 3D plasma frequency by $\omega_\mathrm{P, 2D}^2(\omega) = \omega_\mathrm{P, 3D}^2(\omega) L / 2$ where L is the supercell height.

The plasma frequency defined above determines the intraband response of the 2D metal to external fields. In particular, it determines the dispersion relation of plasmon excitations in the metal. The latter are defined by the condition $\varepsilon^\mathrm{2D}(\omega_{\rm P}) = 1 + 2 \pi q \alpha^{2{\rm D}}(\omega_{\rm P}) = 0$, where q is the plasmon wave vector. Neglecting interband transitions (the effect of which is considered in section 3.2.4), the 2D plasmon dispersion relation becomes

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \omega_\mathrm{P}(q) = \omega_{\mathrm{P,2D}} \sqrt{q} \label{eq:plasmonfrequency}. \nonumber \end{align} \tag{ 9 }$$

The plasma frequencies, $\omega_{\mathrm{P, 2D}}$, for polarisation in the x and y directions, respectively, are calculated for all metals in the C2DB using the linear tetrahedron method [98] to interpolate matrix elements and eigenvalues based on a PBE band structure calculation with a k-point density of $20/\mathring{\rm A}^{-1}$.

The polarisability tensor $\alpha_{ij}$ is defined by

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle P_i(\mathbf{q},\omega)=\sum_j\alpha_{ij}(\mathbf{q},\omega)E_j(\mathbf{q},\omega), \nonumber \end{align} \tag{ 10 }$$

where P_i is the i'th component of the induced polarisation averaged over a unit cell and E_j is the j'th component of the macroscopic electric field. Using that $\newcommand{\e}{{\rm e}} P_i=(D_i-E_i)/(4\pi)=\sum_j(\epsilon_{ij}-\delta_{ij})E_j/(4\pi)$ one observes that $\newcommand{\e}{{\rm e}} \alpha_{ij}=(\epsilon_{ij}-\delta_{ij})/(4\pi)$, where $\newcommand{\e}{{\rm e}} \epsilon_{ij}$ is the dielectric function. In contrast to the dielectric function, whose definition for a 2D material is not straightforward [119], the polarisability allows for a natural generalisation to 2D by considering the induced dipole moment per unit area,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle P_i^{\mathrm{2D}}(\mathbf{q},\omega)=\sum_j\alpha^{\mathrm{2D}}_{ij}(\mathbf{q},\omega)E_j(\mathbf{q},\omega). \nonumber \end{align} \tag{ 11 }$$

Since the P_i is a full unit cell average and $P_i^{\mathrm{2D}}$ is integrated in the direction orthogonal to the slab, we have $P_i^{\mathrm{2D}}=LP_i$ and $\alpha^{\mathrm{2D}}_{ij}=L\alpha_{ij}$, where L is the length of the unit cell in the direction orthogonal to the slab.

In the following, we focus on the longitudinal components of the polarisability and dielectric tensors, which are simply denoted by α and . These are related to the density-density response function, χ, via the relations

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \alpha^{\mathrm{2D}}(\mathbf{q},\omega)=\frac{L}{4\pi}(\epsilon(\mathbf{q},\omega)-1),\label{eq:2Dpol} \nonumber \end{align} \tag{ 12 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \epsilon^{-1}(\mathbf{q},\omega)=1+\langle v_c(\mathbf{q})\chi(\omega)\rangle_\mathbf{q}, \nonumber \end{align} \tag{ 13 }$$

where $v_c$ is the Coulomb interaction and

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:average} \langle v_c\chi(\omega)\rangle_\mathbf{q}=\frac{1}{V}\int_{\rm Cell} d\mathbf{r}d\mathbf{r'}d\mathbf{r''}v_c(\mathbf{r},\mathbf{r'})\chi(\mathbf{r'},\mathbf{r''},\omega)e^{-i\mathbf{q}(\mathbf{r}-\mathbf{r''})}, \nonumber \end{align} \tag{ 14 }$$

where Cell is the supercell with volume $V$. The re- sponse function, χ, satisfies the Dyson equation [132] $\chi=\chi^{\mathrm{irr}} + \chi^{\mathrm{irr}}v_c\chi,$ where $\chi^{\mathrm{irr}}$ is the irreducible density-density response function. In the random phase approximation (RPA) $\chi^{\mathrm{irr}}$ is replaced by the non-interacting response function, $\chi^0$, whose plane wave representation is given by [133, 134]

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:chi0} &\chi^0_{\mathbf{G}\mathbf{G}'}(\mathbf{q}, \omega)=\frac{1}{\Omega}\sum_{\mathbf{k}\in{{\rm BZ}}}\sum_{mn}(\,f_{n\mathbf{k}}-f_{m\mathbf{k+q}})\nonumber \\ &\times\frac{\langle\psi_{n\mathbf{k}}|e^{-i (\mathbf{q} + \mathbf{G})\cdot{\bf r}}|\psi_{m\mathbf{k+q}}\rangle\langle\psi_{m\mathbf{k+q}}|e^{i (\mathbf{q} + \mathbf{G}')\cdot{\bf r}}|\psi_{n\mathbf{k}}\rangle}{\hbar\omega+\varepsilon_{n\mathbf{k}}-\varepsilon_{m\mathbf{k+q}}+i\eta}, \nonumber \end{align} \tag{ 15 }$$

where $\mathbf{G}, \mathbf{G}'$ are reciprocal lattice vectors and Ω is the crystal volume.

For all materials in the database, we calculate the polarisability within the RPA for both in-plane and out-of-plane polarisation in the optical limit ${\rm q}\rightarrow {0}$. For metals, the interband contribution to the polarisability is obtained from equation (15) while the intraband contribution is treated separately as described in section 2.14. The single-particle eigenvalues and eigenstates used in equation (15) are calculated with PBE, a k-point density of $20/\mathring{\rm A}^{-1}$ (corresponding to a k-point grid of $48\times48$ for MoS₂ and $60\times60$ for graphene), and 800 eV plane wave cutoff. The Dyson equation is solved using a truncated Coulomb potential [105, 111] to avoid spurious interactions from neighboring images. We use the tetrahedron method to interpolate the eigenvalues and eigenstates and a peak broadening of $\newcommand{\e}{{\rm e}} \eta=50$ meV. Local field effects are accounted for by including $\mathbf G$-vectors up to 50 eV. For the band summation we include 5 times as many unoccupied bands as occupied bands, which roughly corresponds to an energy cutoff of 50 eV. The calculations are performed without spin–orbit coupling.

In figure 16 we show the real and imaginary part of $\alpha^{\mathrm{2D}}$ for the semiconductor MoS₂. The PBE band gap of this material is 1.6 eV and we see the onset of dissipation at that energy. We also see that the initial structure of $\mathrm{Im}~{\alpha}$ is a constant, which is exactly what would be expected from the density of states in a 2D material with parabolic dispersion. Finally, we note that the static polarisability $\left.\mathrm{Re}~\alpha\right|_{\omega=0}\approx6\mathring{\rm A}$, which can easily be read off the figure. The polarisability is also shown for the metallic 1T-NbS₂ where we display the real part with and without the intraband Drude contribution $\newcommand{\e}{{\rm e}} \omega^2_{\mathrm{P,2D}} /(\hbar\omega+i\eta){}^2$.

Zoom In Zoom Out Reset image size

The power absorbed by a 2D material under illumination of a monochromatic light field with polarisation $\hat{\mathbf{e}}$ is quantified by the dimensionless absorbance:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:abs} {\rm Abs}(\omega)=4\pi\omega\alpha^{\mathrm{2D}}(q\hat{\mathbf{e}}\rightarrow0,\omega) / c, \nonumber \end{align} \tag{ 16 }$$

where c is the speed of light. In section 2.15 we gave a prescription for evaluating $\alpha^{\mathrm{2D}}$ in the RPA. However, absorption spectra of 2D semiconductors often display pronounced excitonic effects, which are not captured by the RPA. The Bethe–Salpeter equation (BSE) is a well-known method capable of describing excitonic effects and has been shown to provide good agreement with experimental absorption spectra for a wide range of materials [135].

For materials with finite band gap and up to four atoms per unit cell, we have calculated the RPA and the BSE absorption spectra for electric fields polarised parallel and perpendicular to the layers. The calculations are performed on top of PBE eigenstates and eigenvalues with spin–orbit coupling included and all unoccupied band energies shifted by a constant in order to reproduce the ${\rm G_0W_0}$ quasiparticle gap (the scissors operator method). If the ${\rm G_0W_0}$ gap is not available we use the GLLBSC gap for non-magnetic materials and the HSE gap for magnetic materials (since GLLBSC is not implemented in GPAW for spin-polarised systems). The screened interaction entering the BSE Hamiltonian is calculated within the RPA using a non-interacting response function evaluated from equation (15) with local field effects (i.e. $\mathbf G$-vectors) included up to 50 eV and 5 times as many unoccupied bands as occupied bands for the sum over states. We apply a peak broadening of $\newcommand{\e}{{\rm e}} \eta=50$ meV and use a truncated Coulomb interaction. The BSE Hamiltonian is constructed from the four highest occupied and four lowest unoccupied bands on a k-point grid of density of $20/\mathring{\rm A}^{-1}$, and is diagonalised within the Tamm–Dancoff approximation. We note that the Tamm–Dancoff approximation has been found to be very accurate for bulk semiconductors [136]. For monolayer MoS₂ we have checked that it reproduces the full solution of the BSE, but its general validity for 2D materials, in particular those with small band gaps, should be more thoroughly tested.

In figure 17 we show the optical absorption spectrum of MoS₂ obtained with the electric field polarised parallel and perpendicular to the layer, respectively. Both RPA and BSE spectra are shown (the in-plane RPA absorbance equals the imaginary part of the RPA polarisability, see figure 16 (left), apart from the factor $4 \pi \omega$ and the scissors operator shift). The low energy part of the in-plane BSE spectrum is dominated by a double exciton peak (the so-called A and B excitons) and is in excellent agreement with experiments [55].

Zoom In Zoom Out Reset image size

In general, calculations of electronic excitations of 2D materials converge rather slowly with respect to k-points due to the non-analytic behavior of the dielectric function in the vicinity of q  =  0 [119, 137, 138]. In figure 18 we show the k-point dependence of the binding energy of the A exciton in MoS₂ obtained as the difference between the direct band gap and the position of the first peak in figure 17. We observe a strong overestimation of the exciton binding energy at small k-point samplings, which converges slowly to a value of  ∼0.5 eV at large k-point samplings. For $48\times 48$ k-points, corresponding to the k-point sampling used for the BSE calculations in the database, the exciton binding energy is 0.53 eV, whereas a $1 / N_k^2$ extrapolation to infinite k-point sampling gives 0.47 eV (see inset in figure 18). In general, the exciton binding energy decreases with increasing k-point sampling, and thus the exciton binding energies reported in the C2DB might be slightly overestimated. However, since ${\rm G_0W_0}$ band gaps also decrease when the k-point sampling is increased (see figure 11) the two errors tend to cancel such that the absolute position of the absorption peak from BSE-${\rm G_0W_0}$ converges faster than the band gap or exciton binding energy alone.

Zoom In Zoom Out Reset image size

The BSE-${\rm G_0W_0}$ method has previously been shown to provide good agreement with experimental photoluminescence and absorption measurements on 2D semiconductors. In table 7 we show that our calculated position of the first excitonic peak agree well with experimental observations for four different TMDCs and phosphorene. Experimentally, the monolayers are typically supported by a substrate, which may alter the screening of excitons. However the resulting decrease in exciton binding energies is largely cancelled by a reduced quasiparticle gap such that the positions of the excitons are only slightly red-shifted as compared with the case of pristine monolayers [139].

Table 7. Comparison between calculated and experimental positions of the first excitonic peak for four different transition metal dichalcogenide monolayers and phosphorene.

Material	Energy of the first bright exciton (eV)
	BSE@PBE-${\rm G_0W_0}$	Experiment
MoS2	2.00	1.83 [140], 1.86 [141], 1.87 [142]
MoSe2	1.62	1.57 [140], 1.57 [143], 1.58 [144]
WS2	2.07	1.96 [141], 2.02 [144]
WSe2	1.71	1.64 [142], 1.66 [143]
P (phosphorene)	1.45	1.45 [145], 1.75 [146]
MAD w.r.t. experiment	0.066	—

In table 8 we list the major classes of materials currently included in the database. The materials are grouped according to their prototype, see section 2.2.2. For each prototype we list the corresponding space group, the total number of materials, and the number of materials satisfying a range of different conditions. The atomic structure of some of the different prototypes were shown in figure 4. The vast majority of the 2D materials that have been experimentally synthesised in monolayer form are contained in the C2DB (the 55 materials in figure 7 in addition to seven metal-organic perovskites). These materials are marked in the database and a literature reference is provided. Additionally, 80 of the monolayers in the C2DB could potentially be exfoliated from experimentally known layered bulk structures [16–19]. These materials are also marked and the ID of the bulk compound in the relevant experimental crystal structure database is provided.

Table 8. Overview of the materials currently in the C2DB. The table shows the number of compounds listed by their crystal structure prototype and selected properties. $E_{\mathrm{gap}} > 0$ and 'direct gap' refer to the PBE values, 'high stability' refers to the stability scale defined in section 2.4.1, and the last three columns refer to the magnetic state, see section 2.1. In this overview, separate magnetic phases of the same structure are considered different materials.

Prototype	Symmetry	Number of materials
		Total	$E_{\mathrm{gap}} > 0$	Direct gap	High stability	NM	FM	AFM
C	P6/mmm	4	4	3	1	4	0	0
CH	P$\overline{3}$m1	8	7	6	1	8	0	0
CH2Si	P3m1	2	2	2	1	2	0	0
BN	P$\overline{3}$m2	10	9	5	1	10	0	0
GaS	P$\overline{3}$m2	125	34	95	8	100	18	7
FeSe	P4/nmm	103	13	90	26	74	18	11
GeSe	P3m1	20	19	5	6	20	0	0
PbSe	P4/mmm	44	6	38	1	33	8	3
P	Pmna	9	9	0	1	9	0	0
MoS2	P$\overline{3}$m2	241	85	176	53	156	85	0
CdI2	P$\overline{3}$m1	315	95	231	90	218	80	17
WTe2	P21/m	75	29	48	34	57	13	5
FeOCl	Pmmn	443	92	385	65	328	63	52
MoSSe	P3m1	9	6	6	5	8	1	0
C3N	P6/mmm	25	1	24	0	25	0	0
YBr3	P6/mmm	57	11	51	0	21	24	12
TiCl3	P$\overline{3}$2m	69	35	51	2	32	23	14
BiI3	P$\overline{3}$m1	123	69	66	15	48	54	21
TiS3	Pmmn	34	8	28	5	31	2	1
MnTe3	P21/m	29	3	27	1	22	4	3
Cr3WS8	Pmm2	35	34	18	8	35	0	0
CrWS4	Pmm2	18	17	7	8	18	0	0
Ti2CO2	P$\overline{3}$m1	28	8	20	12	19	6	3
Ti2CH2O2	P$\overline{3}$m1	13	3	12	3	10	2	1
Ti3C2	P$\overline{3}$m2	12	0	12	0	7	5	0
Ti3C2O2	P$\overline{3}$m2	26	0	26	0	20	6	0
Ti3C2H2O2	P$\overline{3}$m2	14	0	14	0	10	4	0
PbA2I4	P1	27	27	27	0	27	0	0
Sum		1918	626	151	347	1352	416	150

To illustrate how all the materials are distributed in terms of stability, we show the energy above the convex hull plotted against $\tilde{\omega}^2_{\mathrm{min}}$ in figure 19. It can be seen that the structures naturally sort themselves into two clusters according to the dynamic stability. The points have been colored according to the three levels for dynamic stability introduced in section 2.4. The lower panel shows the distribution of the materials in the grey region on a linear scale. While most of the experimentally known materials (red and black dots) have high dynamic stability, a significant part of them fall into the medium stability category. The marginal distributions on the plot show that the more dynamically stable materials are also more thermodynamically stable. The mean energy above the convex hull is 0.12 eV for the materials with high dynamical stability, while it is 0.25 eV for the others.

Zoom In Zoom Out Reset image size

In table 9 we show the key properties of a selected set of stable materials, distributed across a variety of different crystal structure prototypes. To our knowledge, these materials are not experimentally known, and they are therefore promising candidates further study and experimental synthesis.

Table 9. Key properties of selected stable materials in the C2DB, which have not been previously synthesised. The calculated properties are the magnetic state, formation energy, energy above the convex hull, work function, PBE gap and and the nature of the gap (direct or indirect).

Prototype	Formula	Magnetic state	$\Delta H$ (eV)	$\Delta H_{\mathrm{hull}}$ (eV)	Φ (eV)	PBE gap (eV)	Direct gap
BiI3	VI3	FM	−0.51	−0.15	5.3
	CoCl3	NM	−0.65	−0.21		1.13	No
	CoBr3	NM	−0.41	−0.16		0.96	No
	CoI3	NM	−0.14	−0.14		0.53	No
CdI2	FeO2	FM	−1.14	−0.36	7.31
	MnSe2	FM	−0.47	−0.18	5.09
	MnS2	FM	−0.57	−0.12	5.74
	PdO2	NM	−0.40	−0.08		1.38	No
	CaBr2	NM	−2.09	−0.02		4.86	No
FeOCl	RhClO	NM	−0.65	−0.18	5.49
	NiClO	AFM	−0.64	−0.17	6.32
	NiBrO	AFM	−0.52	−0.16	5.78
	ScIS	NM	−1.68	−0.14		1.66	Yes
FeSe	CoSe	FM	−0.27	0.02	4.22
	RuS	NM	−0.38	0.05	4.72
	MnSe	AFM	−0.50	−0.20		0.90	No
	MnS	AFM	−0.64	−0.19		0.78	No
GaS	AlSe	NM	−0.72	−0.02		2.00	No
	AlS	NM	−0.89	0.00		2.09	No
GeSe	GeSe	NM	−0.19	0.04		2.22	No
	GeS	NM	−0.22	0.05		2.45	No
	GeTe	NM	−0.01	0.09		1.47	No
	SnSe	NM	−0.33	0.10		2.15	No
MoS2	VS2	FM	−0.88	−0.02	5.95
	ScBr2	FM	−1.59	−0.40		0.16	No
	YBr2	FM	−1.73	−0.23		0.34	No
	FeCl2	FM	−0.67	−0.16		0.35	Yes
	TiBr2	NM	−1.14	−0.04		0.76	No
	ZrBr2	NM	−1.34	−0.04		0.83	No
Ti2CO2	Zr2CF2	NM	−2.36	−0.08	3.92
	Hf2CF2	NM	−2.26	0.03	3.62
	Y2CF2	NM	−2.50	−0.17		1.12	No
WTe2	NbI2	NM	−0.37	0.04	3.01
	HfBr2	NM	−1.16	−0.18		0.85	No
	OsSe2	NM	−0.17	0.00		0.57	No

In the following sections we present a series of case studies focusing on different properties of 2D materials including band gaps and band alignment, effective masses and mobility, magnetic order, plasmons and excitons. The purpose is not to provide an in-depth nor material specific analysis, but rather to explore trends and correlations in the data and showcase some potential applications of the C2DB. Along the way, we report some of the novel candidate materials revealed by this analysis, which could be interesting to explore closer in the future.

3.2.1. Band gaps and band alignment

The band gaps and band edge positions of all semiconductors and insulators in the C2DB have been calculated with the PBE, HSE06, and GLLBSC xc-functionals while ${\rm G_0W_0}$ calculations have been performed for the  ∼250 simplest materials. The relatively large size of these datasets and the high degree of consistency in the way they are generated (all calculations performed with the same code using same PAW potentials and basis set etc) provide a unique opportunity to benchmark the performance of the different xc-functionals against the more accurate ${\rm G_0W_0}$ method.

Figure 20 compares the size and center of the band gaps obtained with the density functionals to the ${\rm G_0W_0}$ results. Relative to ${\rm G_0W_0}$ the PBE functional underestimates the gaps by 45%, i.e. on average the PBE values must be scaled by 1.83 to reproduce the ${\rm G_0W_0}$ results. The HSE06 band gaps are closer to ${\rm G_0W_0}$ but are nevertheless systematically underestimated by more than 20%. On the other hand, GLLBSC shows very good performance with band gaps only 2% smaller than ${\rm G_0W_0}$ on average. Table 10 shows the mean absolute deviations of the DFT methods relative to ${\rm G_0W_0}$. We note that although GLLBSC provides an excellent description of the ${\rm G_0W_0}$ band gaps on average the spread is sizable with a mean absolute deviation of 0.4 eV.

Table 10. The mean absolute deviation (in eV) of the band gap and band gap center calculated with three different xc-functionals with respect to ${\rm G_0W_0}$.

	PBE	HSE06	GLLBSC
MAD w.r.t. ${\rm G_0W_0}$ (band gap)	1.49	0.82	0.38
MAD w.r.t. ${\rm G_0W_0}$ (gap center)	0.37	0.32	0.76

Zoom In Zoom Out Reset image size

We note a handful of outliers in figure 20 with large HSE band gaps compared to PBE and ${\rm G_0W_0}$. For one of these, namely the ferromagnetic CoBr₂-CdI₂, we obtain the band gaps: 0.30 eV (PBE), 3.41 eV (HSE), and 0.91 eV (G₀W₀). For validation, we have performed GPAW and QuantumEspresso calculations with the norm-conserving HGH pseudopotentials and plane wave cutoff up to 1600 eV. The converged band gaps are 0.49 eV (GPAW-HGH-PBE), 0.51 eV (QE-HGH-PBE) and 3.69 eV (GPAW-HGH-HSE), 3.52 eV (QE-HGH-HSE), which are all in reasonable agreement with the C2DB results. It should be interesting to explore the reason for the anomalous behavior of the HSE band gap in these materials.

Compared to the band gaps, the gap centers predicted by PBE and HSE06 are in overall better agreement with the ${\rm G_0W_0}$ results. This implies that, on average, the ${\rm G_0W_0}$ correction of the DFT band energies is symmetric on the valence and conduction band. In contrast, the GLLBSC predicts less accurate results for the gap center. This suggests that an accurate and efficient prediction of absolute band energies is obtained by combining the GLLBSC band gap with the PBE band gap center.

Next, we consider the band alignment at the interface between different 2D materials. Assuming that the bands line up with a common vacuum level and neglecting hybridisation/charge transfer at the interface, the band alignment is directly given by the VBM and CBM positions relative to vacuum.

We focus on pairs of 2D semiconductors for which the ${\rm G_0W_0}$ band alignment is either Type II ($\Delta E > 0$) or Type III ($\Delta E < 0$), see figure 21 (left). Out of approximately 10 000 bilayers predicted to have Type II band alignment by ${\rm G_0W_0}$, the PBE and HSE06 functionals predict qualitatively wrong band alignment (i.e. Type III) in 44% and 21% cases, respectively (grey shaded areas). In particular, PBE shows a sizable and systematic underestimation of $\Delta E$ as a direct consequence of the underestimation of the band gaps in both monolayers.

Zoom In Zoom Out Reset image size

3.2.2. Effective masses and mobilities

The carrier mobility relates the drift velocity of electrons or holes to the strength of an applied electric field and is among the most important parameters for a semiconductor material. In general, the mobility is a sample specific property which is highly dependent on the sample purity and geometry, and (for 2D materials) interactions with substrate or embedding layers. Here we consider the phonon-limited mobility, which can be considered as the intrinsic mobility of the material, i.e. the mobility that would be measured in the absence of any sample specific- or external scattering sources.

The effective masses of the charge carriers have been calculated both with and without SOC for  ∼600 semiconductors. Figure 22 shows the electron mass plotted against the hole mass. The data points are scattered, with no clear correlation between the electron and hole masses. Overall, the electron masses are generally slightly smaller than the hole masses. The mean electron mass is 0.9 m₀, while the mean hole mass is 1.1 m₀, and 80% of the electron masses are below m₀ while the fraction is only 65% for the holes. This is not too surprising, since, on average, the energetically lower valence band states are expected to be more localised and thus less dispersive than the conduction band states.

Zoom In Zoom Out Reset image size

The right panel of figure 22 shows the effective mass for electrons and holes plotted as a function of the inverse band gap. It can be seen that there is no clear correlation between the two quantities, which is confirmed by calculating the cross-correlation coefficient: for both electrons and holes it is less than 0.02. This provides empirical evidence against the linear relation between effective masses and inverse band gaps derived from $\mathbf{k} \cdot \mathbf{p}$ perturbation theory. The relation is based on the assumption that the perturbative expansion is dominated by the conduction and valence band and that the momentum matrix element between these states, $\newcommand{\bi}{\boldsymbol} \left\langle u_{\mathrm{c}} \big \vert \hat{\mathbf{p}}\big \vert u_{\mathrm{v}}\right\rangle$, does not vary too much as function of the considered parameter (here the type of material). These assumptions clearly do not hold across a large set of different semiconductors. If we focus on a specific class of materials, e.g. sulfides in the MoS₂ structure indicated by the highlighted symbols, we see a slightly improved trend but still with significant fluctuations.

If one assumes energetically isolated and parabolic bands, the intrinsic mobility limited only by scattering on acoustic phonons can be estimated from the Takagi formula [147],

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \mu_i = \frac{e \hbar^3 \rho v^2_i }{k_{\mathrm{B}} T m_i^{*}m^{*}_{\mathrm{d}} D_{i}^2}\label{eq:takagi}. \nonumber \end{align} \tag{ 17 }$$

Here i refers to the transport direction, ρ is the mass density, $v_i$ is the speed of sound in the material, $m_{i}^*$ is the carrier mass, $m_{\mathrm{d}}^*$ is the equivalent isotropic density-of-state mass defined as $m^*_{\mathrm{d}}=\sqrt{m_x^*m_y^*}$, and D_i is the deformation potential. We stress that the simple Takagi formula is only valid for temperatures high enough that the acoustic phonon population can be approximated by the Rayleigh–Jeans law, $n\approx \hbar \omega_{\mathrm{ac}} / k_{\mathrm{B}} T$, but low enough that scattering on optical phonons can be neglected.

For the semiconductors in the C2DB we have found that the denominator of equation (17) varies more than the numerator. Consequently, a small product of deformation potential and effective mass is expected to correlate with high mobility. Figure 23 shows the deformation potential plotted against the carrier mass for the valence and conduction bands, respectively. The shaded area corresponds, somewhat arbitrarily, to the region for which $m_i^* D_i < m_0(1 ~\mathrm{eV})$. The 2D semiconductors which have been synthesised in monolayer form are indicated with orange symbols while those which have been used in field effect transistors are labeled. Consistent with experimental findings, phosphorene (P) is predicted to be among the materials with the highest mobility for both electrons and holes.

Zoom In Zoom Out Reset image size

Interestingly, a number of previously unknown 2D materials lie in this shaded region and could be candidates for high mobility 2D semiconductors. Table 11 lists a few selected materials with high intrinsic mobility according to equation (17), which all have 'high' overall stability (see section 2.4.1). In the future, it will be interesting to explore the transport properties of these candidate materials in greater detail.

Table 11. Key transport properties of selected materials with high intrinsic room-temperature mobility according to equation (17). All the materials shown have 'high' overall stability as defined in section 2.4.1. $\mu_{\mathrm{high}}$ is the larger value of the mobilities in the x or y directions, m^* is the corresponding effective mass, and $\mu_{\mathrm{high}} / \mu_{\mathrm{low}}$ is the ratio of the mobilities in the two directions.

Carrier	Formula	Prototype	PBE gap (eV)	$\mu_{\mathrm{high}}$ (cm2 V−1 s−1)	m* (m0)	$\frac{\mu_{\mathrm{high}}}{\mu_{\mathrm{low}}}$
Holes	PbS2	MoS2	1.39	30 000	0.62	1.4
	OsO2	WTe2	0.17	23 000	0.23	3.0
	ZrCl2	MoS2	0.98	12 000	0.43	1.1
	WSSe	MoSSe	1.40	5500	0.37	1.0
Electrons	PtTe2	CdI2	0.30	9600	0.17	1.3
	GaO	GaS	1.56	7200	0.32	1.0
	NiS2	CdI2	0.58	6000	0.29	1.5
	RuTe2	WTe2	0.64	4600	1.55	8.5

To put the numbers in table 11 to scale, we consider the well studied example of MoS₂. For this material we obtain an electron mobility of 240 cm² V⁻¹ s⁻¹ while a full ab initio calculation found a phonon-limited mobility of 400 cm² V⁻¹ s⁻¹ (in good agreement with experiments on hBN encapsulated MoS₂ [148]), with the acoustic phonon contribution corresponding to a mobility of 1000 cm² V⁻¹ s⁻¹. Similarly, for the series MX₂ (M  =  W, Mo, X  =  S, Se), we calculate room-temperature electron mobilities between 200 cm² V⁻¹ s⁻¹ and 400 cm² V⁻¹ s⁻¹, which are all within 50% of the ab initio results [149]. Presumably, as in the case for MoS₂, the good quantitative agreement is partly a result of error cancellation between an overestimated acoustic phonon scattering and the neglect of optical phonon scattering. Importantly, however, the relative ordering of the mobilities of the four MX₂ monolayers is correctly predicted by equation (17) for all but one pair (MoS₂ and WSe₂) out of the six pairs. These results illustrate that equation (17) should only be used for 'order of magnitude' estimates of the mobility but that relative comparisons of mobilities in different materials are probably reliable.

3.2.3. Magnetic properties

Recently, a single layer of CrI₃ was reported to exhibit ferromagnetic order with a Curie temperature of 45 K [12]. This comprises the first example of a pure 2D material exhibiting magnetic order and there is currently an intense search for new 2D materials with magnetic order persisting above room temperature [150–152].

For 2D materials, magnetic order will only persist at finite temperatures in the presence of magnetic anisotropy (MA). Indeed, by virtue of the Mermin–Wagner theorem, magnetic order is impossible in 2D unless the rotational symmetry of the spins is broken [153]. A finite MA with an out of plane easy axis breaks the assumption behind the Mermin–Wagner theorem and makes magnetic order possible at finite temperature. The critical temperature for magnetic order in 2D materials will thus have a strong dependence on the anisotropy.

The MA originates from spin–orbit coupling and is here defined as the energy difference between in-plane and out-of-plane orientation of the magnetic moments, see equation (4). With our definition, a negative MA corresponds to an out-of-plane easy axis. We note that most of the materials in the C2DB are nearly isotropic in-plane. Consequently, if the easy axis lies in the plane, the spins will exhibit an approximate in-plane rotational symmetry, which prohibits magnetic order at finite temperatures. Since spin–orbit coupling becomes large for heavy elements, we generally expect to find larger MA for materials containing heavier elements. In general the magnitude of the MA is small. For example, for a monolayer of CrI₃ with a Curie temperature of 45 K [12] we find a MA of –0.85 meV per Cr atom in agreement with previous calculations [154]. Although small, the MA is, however, crucial for magnetic order to emerge at finite temperature.

In figure 24 we show the magnitude of the magnetic anisotropy (red triangles) and the magnetic moment per metal atom (blue triangles) averaged over all materials with a given chemical composition. The plot is based on data for around 1200 materials in the medium to high stability categories (see table 2) out of which around 350 are magnetic. It is interesting to note that while the magnetic moment is mainly determined by the metal atom, the MA depends strongly on the non-metal atom. For example, the halides (Cl, Br, I) generally exhibit much larger MAs than the chalcogenides (S, Se, Te). Overall, iodine (I) stands out as the most significant element for a large MA while the 3d metals Cr, Mn, Fe, Co are the most important elements for a large magnetic moment. Since the MA is driven by spin–orbit coupling (SOC) and the spin is mainly located on the metal atom, one would expect a large MA to correlate with a heavy metal atom. However, it is clear from the figure that it is not essential that the spin-carrying metal atom should also host the large SOC. For example, we find large MA for several 3d metal-iodides despite of the relatively weak SOC on the 3d metals. This shows that the MA is governed by a rather complex interplay between the spins, orbital hybridisation and crystal field.

Zoom In Zoom Out Reset image size

A selection of materials predicted to have high overall stability (see section 2.4.1) and high out-of-plane magnetic anisotropy (${\rm MA}< -0.3$ meV/magnetic atom) is listed in table 12. We find several semiconductors with anisotropies comparable to CrI3 and some metals with higher values. If we also look at materials with medium overall stability, we find semiconductors with anisotropies up to 2 meV/atom. It is likely that some of these materials will have Curie temperatures similar to, or even higher than, CrI₃.

Table 12. Selection of magnetic materials with a negative MA per magnetic atom. The prototype, the magnetic moment of the magnetic atom, the energy gap calculated with PBE xc-functional and the magnetic state are also shown. The experimentally synthesised ferromagnetic monolayer CrI₃ is highlighted.

Formula	Prototype	Magnetic moment ($\mu_{\mathrm{B}}$)	PBE gap (eV)	MA (meV/atom)	Magnetic state	$\Delta H_{\mathrm{hull}}$ (eV/atom)
OsI3	BiI3	0.9	0.0	−3.17	FM	0.18
CrTe	FeSe	2.6	0.0	−1.85	AFM	0.15
FeCl3	BiI3	1.0	0.01	−1.84	FM	−0.08
FeTe	FeSe	1.9	0.0	−1.06	FM	0.08
MnTe2	CdI2	2.7	0.0	−0.94	FM	−0.10
FeBr3	BiI3	1.0	0.04	−0.88	FM	−0.04
CrI$_{\mathbf{3}}$	BiI$_{\mathbf{3}}$	3.0	0.90	$-{\bf 0.85}$	FM	$-{\bf 0.21}$
MnTe	FeSe	3.8	0.69	−0.75	AFM	−0.15
NiO	PbSe	1.1	0.0	−0.53	FM	0.05
FeBrO	FeOCl	1.1	0.0	−0.47	FM	−0.05
CrISe	FeOCl	3.0	0.0	−0.45	FM	−0.10
MnSe2	CdI2	2.8	0.0	−0.40	FM	−0.18
CrIS	FeOCl	3.0	0.35	−0.36	FM	−0.10
MnO2	CdI2	3.0	1.13	−0.36	FM	0.02
VCl3	BiI3	2.0	0.0	−0.35	FM	−0.01
MnSe	FeSe	3.7	0.90	−0.31	AFM	−0.20
CrSe	FeSe	2.0	0.0	−0.31	AFM	0.12

In addition to the MA, the critical temperature depends sensitively on the magnetic exchange couplings. We are presently developing a workflow for systematic calculation of exchange coupling constants, which will allow us to estimate the Curie temperature of all the magnetically ordered 2D materials. The database contains several 2D materials with anti-ferromagnetic order. As a note of caution, we mention that the magnetic interactions in AFM materials typically arise from the super-exchange mechanism, which is poorly described by PBE and requires a careful verification using a PBE  +  U scheme [155].

3.2.4. Plasmons

The unique optical properties of 2D materials make them highly interesting as building blocks for nanophotonic applications [156, 157]. Many of these applications involve electron rich components which can capture, focus, and manipulate light via plasmons or plasmon-polaritons. Graphene sheets can host plasmons that are long lived, can be easily tuned via electrostatic or chemical doping, and can be used to confine light to extremely small volumes [158]. However, due to the limited charge carrier density achievable in graphene, its plasmons are limited to the mid-infrared regime. Here we show that some metallic monolayers support plasmons with significantly higher energies than graphene and could potentially push 2D plasmonics into the optical regime.

Figure 25 (left) shows the plasmon dispersion for monolayer NbS₂ in the MoS₂ crystal structure. The effect of interband transitions on the plasmon is significant as can be seen by comparison to the pure intraband plasmon ($\hbar\omega_\mathrm{P}$). The true plasmon energies are obtained from the poles of the (long wave length limit) dielectric function including the interband transitions, $\newcommand{\e}{{\rm e}} \epsilon = 1 + 2 \pi q (\alpha^{2{\rm D}, \mathrm{intra}} + \alpha^{2{\rm D}, \mathrm{inter}})$, yielding $\omega^\mathrm{true}_\mathrm{P} = \omega_\mathrm{P, 2D}q^{1/2}[1 + 2\pi q \alpha^{2{\rm D}, \mathrm{inter}}(\omega_\mathrm{P}^\mathrm{true})]^{-1/2}$. For simplicity we ignore the frequency dependence of the interband polarisability, i.e. we set $\alpha^{2{\rm D}, \mathrm{inter}}(\omega_\mathrm{P}^\mathrm{true})\approx$$\alpha^{2{\rm D}, \mathrm{inter}}(\omega=0)$, which should be valid for small plasmon energies (far from the onset of interband transitions). The validity of this approximation is confirmed by comparing to the full ab initio calculations of Andersen et al (blue dots) which include the full q- and ω-dependence [159]. The figure shows that interband screening generally reduces the plasmon energy and becomes increasingly important for larger q.

Zoom In Zoom Out Reset image size

Figure 25 (right) shows the in-plane averaged true plasmon energy of all metals in the C2DB plotted against the intraband plasmon energy at a fixed plasmon wavelength of $\lambda_\mathrm{P}=50$ nm (corresponding to q₀ at the dashed vertical line in the left panel). For comparison, the plasmon energy of freestanding graphene at this wavelength and for the highest achievable doping level ($E_F=\pm 0.5$ eV relative to the Dirac point) is around 0.4 eV. The data points are colored according to the overall stability level as defined in section 2.4.1. Table 13 shows a selection of the 2D metals with 'high' overall stability (see section 2.4.1) and large plasmon frequencies. We briefly note the interesting band structures of the metals in the FeOCl prototype (not shown) which exhibits band gaps above or below the partially occupied metallic band(s), which is likely to lead to reduced losses in plasmonic applications [160]. A detailed study of the plasmonic properties of the lead candidate materials will be published elsewhere. However, from figure 25 (right) it is already clear that several of the 2D metals have plasmon energies around 1 eV at $\lambda_\mathrm{P}=50$ nm, which significantly exceeds the plasmon energies in highly doped graphene.

Table 13. Selection of 2D metals with high plasmon energies $\omega_\mathrm{P}^\mathrm{true}$ for a plasmon wavelength of $\lambda_\mathrm{P}=50$ nm. The interband screening $\alpha^{2{\rm D}, \mathrm{inter}}$ at $\omega=0$ and the in-plane averaged 2D plasma frequency $\omega_\mathrm{P, 2D}$, which are required to reproduce $\omega_\mathrm{P}^\mathrm{true}$, are also reported.

Prototype	Formula	Magnetic state	$\omega_\mathrm{P}^\mathrm{true}$ (eV)	$\alpha^{2{\rm D}, \mathrm{inter}}$ (Å)	$\omega_\mathrm{P, 2D}$ (eV Å0.5)
TiS3	TaSe3	NM	0.99	12.54	12.48
FeOCl	PdClS	NM	0.93	4.13	9.52
FeOCl	NiClS	NM	0.90	5.60	9.66
CdI2	TaS2	NM	0.82	5.59	8.79
FeOCl	ZrIS	NM	0.75	7.68	8.43
CdI2	NbS2	NM	0.73	8.2	8.42
FeOCl	ZrClS	NM	0.73	13.6	9.35
TiS3	TaS3	NM	0.73	34.22	12.44
PbSe	NiO	FM	0.72	2.9	7.16

3.2.5. Excitons

Two-dimensional materials generally exhibit pronounced many-body effects due to weak screening and strong quantum confinement. In particular, exciton binding energies in monolayers are typically an order of magnitude larger than in the corresponding layered bulk phase and it is absolutely crucial to include excitonic effects in order to reproduce experimental absorption spectra.

The electronic screening is characterised by the in-plane 2D polarisability, see equation (12). For a strictly 2D insulator, the screened Coulomb potential can be written as $\newcommand{\e}{{\rm e}} W^{\mathrm{2D}}(q)=v_c^{\mathrm{2D}}(q)/\epsilon^{\mathrm{2D}}(q)$ with $\newcommand{\e}{{\rm e}} \epsilon^{\mathrm{2D}}(q)=1+2\pi\alpha^{\mathrm{2D}} q$ and $v_c^{\mathrm{2D}}(q)=2\pi/q$ is the 2D Fourier transform of the Coulomb interaction. The q-dependence of $\newcommand{\e}{{\rm e}} \epsilon^{\mathrm{2D}}$ indicates that the screening is non-local, i.e. it cannot be represented by a q-independent dielectric constant, and that Coulomb interactions tend to be weakly screened at large distances (small q vectors) [119, 161, 162]. This is in sharp contrast to the case of 3D semiconductors/insulators where screening is most effective at large distances where its effect can be accounted for by a q-independent dielectric constant. For a two-band model with isotropic parabolic bands, the excitons can be modeled by a 2D Hydrogen-like (Mott–Wannier) Hamiltonian where the Coulomb interaction is replaced by $\newcommand{\e}{{\rm e}} W=1/\epsilon r$ and the electron mass is replaced by a reduced excitonic mass $\mu_{\mathrm{ex}}$ derived from the curvature of conduction and valence bands. It has previously been shown that the excitonic Rydberg series of a 2D semiconductor can be accurately reproduced by this model if the dielectric constant, , is obtained by averaging $\newcommand{\e}{{\rm e}} \epsilon^{\mathrm{2D}}(q)$ over the extent of the exciton in reciprocal space [163]. For the lowest exciton (n  =  1), the binding energy can then be expressed as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:E_B} E_{\mathrm{B}}=\frac{8\mu_{\mathrm{ex}}}{(1+\sqrt{1+32\pi\alpha^{\mathrm{2D}}\mu_{\mathrm{ex}}/3})^2}. \nonumber \end{align} \tag{ 18 }$$

It has furthermore been demonstrated that this expression gives excellent agreement with a numerical solution of the Mott–Wannier model employing the full q-dependent dielectric function, $\newcommand{\e}{{\rm e}} \epsilon^{\mathrm{2D}}(q)=1+2\pi\alpha^{\mathrm{2D}} q$, for 51 transition metal dichalcogenides [163]. We note that the previous calculations were based on LDA and we generally find that the PBE values for $\alpha^{\mathrm{2D}}$ obtained in the present work are 10–20% smaller compared with LDA.

In figure 26, we compare the binding energy of the lowest exciton obtained from BSE-PBE with G₀W₀ scissors operator and the Mott–Wannier model equation (18), respectively. Results are shown for the 194 non-magnetic semiconductors out of the total set of $\sim 250$ materials for which BSE calculations have been performed. We focus on the optically active zero-momentum excitons and compute the exciton masses by evaluating the curvature of the band energies at the direct gap, see section 2.11. For anisotropic materials we average the heavy and light exciton masses as well as the x and y components of the polarisability, $\alpha^{\mathrm{2D}}$, to generate input parameters for the isotropic model equation (18).

Zoom In Zoom Out Reset image size

Although a clear correlation with the BSE results is observed, it is also evident that the Mott–Wannier model can produce significant errors. The mean absolute deviation between BSE and the model is 0.28 eV for all materials and 0.20 eV for the subset of transition metal dichalcogenides (TMDCs). Furthermore, the Mott–Wannier model seems to overestimate $E_{\mathrm{B}}$ for more strongly bound excitons while the opposite trend is seen for weakly bound excitons. As explained below these trends are a consequence of systematic errors in the Mott–Wannier model which can be traced to two distinct sources.

• (i)  
  Weak screening: If $\alpha^{\mathrm{2D}}$ is small (on the order of 1 Å), the exciton becomes strongly localised and the orbital character of the states comprising the exciton plays a significant role. In general, the Mott–Wannier model tends to overestimate the exciton binding energy in this case as can be seen from the relatively large deviation of points with model binding energies  >2.0 eV in figure 26. The overestimated binding energy results from the homogeneous electron and hole distributions implicitly assumed in the Mott–Wannier model. In reality, the short range variation of the electron and hole distributions is determined by the shape of the conduction and valence band states. In general these will differ leading to a reduced spatial overlap of the electron and hole and thus a lower Coulomb interaction. For example, SrCl₂ in the CdI₂ prototype ($\alpha^{\mathrm{2D}}=0.68$ Å) has a BSE binding energy of 2.1 eV and a model binding energy of 3.4 eV. From the PDOS of this material (see the C2DB webpage) it is evident that the electron and hole are mainly residing on the Sr and Cl atoms, respectively.
• (ii)  
  Breakdown of the parabolic band approximation: Materials with small band gaps often exhibit hyperbolic rather than parabolic band structures in the vicinity of the band gap. This typically happens in materials with small band gaps such as BSb in the BN prototype. In figure 26 these materials can be identified as the cluster of points with model binding energies  <0.25 eV and BSE binding energies  >0.25 eV. A similar situation occurs if the conduction and valence bands flatten out away from the band gap region. In both of these cases, the excitons tend to be delocalised over a larger area in the Brillouin zone than predicted by the parabolic band approximation of the Mott–Wannier model. Typically, such delocalisation will result in larger binding energies than predicted by the model. For example, FeI₂ in the CdI₂ prototype exhibits shallow band minima in a ring around the Γ-point and has a BSE binding energy of 1.1 eV and a model binding energy of 0.5 eV because the model assumes that the exciton will be located in the vicinity of the shallow minimum (and thus more delocalised in real space). A detailed inspection reveals that such break down of the parabolic band approximation is responsible for most of the cases where the model underestimates the binding energy.

Other sources of errors come from contributions to the exciton from higher/lower lying bands, i.e. break down of the two-band approximation, and anisotropic exciton masses not explicitly accounted for by equation (18).

Based on this comprehensive and unbiased assessment of the Mott–Wannier model, we conclude that while the model can be useful for understanding trends and qualitative properties of excitons, its quantitative accuracy is rather limited when applied to a broad set of materials without any parameter tuning. For quantitative estimates $\alpha^{\mathrm{2D}}$ should not be too small (certainly not less than 2 Å) and the the validity of the effective mass approximation should be carefully checked by inspection of the band structure.

It has been argued that there should exist a robust and universal scaling between the exciton binding energy and the quasiparticle band gap of 2D semiconductors, namely $E_\mathrm{B} \approx E_{\mathrm{gap}}/4$ [164]. This scaling relation was deduced empirically based on BSE-GW calculations for around 20 monolayers and explained from equation (18) and the relation $\alpha^{\mathrm{2D}}\propto 1/E_{\mathrm{gap}}$ from $\mathbf k \cdot \mathbf p$ perturbation theory. Another work observed a similar trend [165] but explained it from the $1/E_{{\rm gap}}$ dependence of the exciton effective mass expected from $\mathbf k \cdot \mathbf p$ perturbation theory. Based on our results we can completely refute the latter explanation (see figure 22 (right)). In figure 26 (right) we show the exciton binding energy plotted versus the direct PBE and ${\rm G_0W_0}$ band gaps, respectively. While there is a correlation, it is by no means as clear as found in [164].