# Electronic and Excitonic Properties of Semi-Hydrogenated Borophene Sheets

Mohammad Ali Mohebpour  Shobair Mohammadi Mozvashi  Computational Nanophysics
Laboratory (CNL), Department of Physics, University of Guilan, P. O. Box
41335-1914, Rasht, Iran.  Sahar Izadi Vishkayi  School of Physics, Institute
for Research in Fundamental Sciences (IPM), P. O. Box 19395-5531, Tehran,
Iran.  Meysam Bagheri Tagani<EMAIL_ADDRESS>Computational Nanophysics
Laboratory (CNL), Department of Physics, University of Guilan, P. O. Box
41335-1914, Rasht, Iran.

###### Abstract

Borophene has triggered a surge of interest due to its outstanding properties
including mechanical flexibility, polymorphism, and opto-electrical
anisotropy. Very recently, a novel semi-hydrogenated borophene, called
$\alpha^{\prime}$-4H, was synthesized in large-scale freestanding samples,
which exhibits excellent air-stability and semiconducting nature. Herein,
using the density functional theory (DFT) and many-body perturbation theory
(MBPT), we investigate the electronic and excitonic optical properties of
$\alpha^{\prime}$-4H borophene. The DFT results reveal that by breaking the
mirror symmetry and increasing the buckling height of pure $\alpha^{\prime}$
borophene, hydrogenation causes an orbital hybridization and opens an indirect
band gap of 1.49 eV in $\alpha^{\prime}$-4H borophene. This value is corrected
to be 1.98, 2.23, and 2.52 eV under the G0W0, GW0, and GW levels of theory,
respectively. The optical spectrum achieved from solving the Bethe-Salpeter
equation shows an optical band gap of $\sim$ 2.40 eV, which corresponds to a
strongly bound and stable bright exciton with a binding energy of $\sim$ 1.18
eV. More importantly, the excitonic states are robust against tension up to
10%, where the monolayer is dynamically stable. We also design and study the
bilayer $\alpha^{\prime}$-4H borophene with different stackings. For the weak
van der Waals interactions between the layers, the bilayer can preserve most
of the structural and electronic properties of the monolayer. Our study
exposes the underlying physics behind the structural, electronic, and optical
properties of $\alpha^{\prime}$-4H borophene and suggests it as a very
promising candidate for flexible optoelectronic applications.

##  I Introduction

Two-dimensional (2D) boron, also known as borophene, is recently considered as
a game-changing player of future electronics because of its exceptional
properties including polymorphism, electronic and optical anisotropy,
ultrahigh thermal conductance, optical transparency, and phonon-mediated
superconductivity Mannix  _et al._ (2018); Liu  _et al._ (2020a); Penev  _et
al._ (2016); Vishkayi and Tagani (2017); Vajary  _et al._ (2018); Vishkayi and
Tagani (2018). Moreover, short B-B bonds and low atomic mass leads to very
outstanding mechanical characteristics in borophene Mannix  _et al._ (2015).
Numerous studies have studied the applications of borophene in rechargeable
ion batteries, gas storages, catalysts, and sensors Arabha  _et al._ (2020);
Ayodhya and Veerabhadram (2020). Several recent studies have focused on the
use of borophene in flexible electronics, wearable devices, and medicines,
owing to its flexible nature, low mass density, and non-toxicity Hou  _et al._
(2020); Xie  _et al._ (2020).

Boron has three valance electrons and, to compensate its octet completion
deficiency, it tends to make highly-delocalized multicenter bonds, resulting
in a vast structural diversity and polymorphism Entwistle and Marder (2002).
Several motifs of borophene have been proposed and synthesized up to now,
among which are $\alpha$, $\alpha^{\prime}$, $\beta$, and $\chi$ phases, which
have trigonal structures with different numbers and positions of hexagon holes
(HHs) Kong  _et al._ (2019); Powar and Pandav (2019). Different phases of
borophene are defined in terms of parameter $\nu=n/N$, where $n$ and $N$ are
the numbers of HHs and boron atoms in a unit cell, respectively. Most of the
physical properties of borophenes strongly depend on this morphology variable,
thus, there is a wide set of options to reach a borophene phase with desired
properties. Furthermore, unlike other 2D materials, borophene undergoes a
structural phase transition, rather than direct fracturing under tensile
strain Mannix  _et al._ (2015).

Due to its atomic arrangement, borophene cannot be exfoliated mechanically.
Therefore, other methods such as chemical vapor deposition, molecular beam
epitaxy, and liquid phase exfoliation should be used Bhattacharyya  _et al._
(2017). Moreover, because of its electron deficiency, borophene cannot be
stable in a free-standing form and support of a metal substrate is required
Suehara  _et al._ (2010). Several successful syntheses have been reported on
Ag, Cu, Ni, and Au substrates Wu  _et al._ (2019); Kiraly  _et al._ (2019);
Zhong  _et al._ (2017). It is worth mentioning that the Ag, Cu, and Ni
substrates donate electrons to, and Au substrate withdraws electrons from
borophene, hence, the borophenes grow flat on the former substrates and
buckled on the latter one Mannix  _et al._ (2015); Liu  _et al._ (2013).

Because of its high reactivity, most of boron nanostructures suffer from rapid
oxidation and presence of other ambient contaminations Mannix  _et al._
(2015). Surface modification is found to be a capable tool for enhancing the
stability and electronic properties of nanostructures Ryder  _et al._ (2016);
Johns and Hersam (2013). Specifically, hydrogenation is reported widely to
have a major impact on the stability of borophene and other two-dimensional
materials Li  _et al._ (2021); Mohebpour  _et al._ (2020).

Recently, Hou et al Hou  _et al._ (2020) have prepared a novel semi-
hydrogenated borophene phase in large quantities by in-situ decomposition of
NaBH4 under the exposure of hydrogen. This borophene was named
$\alpha^{\prime}$-4H because it is similar to $\alpha^{\prime}$ borophene,
with 4 adsorbed hydrogen atoms in each unit cell. The $\alpha^{\prime}$-4H
borophene has three major advantages which make it very attractive. First, it
has high dynamical stability in freestanding form: by withdrawing electrons
from the boron atoms, hydrogenation acts similar to an Au substrate, resulting
a stable buckled structure. Second, it is stable in the air for a long time:
the screening of hydrogen passivates the surface of borophene and prevents it
from reactions with oxygen or other ambient atoms. Third, unlike other phases,
$\alpha^{\prime}$-4H borophene has an energy band gap of $\sim$2.5 eV, which
brings hopes for its applications in optoelectronics as a semiconductor.

One of the most important drivers of recent advances in condensed matter
physics is the feedback Ping-Pong between theory and experiment. High-quality
theoretical and computational investigations can truly guide the experiments
towards novel materials with desired properties. Moreover, discussing the
underlying physics of quantum systems improves the physical insights and
enhances the theory for future advances. To have a small contribution to this
trend, in our previous work, we systematically investigated the mechanical
flexibility and strength of $\alpha^{\prime}$-4H borophene by first-principle
calculations Mozvashi  _et al._ (2021). We also demonstrated dynamical
stability, mechanical anisotropy, high strain compliance, and tunable band gap
for this monolayer.

In this paper, we investigate the electronic and excitonic optical properties
of $\alpha^{\prime}$-4H borophene using ab-initio many-body calculations. By
density functional theory (DFT) calculations, we initially discuss the
microscopic aspects in which the metallic $\alpha^{\prime}$ borophene turns
into the semiconducting $\alpha^{\prime}$-4H. We show how hydrogenation leads
to an orbital hybridization and opens a wide band gap of 1.49 eV at the PBE
level of theory.

Many-body effects (electron-electron and electron-hole interactions) are
significantly enhanced in low-dimensional systems Cudazzo  _et al._ (2010);
Yang  _et al._ (2009). This is mainly due to the reduced Coulomb screening and
geometrical confinement. Hence, we expect to see the critical influence of
electron-electron and electron-hole interactions on the electronic and optical
properties of $\alpha^{\prime}$-4H borophene. First, we calculate the
quasiparticle band gap to be 1.98, 2.23, and 2.52 eV at the G0W0, GW0, and GW
levels of many-body calculations, respectively. Then, we solve the Bethe-
Salpeter equation on top of quasiparticle energies (the so-called GW+BSE) to
obtain the optical spectrum. It shows that the direct optical excitation is
characterized by an optical band gap of $\sim$ 2.40 eV, which corresponds to a
strongly bound and stable Frenkel exciton with a binding energy of $\sim$ 1.18
eV, making $\alpha^{\prime}$-4H borophene a desirable candidate for exploring
optical applications.

Finally, we design and investigate the bilayer $\alpha^{\prime}$-4H borophene
with different stackings and discuss that for the weak van der Waals
interaction between the layers, the bilayer can preserve most of structural
and electronic properties of the monolayer and can be used for similar
applications. Our results bring insights about physical properties of the
monolayer and bilayer $\alpha^{\prime}$-4H borophene and suggest it as a
promising material for flexible optoelectronic applications.

##  II Theoretical Background

We start the density functional theory (DFT) calculation to find the ground-
state eigenvalues ($E_{nk}$) and eigenfunctions ($\Phi_{nk}$) by solving the
Kohn-Sham equation given below:

$(T+V_{n-e}+V_{H}+V_{xc})\Phi_{nk}=E_{nk}\Phi_{nk},$ (1)

where T is the kinetic energy operator, $V_{n-e}$ the Coulomb potential of the
nuclei, $V_{H}$ the Hartree potential, and $V_{xc}$ the exchange-correlation
potential. The quasiparticle (QP) energies ($E_{nk}^{QP}$) are determined
within the GW approximation as a first-order perturbative correction to the
DFT single-particle energies ($E_{nk}$) Shishkin and Kresse (2007); Hybertsen
and Louie (1986):

$\displaystyle
E_{nk}^{QP,1}=Re[\bra{\Phi_{nk}}T+V_{n-e}+V_{H}+\Sigma_{xc}(G,W;E_{nk})$ (2)
$\displaystyle\ket{\Phi_{nk}}]=E_{nk}+\braket{\Phi_{nk}}{\Sigma^{GW}-V_{xc}}{\Phi_{nk}},$

where $\Sigma=iG_{0}W_{0}$ is the self-energy operator. This equation is
solved non-self-consistently. Therefore, it is the simplest level of theory
(the so-called G0W0) and computationally the most efficient one. It calculates
the QP energies from a single iteration by neglecting all off-diagonal
elements of the self-energy and employing a Taylor expansion of the self-
energy around the DFT energies. The noninteracting Green’s function ($G_{0}$)
is directly calculated from single-particle orbitals and energies as:

$G_{0}(r,r^{\prime},\omega)=\sum\frac{\Phi_{nk}(r)\Phi_{nk}^{*}(r^{\prime})}{\omega-
E_{nk}-i\eta~{}sgn(E_{F}-E_{nk})},$ (3)

where $E_{F}$ is the Fermi energy and $\eta$ a positive infinitesimal
quantity. The noninteracting screened Coulomb potential ($W_{0}$) is obtained
by the inverse dielectric function ($\varepsilon^{-1}$) and bare Coulomb
interaction ($\upsilon$) as below Liu  _et al._ (2020b):

$W_{0}(r,r^{\prime},\omega^{\prime})=\int
dr^{\prime\prime}\varepsilon^{-1}(r,r^{\prime\prime},\omega^{\prime})\upsilon(r^{\prime\prime},r^{\prime}).$
(4)

At higher levels of theory (i.e. GiW0 and GiWi), the updated QP energies
($E_{nk}^{QP,i+1}$) are achieved from the QP energies of the previous
iteration as follows:

$\displaystyle
E_{nk}^{QP,i+1}=E_{nk}^{QP,i}+Z_{nk}Re[\bra{\Phi_{nk}}T+V_{n-e}$ (5)
$\displaystyle+V_{H}+\Sigma_{xc}(G,W;E_{nk})\ket{\Phi_{nk}}-E_{nk}^{QP,i}],$

where the renormalization factor ($Z_{nk}$) is given by Karlicky and Otyepka
(2013); Shahrokhi and Leonard (2016); Shishkin and Kresse (2006):

$Z_{nk}=(1-Re\bra{\Phi_{nk}}\frac{\partial}{\partial\omega}\Sigma(\omega)_{|_{E_{nk}^{QP,i}}}\ket{\Phi_{nk}})^{-1}.$
(6)

At the GW0 level of theory, the Green’s function ($G$) is iteratively updated
while the screened Coulomb potential ($W$) is kept fixed at the initial DFT
value. At the GW level, both of them are reevaluated at each iteration ($i$).
Here, the exact (interacting) Green’s function is linked to its non-
interacting version with the Dyson equation given as, $G=G_{0}+G_{0}\Sigma G$
Blase  _et al._ (2020).

The imaginary part of the macroscopic dielectric function is obtained by
summation over the empty conduction bands within the long-wavelength limit
($q$ $\rightarrow$ 0) Gajdoš  _et al._ (2006):

$\begin{split}Im\varepsilon_{\alpha\beta}(\omega)=\frac{4\pi^{2}e^{2}}{\Omega}\lim_{q\to
0}\frac{1}{|q|^{2}}\sum_{c,v,k}2w_{k}\delta(\varepsilon_{ck}-\varepsilon_{vk}-\omega)\\\
\times\braket{u_{ck+e_{\alpha}q}}{u_{vk}}\braket{u_{ck+e_{\beta}q}}{u_{vk}}^{*},\end{split}$
(7)

where $q$ is the Bloch vector of the incident wave, $w_{k}$ the $k$-point
weight, $\Omega$ the volume of the unit cell, and $u_{ck}$ the cell periodic
part of the wave function. The vectors $e_{\alpha}$ are unit vectors for the
three Cartesian directions. The indices $c$ and $v$ correspond to the
conduction and valence band states, respectively. The electron-hole excited
states are characterized by the expansion Shahrokhi and Leonard (2017):

$\ket{S}=\sum_{c}^{elec}\sum_{v}^{hole}\sum_{k}A_{cvk}^{S}\ket{cvk},$ (8)

where $A^{S}$ is the amplitude of a free electron-hole pair configuration
composed of the electron state $\ket{ck}$ and the hole state $\ket{vk}$. It is
achieved from the diagonalization of the following excitonic equation Rohlfing
and Louie (1998, 2000):

$\begin{split}(E_{ck}^{QP}-E_{vk}^{QP})A_{cvk}^{S}+\sum_{c^{\prime}v^{\prime}k^{\prime}}\braket{cvk}{\Xi^{e-h}}{c^{\prime}v^{\prime}k^{\prime}}A_{c^{\prime}v^{\prime}k^{\prime}}^{S}\\\
=\Omega^{s}A^{s}_{cvk},\end{split}$ (9)

where $E_{ck}^{QP}(E_{vk}^{QP})$ denotes the QP eigenvalues of valance
(conduction) band at a specific $k$-point, kernel $\Xi^{e-h}$ describes the
screened interaction between excited electrons and holes, and $\Omega^{S}$ is
the excitation energy.

##  III Computational Details

The first-principles calculations were performed in the framework of DFT,
using the generalized gradient approximation developed by Perdew-Burke-
Ernzerhof (GGA-PBE) and the projector augmented waves (PAW) pseudopotentials
implemented in the Quantum ESPRESSO package Giannozzi  _et al._ (2009). The
energy cutoff for the plane-wave basis set was set to be 60 Ry. The Brillouin
zone was integrated with a 13$\times$ 13 $\times$ 1 $\Gamma$-centered
$k$-point mesh. The atomic positions and lattice parameters were fully relaxed
until a force convergence of $10^{-3}$ eV/Å was achieved. A vacuum space of
$L_{z}$ = 20 Å was considered to avoid spurious interactions.

The many-body perturbation calculations were carried out using the GW
approximation. The eigenvalues and eigenstates achieved from DFT-PBE served as
input to calculate the QP energies. The calculations were performed under
different levels of self-consistency including single-shot (i.e. G0W0) and
partially self-consistent (i.e. GW0 and GW). A compromised set of input
parameters (e.g. 9 $\times$ 9 $\times$ 1 $k$-point mesh, $L_{z}$ = 25 Å, 40
virtual bands, and 4 iterations) was utilized for our GW calculations.

The QP band structure was interpolated at the G0W0 level using the maximally
localized Wannier functions (MLWFs) implemented in the WANNIER90 code Mostofi
_et al._ (2014). Here, the $sp^{3}$ and $s$ orbitals of B and H atoms were
respectively chosen for the initial projections. The excitonic optical
properties were investigated by solving the Bethe-Salpeter equation (BSE) on
top of GW eigenvalues and wave functions, using the Tamm-Dancoff approximation
(TDA). This approximation takes only the resonant part of the BS Hamiltonian
into account. The 10 highest valence bands and the 10 lowest conduction bands
were considered as the basis for the excitonic eigenstates.

##  IV Results and Discussion

###  IV.1 Structural and electronic properties

The structural configurations of fully relaxed $\alpha^{\prime}$ and
$\alpha^{\prime}$-4H borophenes are displayed in FIG. 1 (a, b) and the
corresponding structural parameters are listed in Table S1. As can be seen,
the pure $\alpha^{\prime}$ borophene is nearly flat ($\Delta$ = 0.37 Å) and
isotropic, with a hexagon hole concentration of $\nu~{}=~{}1/8$. The
hydrogenation leads to an increase in the buckling in $\alpha^{\prime}$-4H
borophene ($\Delta$ = 0.88 Å), owing to the attraction between B and H atoms.
The lattice constants of the two monolayers are similar (5.05 Å) but the
isotropy of the lattice breaks after hydrogenation. As a result, the bond
length increases from isotropic 1.68 Å to anisotropic 1.69 and 1.77 Å.

FIG. 1 (a, b) also indicate the electron density cut planes of
$\alpha^{\prime}$ and $\alpha^{\prime}$-4H borophenes. In $\alpha^{\prime}$
borophene, we can see a “sea of electrons”, which is apparently the cause of
its metallic nature. However, in $\alpha^{\prime}$-4H borophene, an electron
accumulation takes place around H atoms, which lowers the in-plane carrier
distribution and decreases the electronic conductivity.

As mentioned above, the $\alpha^{\prime}$-4H borophene was synthesized
recently Hou  _et al._ (2020). To validate our model with the experiment, we
performed simulations of Raman spectrum and scanning tunneling microscopy
(STM). FIG. 1 (c) represents the calculated Raman spectrum of
$\alpha^{\prime}$-4H borophene. We can see three peaks at about 700, 1100, and
2200 cm-1. The first two peaks are representatives of $E_{g}$ and
$A_{1g}+E_{g}$ modes of the B-B cluster, and the last peak is associated with
the combined modes of the B-H bonds. We also indicated the simulated STM image
of $\alpha^{\prime}$-4H borophene in FIG. 1 (d), which is very similar to the
experimentally obtained image. Our calculated Raman spectrum and simulated STM
image support the experimental data very well Hou  _et al._ (2020), therefore,
we can assure that our simulated model is completely similar to the
synthesized $\alpha^{\prime}$-4H borophene.

Figure 1: (a, b) Structural configurations and electron density cut planes of
$\alpha^{\prime}$ and $\alpha^{\prime}$-4H borophenes. (c) Simulated Raman
spectrum and (d) simulated STM image of $\alpha^{\prime}$-4H borophene. Figure
2: (a, b) Partial density of states (PDOS) of $\alpha^{\prime}$ and
$\alpha^{\prime}$-4H borophenes, showing the significant contributing orbitals
of the labeled atomic sites shown in the insets. Full results of PDOS are also
available in FIG. S1. (c, d) Orbital decomposed band structures of
$\alpha^{\prime}$ and $\alpha^{\prime}$-4H borophenes. Both PDOS and band
structure were calculated at the GGA-PBE level of theory.

FIG. 2 summarizes the structural and electronic properties of
$\alpha^{\prime}$ and $\alpha^{\prime}$-4H borophenes. As shown in the inset
of FIG. 2 (a, b), the atomic sites in $\alpha^{\prime}$ borophene are
classified into 6- and 5-bonded boron atoms, labeled with ⓐ and ⓑ,
respectively. There are two ⓐ and six ⓑ sites in each unit cell and we can see
both inversion and mirror symmetry in $\alpha^{\prime}$ borophene. After
hydrogenation, four ⓑ sites break a B-B bond and make a B-H bond instead,
which leads to novel sites labeled with ⓒ. As a result, hydrogenation breaks
the mirror symmetry and only preserves the inversion symmetry in
$\alpha^{\prime}$-4H borophene. These changes directly affect the electronic
states and should be considered as one of the origins of different properties
in the two borophenes.

FIG. 2 also represents the partial density of states (PDOS) and orbital
decomposed band structures of $\alpha^{\prime}$ and $\alpha^{\prime}$-4H
borophenes at the GGA-PBE level. This information gives valuable insights into
the contribution of each atomic site and orbital in the electronic properties.
As it is clear, in $\alpha^{\prime}$ borophene (left panel), $p_{z}$ orbitals
of the ⓑ sites are dominant. Given that $\alpha^{\prime}$ borophene is almost
flat, $p_{z}$ orbitals are perpendicular to neighboring $p_{x(y)}$ orbitals.
Specifically, no hybridization takes place between these orbitals. In the band
structure, we can see a triple energy degeneracy which is similar to a Dirac
point, pierced by a third band. The triple band touching Dirac fermions was
previously observed in other flat borophenes in experimental and theoretical
studies Ezawa (2017); Feng  _et al._ (2017).

In $\alpha^{\prime}$-4H borophene (right panel), hydrogenation pulls some of ⓑ
atoms, labeled with ⓒ, out of the plane and forms a buckled structure. As a
result, $p_{z}$ orbitals of ⓒ atoms would no longer be perpendicular to
neighboring $p_{x(y)}$ orbitals, which leads to hybrid bonds. This
hybridization can be seen very well in the PDOS and band structure. By
hydrogenation, ⓑ sites lose their dominance around the Fermi level. However, ⓐ
sites rather reserve their contribution, especially in the conduction band,
because hydrogenation do not affect them directly. Subsequently,
$\alpha^{\prime}$-4H borophene has dominated occupied ⓒ-$p_{y}$ states and
unoccupied ⓐ-$p_{z}$ states in the valance and conduction band edges,
respectively, which are separated with a relatively wide band gap.

These results are completely supported by the shapes of the wave functions at
the VBM and CBM as represented in FIG. S2. At the VBM, the wave function is
shaped like a dumbbell, which is distributed along the $y$-axis and centered
on ⓒ atoms, showing the contribution of ⓒ-$p_{y}$ orbitals. At the CBM, the
dumbbell-shaped wave function has components along both $x$-axis and $z$-axis,
showing the hybridization between ⓒ-$p_{x}$ and ⓐ-$p_{z}$ orbitals.

To highlight the role of buckling and orbital hybridization as a consequence,
we also calculated the band structure for a conceptual model of flat
$\alpha^{\prime}$-4H borophene. We re-optimized $\alpha^{\prime}$-4H borophene
with a fixed $z$-component, to have the same buckling as $\alpha^{\prime}$
borophene. The resulting band structure is available in FIG. S3 and shows no
band gap. Therefore, we can conclude that the buckling induced by ⓒ atoms is a
strong cause for opening the band gap. It should be noted that the role of
buckling in opening energy band gap has been widely considered in the
community Radha and Lambrecht (2020); Umeno  _et al._ (2019); Nijamudheen  _et
al._ (2015).

The band gap calculated at the GGA-PBE level is 1.49 eV, which is
underestimated by 1.0 eV compared to the experiment Hou  _et al._ (2020).
Using the HSE03 and HSE06 hybrid functionals and considering a fraction of the
exact exchange, the band gap is corrected to be 1.85 and 1.98 eV, respectively
(see FIG. S4). Eventually, at the PBE0 level, the band gap is found to be 2.54
eV, in excellent agreement with the experimental band gap measured by UV-V
spectrum (2.49 eV).

Although DFT gives adequate qualitative information about the electronic and
optical properties, it does not account for the electron-electron and
electron-hole interactions, which have a crucial role in 2D materials.
Therefore, further investigations should be performed using many-body
perturbation calculations.

The QP band gap predicted by many-body GW theory is very sensitive to input
parameters including the spin-orbit coupling (SOC), $k$-point sampling, vacuum
space, number of virtual bands, energy cutoff, and self-consistent iteration.
Therefore, it is mandatory to converge the band gap. We initially checked the
effects of the SOC interaction. Due to small atomic mass of boron, except a
slight downward shift in the electronic states, no noticeable spin splitting
was found in the energy levels, and the band gap remained almost constant
($\Delta$Eg $\simeq$ 0.001 eV). Hence, due to its excessive computational
costs and negligible impacts, this interaction is excluded from the rest of
the calculations.

Next, we performed a bunch of calculations to converge the QP band gap at the
G0W0 level as a function of the (a) $k$-point sampling of Brillouin zone using
40 virtual bands and $L_{z}$ = 20 Å; (b) vacuum space using 40 virtual bands
and 5 $\times$ 5 $\times$ 1 $k$-point mesh; and (c) the number of virtual
bands using $L_{z}$ = 20 Å and 5 $\times$ 5 $\times$ 1 $k$-point mesh as
illustrated in FIG. S6. We realized that a 9 $\times$ 9 $\times$ 1 $k$-point
mesh is enough for a satisfactory convergence within $\sim$ 0.05 eV. Such a
dense $k$-point mesh is also indispensable for an accurate description of the
optical gap and binding energy. The vacuum space plays a crucial role in the
determination of the QP band gap due to the long-range nature of screened
Coulomb interaction. Generally, it is imposed to be between 20 to 30 Å. To
ensure convergence within $\sim$ 0.05 eV, we set the vacuum space to be 25 Å,
which is the highest accuracy we could computationally afford. Our results
also verified that at least 40 virtual bands are required for our many-body GW
calculations. The corresponding band gap is essentially the same when 80 bands
are included (see FIG. S6). Such rapid convergence concerning the number of
unoccupied bands is mainly attributed to the absence of strongly localized
states characterized by flat bands. Moreover, we tested the convergence of the
QP band gap with respect to the number of iterations in partially self-
consistent GiW0 and GiWi as shown in FIG. S6. Obviously, four iterations are
sufficient to reach the convergence threshold of $\sim$ 0.001 eV. Also, the
energy cutoff for the exchange and correlation parts of the self-energy was
set to be 50 and 10 Ry, respectively.

To ensure the accuracy of our many-body GW calculations, we calculated the QP
band gaps of graphane (CH) and fluorographene (CF) monolayers. For example,
for graphane, the band gap was predicted to be 5.62, 5.81, and 6.17 eV under
the G0W0, GW0, and GW levels of theory, respectively, which are in excellent
agreement with 5.64, 5.89, and 6.28 eV reported in previous works Karlicky and
Otyepka (2013). Also, the band gaps obtained for fluorographene agree very
well with previous values as listed in Table S2.

After validation of the required accuracy, now we turn our attention to the GW
results. The general shape of the QP band structure at the G0W0 level is
similar to that obtained by the PBE functional, shown in FIG. S4. The VBM and
CBM still locate at the same points. The Fermi level is still closer to the
VBM, making $\alpha^{\prime}$-4H borophene a p-type semiconductor. The valence
and conduction band edges are symmetrically dispersed, which indicates the
presence of free charge carriers. However, the inclusion of electron-electron
interaction enlarges the band gap from 1.49 up to 1.98 eV. This value is
larger than that of the HSE03 and exactly equal to that obtained by the HSE06.
In other words, a similar accuracy could be achieved by the G0W0 as by the
HSE06. At the GW0 and GW levels, the QP band gap increases up to 2.23 and 2.52
eV, respectively, which is due to the presence of self-consistency. Overall,
one can say that the GW level offers an accuracy similar to the PBE0
functional, in excellent agreement with the experiment Hou  _et al._ (2020).
Also, the direct QP band gap at the $\Gamma$ point are 3.09, 3.31, and 3.58 eV
for G0W0, GW0, and GW, respectively.

###  IV.2 Optical properties

Excitonic effects play a crucial role in the optical properties of 2D
materials for their weak dielectric screening. Hence, to have a reliable
description of the optical properties, it is indispensable to take the
electron-hole Coulomb interaction into account. In this section, we aim to
present the excitonic optical spectra of $\alpha^{\prime}$-4H borophene
achieved from solving the Bethe-Salpeter equation on top of various levels of
the GW. To fully appreciate the many-body effects, we also calculated the
optical spectra within random-phase approximation (RPA) on top of DFT and GW
levels of theory, where the electron-hole interaction is neglected. FIG. 3
exhibits the imaginary part of the macroscopic dielectric function for light
polarized along the $x$-direction at different levels, namely DFT+RPA,
G0W0+RPA, and G0W0+BSE. As can be seen, the inclusion of electron-electron
correlation completely reshapes the optical spectrum and leads to a
significant blue shift ($\sim$ 1.4 eV), which is due to the self-energy
correction. On the other hand, the electron-hole correlation results in a
redshift ($\sim$ 0.92 eV) in the optical spectrum, owing to the cancellation
effect between self-energy correction and excitonic effects. The electron-hole
correlation also leads to an increase (decrease) of the oscillator strength of
the first (second) peak. However, its main impact is the appearance of some
bound excitons below the G0W0 direct band gap (3.09 eV) at the $\Gamma$ point,
which are missing in the G0W0+RPA spectrum. As the imaginary part of the
dielectric function is directly associated with the absorption spectrum, one
can say that the excitonic effects increase the light absorption in the
visible range (1.5 to 3.5 eV).

Figure 3: The imaginary part of the dielectric function of
$\alpha^{\prime}$-4H borophene for light polarized along the $x$-direction at
different levels of theory: DFT+RPA (e-e and e-h correlations neglected),
G0W0+RPA (e-e correlation included and e-h correlation neglected), and
G0W0+BSE (both e-e and e-h correlations included). Figure 4: The imaginary
part of the dielectric function of $\alpha^{\prime}$-4H borophene including
excitonic effects for light polarized along the $x$-direction under different
degrees of self-consistency in the GW. Four iterations were used to obtain
well-converged results from GiW0 and GiWi calculations. The QP (direct) band
gaps are shown with arrows.

The DFT+RPA predicts the absorption edge to be $\sim$ 1.48 eV, which is very
close to the PBE band gap (1.49 eV). The first two peaks at G0W0+RPA appear at
2.88 and 3.08 eV, respectively. They are near the G0W0 direct band gap (3.09
eV) at the $\Gamma$ point, which is due to the exclusion of excitonic effects.
Most importantly, the first peak at G0W0+BSE ($\sim$ 1.96 eV) is found to lie
between the PBE (1.49 eV) and G0W0 (1.98 eV) band gaps as it includes both the
QP correction and excitonic effects. This peak (the so-called optical gap)
originates from the direct transition (bright exciton) from the VBM to the
conduction band, and corresponds to a strongly bound Frenkel exciton with a
binding energy of $E_{b}$ $\simeq$ 1.13 eV. It is worth mentioning that the
exciton binding energy is determined by the difference between the QP direct
band gap and optical gap. The optical gap obtained from G0W0+BSE is $\sim$
0.53 eV smaller than the experimental gap achieved from photoluminescence
measurement (2.49 eV) Hou  _et al._ (2020). This underestimation is attributed
to the overscreening effects associated with building the susceptibility based
on an underestimated Kohn-Sham band gap. Indeed, the small Kohn-Sham band gap
leads to a spurious enhancement of the screening, and subsequently, to an
underestimated optical gap.

For the considerable dependence of the G0W0 on the Kohn-Sham starting point
and the exclusion of many-body exchange and correlation in the ground-state
eigenvalues, we examined the effects of self-consistency. FIG. 4 indicates the
optical spectra of $\alpha^{\prime}$-4H borophene for light polarized along
the $x$-direction under different degrees of self-consistency in the GW. As
can be seen, iteratively updating the QP eigenvalues in the Green’s function
results in a blue shift ($\sim$ 0.24 eV) in the optical spectrum because it
reduces the overestimation of the screening observed in the G0W0. At this
level of theory (i.e. GW0+BSE), the shape of the spectrum is relatively
preserved, and the optical gap is estimated to be $\sim$ 2.20 eV with an
exciton binding energy of $E_{b}$ $\simeq$ 1.11 eV. Therefore, one can say
that $E_{b}$ is independent of the self-consistency in Green’s function.
Moreover, the oscillator strength is only redistributed to higher energies
with no considerable change.

Updating the screened Coulomb potential leads to a blue shift ($\sim$ 0.20 eV)
as well, implying a larger self-energy correction in the GW compared to those
in the G0W0 and GW0. The highest level of theory (i.e. GW+BSE) predicts the
optical gap to be $\sim$ 2.40 eV, which is slightly smaller than the GW band
gap of $\sim$ 2.52 eV. This value agrees very well with the experimental gap
(2.49 eV) Hou  _et al._ (2020) because GW+BSE incorporates many-body effects
in the ground state eigenvalues and meaningfully reduces the dependence on the
DFT starting point using self-consistency. Here, the $E_{b}$ is obtained
$\sim$ 1.18 eV, suggesting relatively stronger cancellation effects. Such a
huge $E_{b}$ is a signature of a strongly bound Frenkel exciton, showing
ultrahigh stability of excitonic states against thermal dissociation, making
$\alpha^{\prime}$-4H borophene a very promising candidate for exploring room-
temperature optical applications. The final value is larger than those of
phosphorene (0.86 eV), MoS2 (0.96 eV), h-BN (0.80 eV), and many other
semiconducting monolayers Qiu  _et al._ (2013); Kolos and Karlický (2019); Xu
_et al._ (2017); Tran  _et al._ (2014), turning $\alpha^{\prime}$-4H borophene
into a highly competitive candidate. As the GW+BSE gives the best results, the
corresponding oscillator strength is shown in FIG. S7.

Because of the huge depolarization effect in 2D systems for the light
polarization perpendicular to the surface (E$\|z$) Yang  _et al._ (2007), we
only focus on the optical spectra for light polarization parallel to the
surface (E$\|x$, E$\|y$). As a consequence of the structural anisotropy in
$\alpha^{\prime}$-4H borophene, a different optical spectrum is obtained for
light polarized along the $y$-direction. As shown in FIG. S8, at all three
levels of theory, the shape of the spectrum is different from that for light
polarized along the $x$-direction (FIG. 4). Also, the first peak has weaker
intensity. However, the same optical gap is achieved. Such in-plane optical
anisotropy provides an additional degree of freedom for manufacturing
optoelectronic devices.

It is worth mentioning that the QP band gap and exciton binding energy
calculated at all three levels of theory (i.e. G0W0, GW0, and GW) show an
excellent agreement with the linear scaling law given as, $E_{b}$ =
0.21$~{}E_{g}$ \+ 0.40 Choi  _et al._ (2015). The standard deviations are less
than 0.05 eV, signifying well-converged results.

We also investigated the variation of band gap under biaxial tensile strain.
As can be seen from FIG. S9, by increasing the tensile strain up to 10%, the
QP and optical band gaps decrease linearly by $\sim$ 0.4 and 0.7 eV,
respectively. This is due to the reduction of buckling under tension. The
exciton binding energy also experiences a decrease of $\sim$ 0.3 eV, which
means the excitonic states are robust against tension up to 10%. It should be
noted, in our previous work, we had presented the dynamical stability of
$\alpha^{\prime}$-4H borophene under such tension by phonon dispersion,
stress-strain, and energy-strain curves Mozvashi  _et al._ (2021).

###  IV.3 Bilayer $\alpha^{\prime}$-4H borophene

Observations of atomic force microscopy (AFM) have shown that the average
thickness of a typical synthesized $\alpha^{\prime}$-4H borophene is about
0.78 nm Hou  _et al._ (2020), which is higher than the monolayer. To
understand the effects of thickness on the electronic properties of
$\alpha^{\prime}$-4H borophene, we studied the bilayer of this material. For
this, we firstly measured the vertical distance between layers by means of
variation of energy to thickness. As FIG. 5 (a) shows, the optimized thickness
of a bilayer is achieved 0.785 nm, which is in great agreement with
measurements of AFM. After gaining the interlayer distance, we re-optimized
the structure with AA stacking. The electron difference density in FIG. 5 (a &
b) shows that while a high electron extracting takes place by hydrogen atoms,
the two layers do not trade much electrons within themselves, which implies
the absence of strong covalent bonds and presence of weak van der Waals
interactions. This might be due to the screening of H atoms among the boron
layers. Because of the weak force between the layers, the lattice constant and
bond lengths of the bilayer do not have a significant difference with the
monolayer.

Figure 5: Bilayer $\alpha^{\prime}$-4H borophene. (a) Side view and electron
difference density map. (b) Variation of electron difference density with
buckling. (c) Band structure of AA stacking. (d) Variation of bilayer energy
with displacement of the upper layer with respect to AA stacking. (e)
Variation of bilayer band gap with displacement of the upper layer. (f)
Variation of band gap of AA stacking with thickness.

As shown in FIG. 5 (c), the AA-stacked bilayer $\alpha^{\prime}$-4H borophene
is an indirect semiconductor with a band gap of 1.4 eV. Most 2D materials
experience a topological band gap transition from their monolayer to bilayer
or multilayer forms. For instance, from monolayer to bilayer, graphene opens a
bang gap of zero to 0.25 eV Zhang  _et al._ (2009) and, on the contrary,
$\beta$ antimonene closes its band gap from 0.76 eV to zero Wang  _et al._
(2015). However, for the weak van der Waals interlayer interaction in bilayer
$\alpha^{\prime}$-4H borophene, it experiences a very small gap change. The
VBM is located at $\Gamma$ point, and the CBM is within the $\Gamma$-K path.
The bilayer $\alpha^{\prime}$-4H borophene is a p-type semiconductor because
VBM is very close to the Fermi energy. This is due to electron deficiency in
boron-based structures. In the VBM and CBM, the band structure dispersion is
parabolic, which exhibits high mobility of electrons and holes. This inclines
the high potential of $\alpha^{\prime}$-4H borophene in optoelectronic
applications, eVen for thicknesses more than monolayer.

To see the impact of stacking mode on the electronic properties of bilayer, we
slipped the upper layer along the lattice vector A1 as much as percentages of
lattice constant, a. The energy difference between different stackings with
respect to AA stacking is shown in FIG. 5 (d), which is less than 20 meV.
These stackings can form Moiré patterns in borophene with no need to external
strain or requirement of large superlattices. Our investigations show that
stacking mode does not have significant impact on the band gap of borophene.
The band gap remains indirect in all the stackings and, as shown in FIG. 5
(e), it fluctuate mildly around the band gap of the AA stacking. When the
displacement changes as half of the lattice constant (AB stacking), the band
gap becomes equal to the band gap of the AA stacking. Because different
stackings do not have very different band gaps, it is not possible to
determine the stacking mode of bilayer $\alpha^{\prime}$-4H borophene through
absorption spectrum.

To understand the dependence of electronic properties on the vertical distance
between the layers, variation of band gap with interlayer distance is shown in
FIG. 5 (f). Herein, with conserving the atomic arrangements, we tuned the
interlayer distance to have thicknesses between 7.7 to 8 Å. We can see a mild
increase of band gap with the increase of the thickness. In other words, a
high increase in the distance cancels the weak van der Waals interaction
between layers and renders the bilayer into two independent monolayers. On the
contrary, decrease of vertical interlayer distance increases the interactions
and leads to a mild decrease in the band gap. Our investigations show that
interlayer distance does not have any impact on the nature of indirect band
gap of $\alpha^{\prime}$-4H borophene. It can be concluded that for the weak
van der Waals interaction between the layers and existence of light weight
atoms, vertical pressure and slipping of the layers do not make significant
change in optical properties of the structure, therefore, no gap difference
was observed for different thickness of the structure in experiment Hou  _et
al._ (2020).

##  V Conclusion

In summary, we reported the electronic and excitonic optical properties of
$\alpha^{\prime}$-4H borophene using first-principles many-body calculations.
We firstly tracked the mechanism in which the metallic $\alpha^{\prime}$
borophene turns into the semiconducting $\alpha^{\prime}$-4H borophene. It was
found that hydrogenation increases the buckling height and breaks the mirror
symmetry of $\alpha^{\prime}$ borophene, which results in hybridization
between $p_{z}$ and neighboring $p_{x(y)}$ orbitals and opens a wide band gap
in $\alpha^{\prime}$-4H borophene. Then, we calculated the quasiparticle bang
gap to be 1.98, 2.23, and 2.52 eV at the G0W0, GW0, and GW levels of theory,
respectively. In the following, we solved the Bethe-Salpeter equation on top
of the GW. The optical spectrum at GW+BSE shows that the direct optical
excitation is characterized by an optical band gap of $\sim$ 2.40 eV, which is
in excellent agreement with the experiment (2.49 eV). The exciton binding
energy was calculated to be $\sim$ 1.18 eV, which is a signature of a strongly
bound Frenkel exciton, suggesting ultrahigh stability of excitonic states
against thermal dissociation. We also investigated the bilayer
$\alpha^{\prime}$-4H borophenes with different stackings and thicknesses. If
was found that they have similar electronic properties to the monolayer, which
gives them similar potentials. These features together with great air-
stability and mechanical strength suggest $\alpha^{\prime}$-4H borophene as a
desirable material for flexible optoelectronic applications at room
temperature.

## Acknowledgment

We are thankful to the Research Council of the University of Guilan for the
partial support of this research.

## Declaration of Interests

The authors declare that they have no conflict of interest.

## References

  * Mannix  _et al._ (2018) A. J. Mannix, Z. Zhang, N. P. Guisinger, B. I. Yakobson, and M. C. Hersam, Nat. Nanotechnol. 13, 444 (2018).
  * Liu  _et al._ (2020a) L. Liu, Z. Zhang, X. Liu, X. Xuan, B. I. Yakobson, M. C. Hersam, and W. Guo, Nano Lett. 20, 1315 (2020a).
  * Penev  _et al._ (2016) E. S. Penev, A. Kutana, and B. I. Yakobson, Nano Lett. 16, 2522 (2016).
  * Vishkayi and Tagani (2017) S. I. Vishkayi and M. B. Tagani, Phys. Chem. Chem. Phys. 19, 21461 (2017).
  * Vajary  _et al._ (2018) R. A. Vajary, M. B. Tagani, and S. I. Vishkayi, Mod. Phys. Lett. B 32, 1850347 (2018).
  * Vishkayi and Tagani (2018) S. I. Vishkayi and M. B. Tagani, Phys. Chem. Chem. Phys. 20, 10493 (2018).
  * Mannix  _et al._ (2015) A. J. Mannix, X.-F. Zhou, B. Kiraly, J. D. Wood, D. Alducin, B. D. Myers, X. Liu, B. L. Fisher, U. Santiago, and J. R. Guest, Science 350, 1513 (2015).
  * Arabha  _et al._ (2020) S. Arabha, A. H. Akbarzadeh, and A. Rajabpour, Compos. B. Eng. 200, 108260 (2020).
  * Ayodhya and Veerabhadram (2020) D. Ayodhya and G. Veerabhadram, FlatChem 19, 100150 (2020).
  * Hou  _et al._ (2020) C. Hou, G. Tai, J. Hao, L. Sheng, B. Liu, and Z. Wu, Angew. Chem. (2020).
  * Xie  _et al._ (2020) Z. Xie, X. Meng, X. Li, W. Liang, W. Huang, K. Chen, J. Chen, C. Xing, M. Qiu, B. Zhang, G. Nie, N. Xie, X. Yan, and H. Zhang, Research 2020, 2624617 (2020).
  * Entwistle and Marder (2002) C. D. Entwistle and T. B. Marder, Angew. Chem. Int. Ed. 41, 2927 (2002).
  * Kong  _et al._ (2019) L. Kong, L. Liu, L. Chen, Q. Zhong, P. Cheng, H. Li, Z. Zhang, and K. Wu, Nanoscale 11, 15605 (2019).
  * Powar and Pandav (2019) N. Powar and R. Pandav, J. Chem. Rev. 1, 271 (2019).
  * Bhattacharyya  _et al._ (2017) G. Bhattacharyya, A. Mahata, I. Choudhuri, and B. Pathak, J. Phys. D: Appl. Phys. 50, 405103 (2017).
  * Suehara  _et al._ (2010) S. Suehara, T. Aizawa, and T. Sasaki, Phys. Rev. B 81, 085423 (2010).
  * Wu  _et al._ (2019) R. Wu, I. K. Drozdov, S. Eltinge, P. Zahl, S. Ismail-Beigi, I. Božović, and A. Gozar, Nat. Nanotechnol. 14, 44 (2019).
  * Kiraly  _et al._ (2019) B. Kiraly, X. Liu, L. Wang, Z. Zhang, A. J. Mannix, B. L. Fisher, B. I. Yakobson, M. C. Hersam, and N. P. Guisinger, ACS Nano 13, 3816 (2019).
  * Zhong  _et al._ (2017) Q. Zhong, L. Kong, J. Gou, W. Li, S. Sheng, S. Yang, P. Cheng, H. Li, K. Wu, and L. Chen, Phys. Rev. Mater. 1, 021001 (2017).
  * Liu  _et al._ (2013) H. Liu, J. Gao, and J. Zhao, Sci. Rep. 3, 1 (2013).
  * Ryder  _et al._ (2016) C. R. Ryder, J. D. Wood, S. A. Wells, Y. Yang, D. Jariwala, T. J. Marks, G. C. Schatz, and M. C. Hersam, Nat. Chem. 8, 597 (2016).
  * Johns and Hersam (2013) J. E. Johns and M. C. Hersam, Acc. Chem. Res. 46, 77 (2013).
  * Li  _et al._ (2021) Q. Li, V. S. C. Kolluru, M. S. Rahn, E. Schwenker, S. Li, R. G. Hennig, P. Darancet, M. K. Y. Chan, and M. C. Hersam, Science 371, 1143 (2021).
  * Mohebpour  _et al._ (2020) M. A. Mohebpour, S. M. Mozvashi, S. I. Vishkayi, and M. B. Tagani, Sci. Rep. 10, 1 (2020).
  * Mozvashi  _et al._ (2021) S. M. Mozvashi, M. A. Mohebpour, S. I. Vishkayi, and M. B. Tagani, Sci. Rep. 11, 1 (2021).
  * Cudazzo  _et al._ (2010) P. Cudazzo, C. Attaccalite, I. V. Tokatly, and A. Rubio, Phys. Rev. Lett. 104, 226804 (2010).
  * Yang  _et al._ (2009) L. Yang, J. Deslippe, C.-H. Park, M. L. Cohen, and S. G. Louie, Phys. Rev. Lett. 103, 186802 (2009).
  * Shishkin and Kresse (2007) M. Shishkin and G. Kresse, Phys. Rev. B 75, 235102 (2007).
  * Hybertsen and Louie (1986) M. S. Hybertsen and S. G. Louie, Phys. Rev. B 34, 5390 (1986).
  * Liu  _et al._ (2020b) C. Liu, J. Kloppenburg, Y. Yao, X. Ren, H. Appel, Y. Kanai, and V. Blum, J. Chem. Phys. 152, 044105 (2020b).
  * Karlicky and Otyepka (2013) F. Karlicky and M. Otyepka, J. Chem. Theory Comput. 9, 4155 (2013).
  * Shahrokhi and Leonard (2016) M. Shahrokhi and C. Leonard, J. Alloys Compd. 682, 254 (2016).
  * Shishkin and Kresse (2006) M. Shishkin and G. Kresse, Phys. Rev. B 74, 035101 (2006).
  * Blase  _et al._ (2020) X. Blase, I. Duchemin, D. Jacquemin, and P.-F. Loos, J. Phys. Chem. Lett. 11, 7371 (2020).
  * Gajdoš  _et al._ (2006) M. Gajdoš, K. Hummer, G. Kresse, J. Furthmüller, and F. Bechstedt, Phys. Rev. B 73, 045112 (2006).
  * Shahrokhi and Leonard (2017) M. Shahrokhi and C. Leonard, J. Alloys Compd. 693, 1185 (2017).
  * Rohlfing and Louie (1998) M. Rohlfing and S. G. Louie, Phys. Rev. Lett. 80, 3320 (1998).
  * Rohlfing and Louie (2000) M. Rohlfing and S. G. Louie, Phys. Rev. B 62, 4927 (2000).
  * Giannozzi  _et al._ (2009) P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni, and I. Dabo, J. Phys.: Condens. Matter 21, 395502 (2009).
  * Mostofi  _et al._ (2014) A. A. Mostofi, J. R. Yates, G. Pizzi, Y.-S. Lee, I. Souza, D. Vanderbilt, and N. Marzari, Comput. Phys. Commun. 185, 2309 (2014).
  * Ezawa (2017) M. Ezawa, Phys. Rev. B 96, 035425 (2017).
  * Feng  _et al._ (2017) B. Feng, O. Sugino, R.-Y. Liu, J. Zhang, R. Yukawa, M. Kawamura, T. Iimori, H. Kim, Y. Hasegawa, and H. Li, Phys. Rev. Lett. 118, 096401 (2017).
  * Radha and Lambrecht (2020) S. K. Radha and W. R. L. Lambrecht, Phys. Rev. B 101, 235111 (2020).
  * Umeno  _et al._ (2019) Y. Umeno, M. Sato, M. Sato, and H. Shima, Phys. Rev. B 100, 155116 (2019).
  * Nijamudheen  _et al._ (2015) A. Nijamudheen, R. Bhattacharjee, S. Choudhury, and A. Datta, J. Phys. Chem. C 119, 3802 (2015).
  * Qiu  _et al._ (2013) D. Y. Qiu, H. Felipe, and S. G. Louie, Phys. Rev. Lett. 111, 216805 (2013).
  * Kolos and Karlický (2019) M. Kolos and F. Karlický, Phys. Chem. Chem. Phys. 21, 3999 (2019).
  * Xu  _et al._ (2017) L. Xu, M. Yang, S. J. Wang, and Y. P. Feng, Phys. Rev. B 95, 235434 (2017).
  * Tran  _et al._ (2014) V. Tran, R. Soklaski, Y. Liang, and L. Yang, Phys. Rev. B 89, 235319 (2014).
  * Yang  _et al._ (2007) L. Yang, C. D. Spataru, S. G. Louie, and M.-Y. Chou, Phys. Rev. B 75, 201304 (2007).
  * Choi  _et al._ (2015) J.-H. Choi, P. Cui, H. Lan, and Z. Zhang, Phys. Rev. Lett. 115, 066403 (2015).
  * Zhang  _et al._ (2009) Y. Zhang, T.-T. Tang, C. Girit, Z. Hao, M. C. Martin, A. Zettl, M. F. Crommie, Y. R. Shen, and F. Wang, Nature 459, 820 (2009).
  * Wang  _et al._ (2015) G. Wang, R. Pandey, and S. P. Karna, ACS Appl. Mater. Interfaces 7, 11490 (2015).