# Blue phosphorene bilayer is a two-dimensional metal - and an unambiguous
classification scheme for buckled hexagonal bilayers

Jessica Arcudia  Departamento de Física Aplicada, Centro de Investigación y de
Estudios Avanzados, Unidad Mérida, Mérida, Yucatán, México  Roman Kempt
Technische Universität Dresden, Fakultät für Chemie und Lebensmittelchemie,
Bergstraße 66c, 01062 Dresden, Germany  Miguel Eduardo Cifuentes-Quintal
Departamento de Física Aplicada, Centro de Investigación y de Estudios
Avanzados, Unidad Mérida, Mérida, Yucatán, México  Thomas Heine
<EMAIL_ADDRESS>Technische Universität Dresden, Fakultät für Chemie
und Lebensmittelchemie, Bergstraße 66c, 01062 Dresden, Germany  Helmholtz
Zentrum Dresden-Rossendorf, Leipzig Research Branch, Permoserstr. 15, 04318
Leipzig, Germany  Department of Chemistry, Yonsei University, Seodaemun-gu,
Seoul 120-749, Republic of Korea  Gabriel Merino<EMAIL_ADDRESS>Departamento de Física Aplicada, Centro de Investigación y de Estudios
Avanzados, Unidad Mérida, Mérida, Yucatán, México

###### Abstract

High-level first-principles computations predict blue phosphorene bilayer to
be a two-dimensional metal. This structure has not been considered before and
was identified by employing a block-diagram scheme that yields the complete
set of five high-symmetry stacking configurations of buckled honeycomb layers,
and allows their unambiguous classification. We show that all of these
stacking configurations are stable or at least metastable configurations both
for blue phosphorene and gray arsenene bilayers. For blue phosphorene, the
most stable stacking configuration has not yet been reported, and surprisingly
it is metallic, while all other arrangements are indirect band gap
semiconductors. As it is impossible to interchange the stacking configurations
by translations, all of them should be experimentally accessible via the
transfer of monolayers. The metallic character of blue phosphorene bilayer is
caused by its short interlayer distance of 3.01 Å and offers the exceptional
possibility to design single elemental all-phosphorus transistors.

two-dimensional metal, phosphorene, arsenene, two-dimensional materials,
stackings

††preprint: APS/123-QED

##  I Introduction

Since the first exfoliation of graphene and the discovery of its remarkable
properties [1], many other 2D materials came in the focus of interest,
including its isoelectronic hexagonal congener boron nitride [2], and the wide
variety of transition metal dichalcogenides [3, 4, 5]. However, it was not
before the exfoliation black phosphorus [6, 7], that the attention for 2D
pnictogens emerged. Black phosphorus is the most stable allotrope of
phosphorus and crystallizes in an orthorhombic structure. But a single layer
of phosphorus, so-called phosphorene, can also crystallize in other forms. The
honeycomb structure predicted by Zhu and Tománek [8] is just 2 meV/atom higher
in energy than the black counterpart, and given that its band gap value is
slightly above the photon energy of visible blue light, it was named blue
phosphorene (in the following bluP). This prediction was materialized just a
couple of years later by Zhang et al., who successfully synthesized monolayer
(ML) bluP by epitaxial growth on an Au(111) substrate [9]. The next pnictogen
of interest is arsenic, mostly found in its bulk form as gray arsenic, and
whose layered rhombohedral structure makes it an excellent candidate for
exfoliation. The first studies on the stability and the properties of a single
arsenic layer, i.e., gray arsenene (in the following grAs), predicted by Kamal
and Ezawa [10] and Kou et al. [11], encouraged experimental groups to
exfoliate few-layers arsenic [12].

Both bluP and its arsenic congener grAs have a structure similar to graphene,
which is a hexagonal lattice with their atoms alternatingly being displaced
out of the 2D plane (Figure 1a). These monolayers show exciting electronic
properties with an indirect band gap in the range of 1.5-2.0 eV [8, 9, 10,
12]. BluP is a $p$-type semiconductor with high carrier mobility [13]. On the
other hand, grAs could be used for transistors or mechanical sensors, due to
its indirect-direct band gap transitions [11, 14], semiconductor-metal
transitions [8, 10], and topological phase transitions under strain [15, 16].

It is well-known that interlayer interactions can significantly alter the
properties of 2D materials, most strikingly discussed recently for
superconductivity in bilayer (BL) graphene [17, 18], but also for band gap
nature and valleytronics in transition-metal dichalcogenides (for a review,
see Ref. [19]), or for metal-insulator transitions in noble metal
chalcogenides (for a review, see Ref. [20]). The properties crucially depend
on the stacking type and twist angle, which can be controlled using various
transfer techniques [21].

Only a few studies addressed interlayer effects in bluP and grAs, and no
details on the stacking order have been reported experimentally up to now [9,
12]. Given the corrugation of the monolayer, more than the high-symmetry AA
and AB stacking orders are expected. In how many ways is it possible to stack
them? Even for the simplest case of two buckled honeycomb layers, this
question has not been answered yet. Herein, we suggest a new approach based on
a simple graphical analysis, which we will call a ‘block diagram’ to ascertain
it.

Using block diagrams, we identify all possible high-symmetry stacking
configurations of BL buckled honeycomb lattices, which are expected to yield
all low-energy forms for bluP and grAs bilayers.

Employing first principles calculations (density-functional theory (DFT), the
Random Phase Approximation (RPA), and single-particle Green’s function
approach $G_{0}W_{0}$), we calculated their structures, thermodynamic and
kinetic stabilities, and explore their electronic structures. We discovered
that the lowest-energy bluP BL was not yet reported to date and, surprisingly,
it is metallic.

##  II Results and Discussion

Because both bluP and grAs crystallize in the trigonal lattice, we can
describe a single layer in the subperiodic layer group $\mathbf{P}$3 (#65)
with the unique sites $1a=(0,0,z)$, $1b=(1/3,2/3,z_{1})$ and
$1c=(2/3,1/3,z_{2})$, where two sites are occupied with $z_{1}=-z_{2}$ [22].

The lattice vectors ($\vec{a}_{1}$ and $\vec{a}_{2}$) delimiting the two-
dimensional lattice are defined as
$\vec{a}_{1}=\frac{1}{2}a\hat{x}-\frac{\sqrt{3}}{2}a\hat{y}$,
$\vec{a_{2}}=\frac{1}{2}a\hat{x}+\frac{\sqrt{3}}{2}a\hat{y}$, where $a$ is the
lattice constant. Then the unit cells of any trigonal monolayer with a
buckling can be represented as in Figure 1b, where the arrows indicate the
out-of-plane displacements. We conveniently label these forms as A1, A-1, B1,
and B-1, where the negative sign in the subscripts denotes a displacement
change of the atoms with respect to the plane (up or down). These four
configurations are symmetry-equivalent and can be used for the construction of
the BL forms. For a better understanding, we design a simple graphical method,
which we call ‘block diagram’ to facilitate the visualization.

Figure 1: (a) Front and side view of a buckled honeycomb monolayer. The unit
cell with its lattice vectors and high symmetry sites 1a, 1b, and 1c. (b) Top:
Four different representations of a unit cell (arrows indicate the out of
plane displacement of the atoms), bottom: block diagrams to represent each
monolayer.

This block diagram is divided into three parts, and each of them specifies a
site of the unit cell (1a, 1b, or 1c). The atoms occupy two of the three sites
and the direction of the arrows is related to the buckling. So, to build a
bilayer, one puts another block on the top. In this way, we identify the total
number of different high-symmetry forms of corrugated honeycomb bilayer as
five, where the stacking sequence of two of them is of AA-type, and for the
remaining three it is AB-type (see block diagrams in Figure 2 and 3). These
five stacking configurations are reported for silicene, whose layers are
covalently bound [23] .

For BL bluP, only four different stacking configurations have been studied to
date [24, 25]. First principles calculations indicated the A1A1 configuration
to be the most stable one, followed by A1B1, A-1B1, and A1A-1. For grAs, there
is an on-going debate whether the A1B-1 or the A1A1 form is the most stable
one [26, 25, 27, 28].

Figure 2: Top view of an AA-type bilayer. The block diagram, atomistic
structure, and phonon and electronic band structures of its two
configurations: A1A1, and A1A-1 Figure 3:  Top view of an AB-type bilayer. The
block diagram, atomistic structure, and phonon and electronic band structures
of its three configurations: A1B1, A-1B1, and A1B-1. The second minimum of
A1B-1 (2-A1B-1) is included. The most stable stacking configuration is
indicated by a frame.

For all stable and metastable stacking configurations discussed in this work,
full geometry optimizations yield similar structural parameters, independent
on the density-functional, the choice of London dispersion correction scheme,
the orbital representation (local basis functions vs. plane waves), or the
underlying code (see Methods for details). In all but one case a single local
minimum per stacking configuration is found. Only for the A1B-1 configuration
of bluP we found two local minima: one minimum corresponds to a structure
where the interlayer distance $d$ is small (3.01 Å) and the structure shows
small corrugation, corresponding to a smaller buckling height $\Delta z$. We
call this structure A1B-1. The second structure (2-A1B-1) shows a larger
interlayer distance of 4.93 Å, corresponding to weakly interacting layers, and
larger buckling. Both for bluP and grAs BL the closer interlayer distance in
A1B-1 is accompanied by a lattice constant increase of 0.1 Å and by a buckling
reduction of $\Delta z\approx 0.07$ Å. The remainder of the structures has
almost identical lattice constants and buckling heights as the monolayer
($a$=3.26 Å and $\Delta z$=1.24 Å, Table 1) [9, 29]. Structural parameters for
all stacking configurations of BL bluP and grAs are summarized in Table 1.

Table 1: The lattice constant ($a$), buckling height ($\Delta z$), and interlayer distance ($d$), computed at the PBE+MBD level for blue phosphorene (bluP) and gray arsenene (grAs) bilayers. | bluP BL | grAs BL
---|---|---
System |  $a$ (Å) |  $\Delta z$ (Å) |  $d$ (Å) |  $a$ (Å) |  $\Delta z$ (Å) |  $d$ (Å)
A1B-1 | 3.36 | 1.17 | 3.01 | 3.69 | 1.35 | 3.66
A1A1 | 3.26 | 1.24 | 4.66 | 3.60 | 1.40 | 4.58
A1B1 | 3.27 | 1.24 | 4.68 | 3.61 | 1.40 | 4.38
2-A1B-1 | 3.27 | 1.24 | 4.93 | - | - | -
A-1B1 | 3.26 | 1.24 | 5.36 | 3.59 | 1.40 | 5.43
A1A-1 | 3.26 | 1.24 | 5.42 | 3.59 | 1.40 | 5.48
Monolayer | 3.26 | 1.24 | - | 3.60 | 1.40 | -

All investigated BL systems are significantly more stable than their ML
counterparts and are unlikely to exfoliate without severe intrusion. For both
bluP and grAs, the A1B-1 stacking configuration was found to be the most
stable one. In both cases, this structure has distinct features making it
quite different compared to the other configurations: It has the smallest
interlayer distance and the smallest corrugation. The interlayer binding
energy, defined as $E_{ib}=\frac{E_{BL}-2E_{ML}}{N},$ exceeds 180 meV per atom
in both cases (Table 2). It is important to note that local and hybrid
density-functionals give different stacking orders for the stacking
configurations, which result in disagreement on the most stable form.
Substantiation at the RPA level, independent if starting from a PBE or PBE0
calculation, result in the same stacking order and clearly identify A1B-1 to
be the most stable form.

For bluP A1B-1 BL, we show an abrupt change in buckling height $\Delta z$ and
the distance between layers $d$ upon variation of lattice constant $a$,
corresponding to the phase transition from the metallic configuration A1B-1 to
the semiconducting configuration 2-A1B-1 (see Fig. S1a-b in the Supplementary
Material). In Fig. S1c we plot the buckling height $\Delta z$ vs. the energy
difference $\Delta E$, indicating that A1B-1 is energetically favorable as it
is less corrugated. Fig. S1d indicates that A1B-1 is extremely sensitive to
variations in the interlayer distance $d$. 2-A1B-1 shows weaker interlayer
interactions and is located in the second, very shallow minimum.

While for BL bluP the configuration A1B-1 is by far the most stable one, for
BL grAs two configurations are energetically competitive: A1A1 and A1B1 are
only less than 26 meV per atom higher in energy. Relative energies with
respect to the monolayers and with respect to the most stable bilayer forms
are given in Table 2.

Table 2: The relative energies ($\Delta E$) and the interlayer binding energy ($E_{ib}$) for all hexagonal blu phosphorene and gray arsenene stacking configurations, computed at the RPA+rSE@PBE0 level (including the ZPE correction). Units are in meV/atom. | bluP BL | grAs BL
---|---|---
System |  $\Delta$E | Eib |  $\Delta$E | Eib
A1B-1 | 0.0 | -180.2 | 0.0 | -185.2
A1A1 | 83.9 | -96.2 | 25.7 | -159.5
A1B1 | 86.8 | -93.4 | 16.1 | -169.0
2-A1B-1 | 104.1 | -76.1 | - | -
A-1B1 | 126.2 | -54.0 | 108.9 | -76.3
A1A-1 | 128.4 | -51.8 | 96.5 | -88.6

The phonon dispersion of the bluP and grAs MLs show the typical features of 2D
materials, with two linear and one out-of-plane quadratic branches of the
acoustic modes, and a clear energetic separation of optical and acoustic
branches (Figure S2). For bluP and grAs BL systems (Figure 2 and 3), the
phonon dispersion significantly differs between the most stable A1B-1 and the
other structures: While for all systems the three acoustic modes are in the
same energy range as three low-energy optical modes (emerging from the 2nd
bilayer system and indicating weaker interactions), for A1B-1 the stronger
interlayer interaction lifts the low-energy optical modes towards higher
energies and these branches show significantly less dispersion compared to the
other stacking configurations. These shifted vibrational modes are of Eg and
Ag symmetry and Raman active, and thus could serve for the characterization of
this BL, similar as it has been demonstrated for graphene and silicene
experimentally [30]. For the high-energy arrangements A-1B1 and A1A-1 of both
bluP and grAs BLs we found small imaginary frequencies which are due to the
limitations of the numerical approach used in the phonon calculations.

To substantiate the electronic band gaps, the electronic structures have been
recalculated within the quasiparticle approach, using the single-shot
$G_{0}W_{0}$ approximation on top of PBE Kohn-Sham bands. Spin-orbit coupling
(SOC) was included in all calculations for the grAs systems (see Table 3).
Besides the well-known band gap underestimation, the limitations of local DFT
calculations also include the incorrect location of the valence band maximum
(VBM) for bluP (see also Figure S3) [31].

Most importantly, both DFT and $G_{0}W_{0}$ identify bluP BL in the most
stable configuration A1B-1 as a metal. All other systems are indirect band gap
semiconductors, with the $G_{0}W_{0}$ band gaps being considerably larger
compared to those calculated by PBE. For BL bluP, electronic band gaps range
from 2.39 to 2.82 eV, values that are lower than those of the ML (3.25 eV).
For grAs, the A1B-1 low-energy form has a remarkably small band gap of 0.58
eV, while the other systems range from 1.38 eV to 2.04 eV, thus being somewhat
narrower than that of the ML (2.30 eV). For all semiconducting bluP and grAs
BLs we found an interesting competition of a parabolic band climaxing at the
$\Gamma$ point and a Mexican hat structure around $\Gamma$ point to become the
VBM. The latter is obtained for bluP BL in A1A1 and A1B1 configuration
(somewhat less pronounced for grAs BL). A similar Mexican hat structure has
been reported for GaSe and other III-VI monolayers, and the small dispersion
may be useful for small FET structures and could give rise to Landau levels
[32, 33].

Table 3: Indirect band gaps (Egap) of ML and BL stacking configurations of bluP and grAs calculated at PBE and $G_{0}W_{0}$ levels of theory. Units are in eV. Spin-orbit coupling (SOC) was considered for the grAs systems. | Egap-bluP | Egap-grAs
---|---|---
System | PBE | $G_{0}W_{0}$ | PBE-SOC |  $G_{0}W_{0}$-SOC
A1B-1 | - | - | 0.12 | 0.58
A1A1 | 1.23 | 2.39 | 0.84 | 1.64
A1B1 | 1.25 | 2.40 | 0.63 | 1.38
2-A1B-1 | 1.41 | 2.59 | - | -
A-1B1 | 1.58 | 2.80 | 1.27 | 2.04
A1A-1 | 1.60 | 2.82 | 1.22 | 1.99
Monolayer | 1.96 | 3.25 | 1.43 | 2.30

Finally, we return to the remarkable 2D metal A1B-1 bluP. The band structure
(Figure 3) includes metallic bands crossing the Fermi level near the $\Gamma$
point, characterizing it as conventional metal. In addition, there is a basin
of charge carriers close to $M$ point. This interesting band structure
suggests anisotropic electronic properties that shall not be further explored
at this stage.

##  III Summary and Outlook

We present a scheme to identify and to label all symmetrically distinct
stacking configurations of corrugated honeycomb bilayers and investigate the
resulting structures blue phosphorene (bluP) and gray arsenene (grAs). We
discovered a new, and a the same time the most stable, configuration of bluP,
which has a small interlayer distance, quenched corrugation and is a new
member of the exclusive group of two-dimensional metals.

Besides the most stable A1B-1 form, also A1A1 and A1B1 stacking configurations
could be obtained experimentally by layer transfer techniques, as their large
interlayer interaction and symmetry constraints prevent interconversion or
relaxation into other structures. Lower-symmetry configurations as they are
known for BL graphene are unlikely here due to the surface corrugation. This
holds both for bluP and grAs. Except for the metallic bluP A1B-1 form, all
investigated bluP and grAs bilayers are indirect band gap semiconductors and
resemble the Mexican hat-type feature of the electronic bands near the valence
band maximum which have been intensively discussed for III-VI-group mono- and
multilayer systems.

BluP shows a semiconductor-metal transition when going from the ML to the most
stable BL. Such transitions have been found for other 2D materials (e.g. noble
metal chalcogenides [34] and GeP3 [35]) and could be used to design single-
material transistors with low Schottky barrier between the electrode and
semiconducting scattering region [36]. In a similar vein, we think it is a
remarkable challenge to construct the first single-element transistor based on
blue phosphorus, with its metallic bilayer configuration serving as electrode
material.

##  IV Methods

To identify all potentially existing isomeric structures, we started geometry
optimization with different starting structures, varying interlayer distances
and buckling. All geometries have been fully optimized by means of DFT within
the framework of the projector-augmented wave (PAW) method [37, 38] as
implemented in the Vienna ab initio simulation package (VASP 5.4.4) [39, 40].

The valence states were expanded in plane-waves with an energy cutoff of 400
eV. The Perdew, Burke, and Ernzerhof (PBE) exchange-correlation functional was
employed [41], with London dispersion interactions taken into account by the
many-body dispersion (MBD) correction as suggested by Tkatchenko and coworkers
[42]. For comparison, we also used the Barone and Adamo’s hybrid PBE0
functional [43, 44] which does not affect the geometries of the bilayers. The
first Brillouin zone was sampled using a $\Gamma$-centered k-grid of 15x15x1
and 42x42x1 k-points for hexagonal phosphorene and arsenene bilayers,
respectively, and Gaussian smearing of 0.10 eV. These numerical parameters
were chosen to ensure a convergence criterion of 10-8 eV in total energy, and
10-3 eV/Å in atomic forces. A vacuum space of 25Å was considered in order to
avoid interactions with the repeated images. The lattice parameters of each
system were obtained by direct minimization of the total energy, with the
atomic positions fully optimized until the interatomic forces were less than
10-3 eV/Å.

The same level was employed for the phonon calculations, for which we used the
small displacement method as implemented in the Phon code [45]. The force
constant matrix was computed using central differences within atomic
displacements of 0.02Å in 9x9 supercells.

Band structures have been calculated with $G_{0}W_{0}$ [46, 47], considering
spin-orbit coupling (SOC) for all arsenic systems.

Relative stabilities have been calculated from single-point energies of
optimized structures on the grounds of PBE-MBD and PBE0-MBD, which give
different stacking orders and required calculations beyond DFT (see Table S1).
We employed the Random-Phase Approximation (RPA) with the renormalized single
excitation correction (rSE) [48] based on Kohn-Sham orbitals from the PBE and
PBE0 level as implemented in FHI-AIMS [49] on tight tier 1 numeric atom-
centered orbitals with added auxiliary diffuse basis functions on 12x12x1
$k$-grids to determine the most stable form. The RPA stacking orders are
independent on the choice of the underlying density-functional (for
comparison, see Supplementary Material), herein we discuss only the
RPA+rSE$@$PBE0 variant. Relative energies $\Delta E$ (in meV per atom) are
given with respect to the most stable stacking configuration, and interlayer
binding energies $E_{ib}$ are defined with respect to the monolayers.

##  V Acknowledgements

This work was funded in Mexico by Cinvestav (Grant SEP-Cinvestav-2018-57). The
CGSTIC (Xiuhcoatl) is acknowledged for allocation of computational resources.
J. A. thanks Conacyt for her Ph.D. fellowship. We thank the Center for
Information Services and High Performance Computing (ZIH) at TU Dresden for
generous allocation of computer time. We thank BMBF (NobleNEMS) for financial
support. We thank Andras Kis (EPFL) for fruitful discussions on 2D
heterostructure manufacturing techniques.

## References

  * Novoselov  _et al._ [2004] K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, Y. Zhang, S. V. Dubonos, I. V. Grigorieva, and A. A. Firsov, Science 306, 666 (2004).
  * Britnell  _et al._ [2012] L. Britnell, R. V. Gorbachev, R. Jalil, B. D. Belle, F. Schedin, M. I. Katsnelson, L. Eaves, S. V. Morozov, A. S. Mayorov, N. M. Peres, _et al._ , Nano Lett. 12, 1707 (2012).
  * Miró  _et al._ [2014a] P. Miró, M. Audiffred, and T. Heine, Chem. Soc. Rev. 43, 6537 (2014a).
  * An  _et al._ [2018] S.-J. An, Y. H. Kim, C. Lee, D. Y. Park, and M. S. Jeong, Sci. Rep. 8, 1 (2018).
  * Haastrup  _et al._ [2018] S. Haastrup, M. Strange, M. Pandey, T. Deilmann, P. S. Schmidt, N. F. Hinsche, M. N. Gjerding, D. Torelli, P. M. Larsen, A. C. Riis-Jensen, _et al._ , 2D Mater. 5, 042002 (2018).
  * Liu  _et al._ [2014] H. Liu, A. T. Neal, Z. Zhu, Z. Luo, X. Xu, D. Tománek, and P. D. Ye, ACS Nano 8, 4033 (2014).
  * Li  _et al._ [2014] L. Li, Y. Yu, G. J. Ye, Q. Ge, X. Ou, H. Wu, D. Feng, X. H. Chen, and Y. Zhang, Nat. Nanotechnol. 9, 372 (2014).
  * Zhu and Tománek [2014] Z. Zhu and D. Tománek, Phys. Rev. Lett. 112, 176802 (2014).
  * Zhang  _et al._ [2016] J. L. Zhang, S. Zhao, C. Han, Z. Wang, S. Zhong, S. Sun, R. Guo, X. Zhou, C. D. Gu, K. D. Yuan, _et al._ , Nano Lett. 16, 4903 (2016).
  * Kamal and Ezawa [2015] C. Kamal and M. Ezawa, Phys. Rev. B 91, 085423 (2015).
  * Kou  _et al._ [2015] L. Kou, Y. Ma, X. Tan, T. Frauenheim, A. Du, and S. Smith, J. Phys. Chem. C 119, 6918 (2015).
  * Beladi-Mousavi  _et al._ [2019] S. M. Beladi-Mousavi, A. M. Pourrahimi, Z. Sofer, and M. Pumera, Adv. Funct. Mater. 29, 1807004 (2019).
  * Xiao  _et al._ [2015] J. Xiao, M. Long, X. Zhang, J. Ouyang, H. Xu, and Y. Gao, Sci. Rep. 5, 9961 (2015).
  * Zhang  _et al._ [2015a] S. Zhang, Z. Yan, Y. Li, Z. Chen, and H. Zeng, Angew. Chem. Int. Ed. 127, 3155 (2015a).
  * Zhang  _et al._ [2015b] H. Zhang, Y. Ma, and Z. Chen, Nanoscale 7, 19152 (2015b).
  * Wang  _et al._ [2016] Y.-p. Wang, C.-w. Zhang, W.-x. Ji, R.-w. Zhang, P. Li, P.-j. Wang, M.-j. Ren, X.-l. Chen, and M. Yuan, J. Phys. D Appl. Phys. 49, 055305 (2016).
  * Cao  _et al._ [2018a] Y. Cao, V. Fatemi, S. Fang, K. Watanabe, T. Taniguchi, E. Kaxiras, and P. Jarillo-Herrero, Nature 556, 43 (2018a).
  * Cao  _et al._ [2018b] Y. Cao, V. Fatemi, A. Demir, S. Fang, S. L. Tomarken, J. Y. Luo, J. D. Sanchez-Yamagishi, K. Watanabe, T. Taniguchi, E. Kaxiras, _et al._ , Nature 556, 80 (2018b).
  * Kuc and Heine [2015] A. Kuc and T. Heine, Chem. Soc. Rev. 44, 2603 (2015).
  * Kempt  _et al._ [2020] R. Kempt, A. Kuc, and T. Heine, Angew. Chem. Int. Ed. (2020).
  * Novoselov  _et al._ [2016] K. Novoselov, A. Mishchenko, A. Carvalho, and A. C. Neto, Science 353 (2016).
  * Aroyo  _et al._ [2006] M. I. Aroyo, A. Kirov, C. Capillas, J. Perez-Mato, and H. Wondratschek, Acta Cryst. A62, 115 (2006).
  * Padilha and Pontes [2015] J. Padilha and R. B. Pontes, J. Phys. Chem. C 119, 3818 (2015).
  * Pontes  _et al._ [2018] R. B. Pontes, R. H. Miwa, A. J. da Silva, A. Fazzio, and J. E. Padilha, Phys. Rev. B 97, 235419 (2018).
  * Ahn  _et al._ [2018] J. Ahn, I. Hong, Y. Kwon, R. C. Clay, L. Shulenburger, H. Shin, and A. Benali, Phys. Rev. B 98, 085429 (2018).
  * Cao  _et al._ [2015] H. Cao, Z. Yu, and P. Lu, Superlattice Microst. 86, 501 (2015).
  * Kecik  _et al._ [2016] D. Kecik, E. Durgun, and S. Ciraci, Phys. Rev. B 94, 205409 (2016).
  * Kadioglu  _et al._ [2018] Y. Kadioglu, J. A. Santana, H. D. Özaydin, F. Ersan, O. Ü. Aktürk, E. Aktürk, and F. A. Reboredo, J. Chem. Phys. 148, 214706 (2018).
  * Wang  _et al._ [2018] B.-J. Wang, X.-H. Li, R. Zhao, X.-L. Cai, W.-Y. Yu, W.-B. Li, Z.-S. Liu, L.-W. Zhang, and S.-H. Ke, J. Phys. Chem. A 6, 8923 (2018).
  * Liang  _et al._ [2017] L. Liang, J. Zhang, B. G. Sumpter, Q.-H. Tan, P.-H. Tan, and V. Meunier, ACS Nano 11, 11777 (2017).
  * Ikiyanat  _et al._ [2019] F. Ikiyanat, E. Torun, R. T. Senger, and H. Sahin, Phys. Rev. B 100, 125423 (2019).
  * Aziza  _et al._ [2018] Z. B. Aziza, V. Zólyomi, H. Henck, D. Pierucci, M. G. Silly, J. Avila, S. J. Magorrian, J. Chaste, C. Chen, M. Yoon, _et al._ , Phys. Rev. B 98, 115405 (2018).
  * Kuc  _et al._ [2017] A. Kuc, T. Cusati, E. Dib, A. F. Oliveira, A. Fortunelli, G. Iannaccone, T. Heine, and G. Fiori, Adv. Electron. Mater. 3, 1600399 (2017).
  * Miró  _et al._ [2014b] P. Miró, M. Ghorbani-Asl, and T. Heine, Angew. Chem. Int. Ed. 53, 3015 (2014b).
  * Jing  _et al._ [2017] Y. Jing, Y. Ma, Y. Li, and T. Heine, Nano Lett. 17, 1833 (2017).
  * Ghorbani-Asl  _et al._ [2016] M. Ghorbani-Asl, A. Kuc, P. Miró, and T. Heine, Adv. Mater. 28, 853 (2016).
  * Blöchl [1994] P. E. Blöchl, Phys. Rev. B 50, 17953 (1994).
  * Kresse and Joubert [1999] G. Kresse and D. Joubert, Phys. Rev. B 59, 1758 (1999).
  * Kresse and Furthmüller [1996a] G. Kresse and J. Furthmüller, Comput. Mater. Sci. 6, 15 (1996a).
  * Kresse and Furthmüller [1996b] G. Kresse and J. Furthmüller, Phys. Rev. B 54, 11169 (1996b).
  * Perdew  _et al._ [1996a] J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996a).
  * Tkatchenko  _et al._ [2012] A. Tkatchenko, R. A. DiStasio Jr, R. Car, and M. Scheffler, Phys. Rev. Lett. 108, 236402 (2012).
  * Perdew  _et al._ [1996b] J. P. Perdew, M. Ernzerhof, and K. Burke, J. Chem. Phys. 105, 9982 (1996b).
  * Adamo and Barone [1998] C. Adamo and V. Barone, Chem. Phys. Lett. 298, 113 (1998).
  * Alfè [2009] D. Alfè, Comput. Phys. Commun. 180, 2622 (2009).
  * Hedin [1965] L. Hedin, Phys. Rev. 139, A796 (1965).
  * Shishkin and G. [2006] M. Shishkin and K. G., Phys. Rev. B 74, 035101 (2006).
  * Ren  _et al._ [2013] X. Ren, P. Rinke, G. E. Scuseria, and M. Scheffler, Phys. Rev. B 88, 035120 (2013).
  * Blum  _et al._ [2009] V. Blum, R. Gehrke, F. Hanke, P. Havu, V. Havu, X. Ren, K. Reuter, and M. Scheffler, Comput. Phys. Commun. 180, 2175 (2009).