# Tunable nonlinear excitonic optical response in biased bilayer graphene

M. F. C. Martins Quintela1,2,3<EMAIL_ADDRESS>N. M. R. Peres1,2,4 T.
Garm Pedersen3 1Department of Physics and Physics Center of Minho and Porto
Universities (CF–UM–UP), Campus of Gualtar, 4710-057, Braga, Portugal
2International Iberian Nanotechnology Laboratory (INL), Av. Mestre José Veiga,
4715-330, Braga, Portugal  3Department of Materials and Production, Aalborg
University, 9220 Aalborg Øst, Denmark  4POLIMA–Center for Polariton–driven
Light–Matter Interactions, University of Southern Denmark, Campusvej 55,
DK-5230 Odense M, Denmark

###### Abstract

Biased bilayer graphene (BBG) is an important system for studies of excitonic
effects in graphene–based systems, with its easily tunable bandgap. This
bandgap is governed by an external gate voltage, allowing one to tune the
optical response of the system. In this paper, we study the excitonic linear
and nonlinear optical response of Bernal stacked BBG as a function of the gate
voltage, both for in–plane (IP) and out–of–plane (OOP) directions. Based on a
semi-analytical model of the electronic structure of BBG describing the
influence of gate voltage on excitonic binding energies, we focus our
discussion on both the IP and OOP excitonic response. Both linear and second
harmonic generation (SHG) nonlinear responses are shown to be very sensitive
to the gate voltage, as both the interband momentum matrix elements and the
bandgap of the system will vary greatly with bias potential.

##  I Introduction

Since the first mechanical isolation of graphene[1], many different
two–dimensional (2D) materials have been studied[2], such as hexagonal Boron
Nitride (hBN)[3] and transition metal dichalcogenides (TMDs)[4]. This interest
then recently transitioned into layered materials with broken vertical
symmetry, which have been considered both from theoretical and experimental
perspectives. These materials include, but are not limited to, buckled
monolayers [5, 6, 7, 8], Janus materials [9, 10, 11, 12, 13], hetero- and
biased homobilayers [14, 15, 16, 17]. In stark contrast to its monolayer
counterpart, BBG has a tunable gap which can be as large as
$150\,\mathrm{meV}$[18]. The magnitude of this bandgap can be controlled via
an external electric field and dielectric environment, and allows the
formation of tunable excitons, which have already been both measured
experimentally[19] and described theoretically[14, 20, 21, 22, 23].

As BBG presents an intrinsic OOP asymmetry due to the interlayer bias
potential, one can directly study the effects of this asymmetry on the
excitonic response. Tuning the interlayer bias potential, one can alter both
the electronic structure of BBG and the interlayer asymmetry, probing the
effects on both IP and OOP optical response[24]. The probing of the OOP
response, namely the OOP non–linearities, leads to additional degrees of
freedom useful for vertical photonics structures[25, 26], allowing for novel
approaches in the design of ultrafast optical devices[12], such as non–linear
holograms[27], broadband ultrafast frequency converters[28, 29], miniaturized
logic gates[30, 31], among others. Furthermore, the presence of various
non–zero tensor components which include non–diagonal components also provides
a greater freedom when experimentally probing the optical response of BBG.

Recent theoretical works[32] have reviewed the effects of gate voltage on the
free–carrier injection rate and shift current in AA- and AB- (Bernal) stacked
bilayer graphene, with many other theoretical works focusing on
single–particle properties of bilayer graphene structures[33, 34, 35, 36, 37,
38], especially when an interlayer twist angle is introduced[39, 40, 41, 42,
43]. However, and as in many layered materials, the optical response of BBG is
dominated by excitons[44] and their effects are, therefore, fundamental in the
study of the optical response of this material. In this paper, we focus on the
effects of gate voltage on the excitonic optical response of AB-stacked
bilayer graphene, considering both linear and SHG nonlinear response.

This paper is organized as follows: we begin, in Section II, by discussing the
electronic structure and its dependence on the external gate voltage. This is
then followed by the computation of the momentum and Berry connection matrix
elements, both for IP and OOP directions. Knowing the momentum matrix
elements, we compute the free–carrier linear response for both the full
Hamiltonian and a reduced two–band model. This allows us to identify the
frequency ranges where two bands are enough to describe the optical response
of the system as a function of the gate voltage. Afterwards, in Section III,
we proceed to the computation of the excitonic states through solving the
Bethe–Salpeter equation. This involves an approximation of the effective
screening length in the electron–hole electrostatic potential, followed by a
discussion on the distinction between intralayer and interlayer phenomena.
Finally, in Section IV, we compute the linear and SHG nonlinear excitonic
optical response of BBG. After obtaining the IP and OOP optical selection
rules for both linear and nonlinear response, we discuss the effects of gate
voltage, as well as the relative amplitudes of the different elements of the
nonlinear conductivity tensor.

##  II Single Particle Biased Bilayer Graphene Hamiltonian

In the basis
$\left\\{\left|1,b\right\rangle,\left|2,b\right\rangle,\left|1,t\right\rangle,\left|2,t\right\rangle\right\\}$,
where $b/t$ denotes the bottom/top layer and $1/2$ denotes the different
sublattices in each layer, the Hamiltonian of Bernal stacked BBG[45] in the
Dirac approximation around the Dirac points $\left(K_{x},\tau K_{y}\right)$,
with $\tau=\pm 1$ the valley index, is given by

$\mathcal{H}=\left(\begin{array}[]{cccc}V&-\gamma_{0}f&0&-i\gamma_{3}f^{*}\\\
-\gamma_{0}f^{*}&V&\gamma_{z}&0\\\ 0&\gamma_{z}&-V&-\gamma_{0}f\\\
i\gamma_{3}f&0&-\gamma_{0}f^{*}&-V\end{array}\right),$ (1)

where $V$ is the bias potential applied to the bilayer system,
$\gamma_{0}=3.12\,\mathrm{eV}$ is the intralayer nearest–neighbour hopping
parameter, and the interlayer parameters $\gamma_{z}=0.377\,\mathrm{eV}$
(interlayer direct vertical coupling) and $\gamma_{3}=0.377\,\mathrm{eV}$
(skew coupling parameter). Restricting our analysis to the Dirac cones
$K,K^{\prime}$, we consider $\mathbf{k}$ measured from the Dirac points and
write

$\displaystyle f\left(\mathbf{k}\right)$
$\displaystyle\approx\frac{\sqrt{3}}{2}a\left(k_{x}-i\tau k_{y}\right)$

with $\hbar v_{F}=\frac{\sqrt{3}}{2}a\gamma_{0}$ the Fermi velocity, the polar
angle of $\mathbf{k}$ defined as
$\theta=\tan^{-1}\left(\frac{k_{y}}{k_{x}}\right)$, and $a=2.46\,\text{\AA}$
the lattice parameter. The interlayer bias potential $2V$ corresponds to
placing the system in an electric field $\frac{2V}{ed}$, with $d$ the
interlayer separation. In this paper we will mainly consider
$V=55\,\mathrm{meV}$ as it is at this bias that the dominant peak in the IP
linear response matches with recent experimental results[19].

Figure 1: (Left) Band structure of BBG with a large bias
$(V=200\,\mathrm{meV})$ near the Dirac points. Coloured dashed lines
correspond to the band structure under zero external bias, vertical arrows
represent the three direct bandgaps. (Right) Evolution of the three distinct
direct bandgaps in BBG with increasing external bias near the Dirac points.
Vertical dashed lines represent the four different values of $V$ applied in
numerical studies.

The band structure under the Dirac approximation is given by

$\displaystyle E_{\lambda_{\eta}}$
$\displaystyle=\lambda\sqrt{\hbar^{2}v_{F}^{2}k^{2}+V^{2}+\frac{1}{2}\gamma_{z}^{2}+\left(2\eta-3\right)\Lambda}$
$\displaystyle\Lambda$
$\displaystyle=\sqrt{\frac{1}{4}\gamma_{z}^{4}+\hbar^{2}v_{F}^{2}k^{2}\left(4V^{2}+\gamma_{z}^{2}\right)},$
(2)

where $\lambda=\pm 1$ represents the conduction and valence bands
($\lambda=+1$ will be conduction (c) bands, $\lambda=-1$ valence (v) bands),
and $\eta=1,2$ identifies the band closest ($\eta=1$) and furthest away
($\eta=2$) from the bandgap. Direct substitution of $\eta=1,2$ into Eq. (2)
leads to a factor of $\pm 1$, as expected for the band structure of this
system. This band structure is presented in Fig. 1. Focusing on direct
bandgaps, we write the generic bandgap as

$E_{g}^{(\eta,\eta^{\prime})}=E_{c_{\eta}}-E_{v_{\eta^{\prime}}}.$ (3)

Due to the nature of the BBG band structure, the minimum bandgap between the
two $\eta=1$ bands will not occur at $k=0$ (see left panel of Fig. 1), instead
occurring at

$k=\sqrt{2}\frac{V}{\hbar
v_{F}}\sqrt{\frac{2V^{2}+\gamma_{z}^{2}}{4V^{2}+\gamma_{z}^{2}}}$ (4)

measured from the Dirac points. For instance,

$\displaystyle E_{g}^{(11)}$
$\displaystyle=\frac{2V\gamma_{z}}{\sqrt{4V^{2}+\gamma_{z}^{2}}}\approx
106\,\mathrm{meV}$ (5)

for the considered bias potential $V=55\,\mathrm{meV}$. The two larger direct
bandgaps, namely $E_{g}^{(12)}$ and $E_{g}^{(22)}$, will occur at $k=0$ and
are significantly larger. For $V=55\,\mathrm{meV}$, these take the values

$\displaystyle E_{g}^{(12)}$
$\displaystyle=V+\sqrt{V^{2}+\gamma_{z}^{2}}\approx 436\,\mathrm{meV},$
$\displaystyle E_{g}^{(22)}$
$\displaystyle=2\sqrt{V^{2}+\gamma_{z}^{2}}\approx 762\,\mathrm{meV},$ (6)

and $E_{g}^{(21)}=E_{g}^{(12)}.$ The three distinct bandgaps are plotted in
the right panel of Fig. 1 for external bias between $0$ and
$110\,\mathrm{meV}$. The large difference between the three bandgaps suggests
that, as we will see ahead, the two $\eta=1$ bands will dominate the
low–energy response of the system.

###  II.1 In–Plane Matrix Elements

For simplicity, we will discard the contributions from skew coupling to the
eigenvectors, considering it instead as a perturbation when computing dipole
matrix elements. In a compact notation, we define the normalized eigenvectors
for future computations as

$\left|\lambda_{\eta}\left(\mathbf{k}\right)\right\rangle=\left[\begin{array}[]{c}e^{-2i\tau\theta}\alpha_{1,\lambda_{\eta}}\left(k\right)\\\
e^{-i\tau\theta}\alpha_{2,\lambda_{\eta}}\left(k\right)\\\
e^{-i\tau\theta}\alpha_{3,\lambda_{\eta}}\left(k\right)\\\
\alpha_{4,\lambda_{\eta}}\left(k\right)\end{array}\right].$ (7)

This gauge choice will carry a pseudo–spin[18, 21, 46, 47, 48, 49] factor of
$m_{s}=-2$ in the $\tau=1$ valley. For compactness in the expressions for the
momentum matrix elements and Berry connections, the $k$ dependence in the
eigenvector components will be suppressed.

To compute the IP interband momentum matrix elements, chosen without loss of
generality in the $x$ direction, we must analyze

$P_{v_{\eta}c_{\eta^{\prime}}\mathbf{k}}^{x}=\frac{m_{0}}{\hbar}\left<v_{\eta}\middle|\frac{\partial\mathcal{H}}{\partial
k_{x}}\middle|c_{\eta^{\prime}}\right>,$ (8)

where $m_{0}$ is the free electron mass. The operator
$\frac{\partial\mathcal{H}}{\partial k_{x}}$ can be readily computed and is
given by

$\displaystyle\frac{\partial\mathcal{H}}{\partial k_{x}}$
$\displaystyle=\left(\begin{array}[]{cccc}0&-\hbar
v_{F}&0&-i\frac{\gamma_{3}}{\gamma_{0}}\hbar v_{F}\\\ -\hbar v_{F}&0&0&0\\\
0&0&0&-\hbar v_{F}\\\ i\frac{\gamma_{3}}{\gamma_{0}}\hbar v_{F}&0&-\hbar
v_{F}&0\end{array}\right).$

Yet, obtaining the matrix elements taking into account skew coupling due to
finite $\gamma_{3}$ in the eigenvectors becomes extremely unwieldy. As a
consequence, we consider an approximation where $\gamma_{3}$ is taken as a
small perturbation to the Dirac BBG system and, therefore, its only effect
will be introducing additional terms in the momentum matrix elements
proportional to $\frac{\gamma_{3}}{\gamma_{0}}\hbar v_{F}$, as presented in
$\frac{\partial\mathcal{H}}{\partial k_{x}}$. This follows the approach
previously outlined in a similar study of the linear response of BBG[21],
where experimental results[19] where accurately predicted under this same
perturbative approach. This also means that skew coupling will play no part in
the calculation of Berry connections, as it is only dependent on the
eigenvectors defined in Eq. (7).

Focusing only on the $\eta=1$ bands, the exciton momentum matrix element will
be given explicitly as

$\displaystyle
P_{v_{\eta}c_{\eta^{\prime}}\mathbf{k}}^{x}=m_{0}v_{F}\left[e^{i\tau\theta}\left(\alpha_{1,v_{\eta}}^{*}\alpha_{2,c_{\eta^{\prime}}}+\alpha_{3,v_{\eta}}^{*}\alpha_{4,c_{\eta^{\prime}}}\right)\right.$
$\displaystyle\,+e^{-i\tau\theta}\left(\alpha_{2,v_{\eta}}^{*}\alpha_{1,c_{\eta^{\prime}}}+\alpha_{4,v_{\eta}}^{*}\alpha_{3,c_{\eta^{\prime}}}\right)$
(9)
$\displaystyle\,\left.-i\frac{\gamma_{3}}{\gamma_{0}}\left(e^{2i\tau\theta}\alpha_{1,v_{\eta}}^{*}\alpha_{4,c_{\eta^{\prime}}}-e^{-2i\tau\theta}\alpha_{4,v_{\eta}}^{*}\alpha_{1,c_{\eta^{\prime}}}\right)\right].$

In this form, it is also clear from the phase factors that, upon
multiplication by the conjugate and angular integration, the lowest
contribution from skew coupling to the linear response will be quadratic in
$\left|\gamma_{3}/\gamma_{0}\right|\approx 0.1$. As
$\left|\gamma_{3}/\gamma_{0}\right|\ll 1$, we will keep only the dominant term
in each expansion. The validity of this perturbative approach has also been
checked against a full tight–binding numerical calculation.

Comparing both the single–particle Hamiltonian as in Eq. (1) and the momentum
matrix element as given in Eq. (9) against those used in Ref. [21], one can
see that Ref. [21] is missing the $i$ factor in the skew coupling contribution
for both the Hamiltonian and the momentum matrix element. Although the lack of
this term is of no consequence for the linear response discussed in Ref. [21]
due to the proportionality to $\left|\gamma_{3}/\gamma_{0}\right|^{2}$, this
factor is of importance for the IP nonlinear response as it will be clear when
the excitonic matrix elements are discussed.

A different contribution, namely trigonal warping, consists of considering the
quadratic term in the series expansion of $f$, leading to[50]

$f\left(\mathbf{k}\right)\approx\frac{\sqrt{3}a}{2}\left[\left(k_{x}-i\tau
k_{y}\right)+i\zeta_{\mathrm{TW}}a\left(k_{x}-i\tau k_{y}\right)^{2}\right],$
(10)

where $\zeta_{\mathrm{TW}}=\frac{\sqrt{3}}{12}$ is a fixed numerical factor
resulting from the series expansion of $f$. Careful analysis of the effects of
this quadratic term on the momentum matrix elements via Eq. (8) reveals that
the phase factor associated will be $e^{\pm 2i\tau\theta}$, identical to the
contribution from skew coupling in Eq. (9). Furthermore, upon computation of
the excitonic matrix elements, this quadratic correction proves to be
negligible when compared against that originating from skew coupling and,
therefore, it will be ignored.

In the $k_{x}$ direction, the matrix elements of the Berry connection between
bands $n_{\eta}$ and $m_{\eta^{\prime}}$ read

$\displaystyle\Omega_{n_{\eta}m_{\eta\prime}\mathbf{k}}^{x}$
$\displaystyle=i\left<n_{\eta}\middle|\frac{\partial}{\partial
k_{x}}\middle|m_{\eta\prime}\right>.$ (11)

The exact form of the Berry connection will be needed when computing
generalized derivatives[51] and
$\Omega_{c_{\eta}c_{\eta^{\prime}}\mathbf{k}}^{x}$ will be given in terms of
the eigenvector components as

$\displaystyle\Omega_{c_{\eta}c_{\eta^{\prime}}\mathbf{k}}^{x}$
$\displaystyle=i\cos\theta\left(\alpha_{1,c_{\eta}}^{*}\frac{\partial\alpha_{1,c_{\eta^{\prime}}}}{\partial
k}+\alpha_{2,c_{\eta}}^{*}\frac{\partial\alpha_{3,c_{\eta^{\prime}}}}{\partial
k}+\alpha_{3,c_{\eta}}^{*}\frac{\partial\alpha_{3,c_{\eta^{\prime}}}}{\partial
k}+\alpha_{4,c_{\eta}}^{*}\frac{\partial\alpha_{4,c_{\eta^{\prime}}}}{\partial
k}\right)$
$\displaystyle\quad-\frac{\tau}{k}\sin\theta\left(2\alpha_{1,c_{\eta}}^{*}\alpha_{1,c_{\eta^{\prime}}}+\alpha_{2,c_{\eta}}^{*}\alpha_{1,c_{\eta^{\prime}}}+\alpha_{3,c_{\eta}}^{*}\alpha_{3,c_{\eta^{\prime}}}\right),$
(12)

and $\Omega_{v_{\eta}v_{\eta^{\prime}}\mathbf{k}}^{x}$ will be defined
analogously for the valence bands.

###  II.2 Out–of–Plane Matrix Elements

Before computing the single particle linear optical response, we will discuss
the OOP momentum and Berry connection. To this end, we begin by defining the
matrix elements of $z$ from the extension of $\frac{d}{2}\sigma_{z}$[52] to a
bilayer system, with $\sigma_{z}$ the diagonal Pauli matrix, as

$\displaystyle z_{n_{\eta}m_{\eta\prime}\mathbf{k}}$
$\displaystyle=\frac{1}{2}\left\langle
n_{\eta}\left|\mathrm{diag}\left[d,d,-d,-d\right]\right|m_{\eta^{\prime}}\right\rangle.$
(13)

The interlayer separation $d$ has been previously discussed in the
literature[53, 54, 55, 56], with reported values for $d\approx
3.46\,\text{\AA}$. The momentum matrix element then reads

$\displaystyle P_{v_{\eta}c_{\eta^{\prime}}\mathbf{k}}^{z}$
$\displaystyle=\frac{m_{0}}{i\hbar}E_{c_{\eta^{\prime}}v_{\eta}\mathbf{k}}z_{v_{\eta}c_{\eta^{\prime}}\mathbf{k}}.$
(14)

Explicitly computing $z_{v_{\eta}c_{\eta^{\prime}}\mathbf{k}}$, we obtain

$\displaystyle P_{v_{\eta}c_{\eta^{\prime}}\mathbf{k}}^{z}$
$\displaystyle=\frac{m_{0}d}{2i\hbar}E_{c_{\eta^{\prime}}v_{\eta}\mathbf{k}}\left[\alpha_{1,v_{\eta}}^{*}\alpha_{1,c_{\eta^{\prime}}}+\alpha_{2,v_{\eta}}^{*}\alpha_{2,c_{\eta^{\prime}}}\right.$
$\displaystyle\quad\left.-\alpha_{3,v_{\eta}}^{*}\alpha_{3,c_{\eta^{\prime}}}-\alpha_{4,v_{\eta}}^{*}\alpha_{4,c_{\eta^{\prime}}}\right].$
(15)

For the nonlinear response, we must also compute the matrix elements
$z_{c_{\eta}c_{\eta^{\prime}}\mathbf{k}}$,
$z_{v_{\eta}v_{\eta^{\prime}}\mathbf{k}}$. Similarly to
$z_{v_{\eta}c_{\eta^{\prime}}\mathbf{k}}$, we write these in terms of the
radial components of the eigenvectors as

$\displaystyle z_{c_{\eta}c_{\eta^{\prime}}\mathbf{k}}$
$\displaystyle=\frac{d}{2}\left[\alpha_{1,c_{\eta}}^{*}\alpha_{1,c_{\eta^{\prime}}}+\alpha_{2,c_{\eta}}^{*}\alpha_{2,c_{\eta^{\prime}}}\right.$
$\displaystyle\quad\left.-\alpha_{3,c_{\eta}}^{*}\alpha_{3,c_{\eta^{\prime}}}-\alpha_{4,c_{\eta}}^{*}\alpha_{4,c_{\eta^{\prime}}}\right],$

and analogous for $z_{v_{\eta}v_{\eta^{\prime}}}$. Hence, we obtain the Berry
connection as[51]

$\displaystyle\Omega_{c_{\eta}c_{\eta^{\prime}}\mathbf{k}}^{z}-\Omega_{v_{\eta}v_{\eta^{\prime}}\mathbf{k}}^{z}=z_{c_{\eta}c_{\eta^{\prime}}\mathbf{k}}-z_{v_{\eta}v_{\eta^{\prime}}\mathbf{k}}.$
(16)

###  II.3 Single Particle Linear Optical Response

In a clean semiconductor at zero temperature, the diagonal optical response
reads [57, 51, 58, 59, 60]

$\sigma_{\beta\beta}\left(\omega\right)=\frac{e^{2}\hbar^{2}\omega}{i\pi^{2}m_{0}^{2}}\sum_{c,v}\int\frac{\left|P_{cv\mathbf{k}}^{\beta}\right|^{2}}{E_{cv\mathbf{k}}\left(E_{cv\mathbf{k}}^{2}-\hbar^{2}\omega^{2}\right)}d^{2}\mathbf{k},$
(17)

where $\sum_{c,v}$ sums over all conduction and valence bands. Within the
two–band approximation, the sum over valence and conduction bands is dropped.
With the expansion near the Dirac cones, integration will now be over the
infinite Dirac cones and a sum over valleys must also be made.

Figure 2: Comparison between full four–band calculation and calculation with
only the two lowest bands of both IP (top) and OOP (bottom) free–carrier
linear optical response for biased bilayer graphene with a bias potential
$V=55\,\mathrm{meV}$ and broadening $\hbar\Gamma=5\,\mathrm{meV}$. Real part
represented in blue, denoted $\Re_{n}$, and imaginary part in orange, denoted
$\Im_{n}$, for the $n-$band calculation. Vertical dashed lines represent the
various bandgaps of the system, from left to right $E_{g}^{(11)}$,
$E_{g}^{(12)}$, $E_{g}^{(22)}$.

Inspecting Eq. (9), one sees that the angular integration is trivial. Although
the radial integration proves much more complicated, one can easily perform it
numerically and, introducing a broadening parameter via the transformation
$\hbar\omega\rightarrow\hbar\omega+i\hbar\Gamma$ with
$\hbar\Gamma=5\,\mathrm{meV}$, obtain both the free–carrier IP (top) and OOP
(bottom) linear optical response as depicted in Fig. 2, with
$\sigma_{0}=e^{2}/4\hbar$.

While the free–carrier response can be easily computed for the full four–band
system, the same is not true when one wishes to obtain the excitonic response.
The four possible combinations of valence and conduction bands would lead to a
drastic increase in the computational complexity of both solving the
Bethe–Salpeter equation and of obtaining the excitonic momentum matrix
elements[49, 21, 61]. Hence, it becomes crucial to reduce the system to a
two–band problem when one wants to consider its excitonic properties. As seen
in Fig. 2, and as expected from the large separation of the higher energy
bands, we observe minimal differences in the real part of the linear response
between the four– and two–band responses if $\hbar\omega\lesssim
400\,\mathrm{meV}$, consistent with previous results [21, 62]. Hence, we will
consider only the two $\eta=1$ bands from this point onward. This will,
however, introduce issues when the imaginary part of $\sigma$ is relevant
(i.e., when considering the magnitude of the various components of the
conductivity tensor). Additionally, as only the two $\eta=1$ bands have an
effect on the real part of the conductivity up to $\hbar\omega\approx
E_{g}^{(12)}$ (for $V=55\,\mathrm{meV}$, up to $\hbar\omega\approx
400\,\mathrm{meV}$), one can expect minimal effects on the real part of the
SHG nonlinear conductivity up to $\hbar\omega\approx E_{g}^{(12)}/2$. The
imaginary part will, however, be significantly different from the four–band
calculation much earlier due to both the longer range of the peaks in the
imaginary part and the presence of resonances at $2\hbar\omega$.

##  III Bethe–Salpeter Equation

As discussed in the preceding section, we will focus on the dominant response
from the two bands closest to the bandgap $c_{1},\,v_{1}$, which we denote
$c,v$. Hence, we are only interested in the eigenvalues $E_{c/v}$ and their
respective eigenvectors, $\left|c/v\right\rangle$. Additionally, we will
ignore contributions to eigenvectors from the skew coupling parameter, as this
will allow us to transform the Bethe–Salpeter equation (BSE) from a 2D
integral equation into a 1D integral equation.

Figure 3: Evolution of the effective screening length $r_{0}$ as a function of
external bias $V$.

Before discussing the excitonic conductivity, we must first compute the
excitonic states. To compute the excitonic wave functions and binding
energies, we will solve the BSE [63, 60, 50, 48, 64], given in momentum space
by

$\begin{aligned}
E_{n}\psi_{cv\mathbf{k}}^{\left(n\right)}&=E_{cv\mathbf{k}}\psi_{cv\mathbf{k}}^{\left(n\right)}+\\\
&+\sum_{\mathbf{q}}V\left(\left|\mathbf{k}-\mathbf{q}\right|\right)\left\langle
c_{\mathbf{k}}\middle|c_{\mathbf{q}}\right\rangle\left\langle
v_{\mathbf{q}}\middle|v_{\mathbf{k}}\right\rangle\psi_{cv\mathbf{q}}^{\left(n\right)}\end{aligned},$
(18)

where $E_{n}$ is the exciton energy of state $n$, $V\left(k\right)$ is the
attractive electrostatic potential coupling the electron and the hole, and
$\psi_{cv\mathbf{k}}^{\left(n\right)}$ is the wave function of the exciton.
For our system, we consider $V\left(k\right)$ to be the Rytova–Keldysh
potential[65, 66], given in momentum space by

$V\left(k\right)=-2\pi\hbar c\alpha\frac{1}{k\left(\epsilon+r_{0}k\right)},$
(19)

with $\alpha$ the fine–structure constant, $\epsilon$ the mean dielectric
constant of the media surrounding the bilayer, and $r_{0}$ as an in–plane
screening length[67] related to the polarizability of the material and usually
obtained from DFT calculations[68]. The effective screening length can be
taken as a first approximation from the dipolar transition amplitudes and is
given by[67]

$\displaystyle r_{0}$ $\displaystyle=\frac{\hbar^{3}c\alpha}{\pi
m_{0}^{2}}\int\frac{\left|\left\langle
c_{\mathbf{k}}\left|P_{cv\mathbf{k}}^{x}\right|v_{\mathbf{k}}\right\rangle\right|^{2}}{E_{cv\mathbf{k}}^{3}}d^{2}\mathbf{k}.$
(20)

This effective screening length is very sensitive to the external bias,
falling quickly as $V$ increases, as can be seen in Fig. 3. For the considered
bias potential, $r_{0}\approx 103\,\text{\AA}$. Similarly to previous
studies[21], the dielectric constant is set at $\epsilon=6.9$ which
corresponds to the the case of BBG encapsulated in hBN at the zero–frequency
limit[69].

We then consider the excitonic wave function to have a well–defined angular
momentum $\ell_{n}$, writing it as

$\displaystyle\psi_{cv\mathbf{k}}^{(n)}=f_{cvk}^{(n)}e^{i\ell_{n}\theta_{k}}.$
(21)

As the bands structure is isotropic, we replace $E_{cv\mathbf{k}}\rightarrow
E_{cvk}$, allowing us to rewrite the Bethe–Salpeter equation converting the
sum into an integral as

$\displaystyle E_{n}f_{cvk}^{(n)}$
$\displaystyle=E_{cvk}f_{cvk}^{(n)}+\frac{1}{4\pi^{2}}\sum_{\ell=-2}^{2}\int_{0}^{+\infty}\int_{0}^{2\pi}V\left(\left|\mathbf{k}-\mathbf{q}\right|\right)\mathcal{A}_{\ell}\left(k,q\right)e^{i\ell\tau\varphi}f_{cvq}^{(n)}e^{i\ell_{n}\varphi}d\varphi\,qdq,$
(22)

where $\varphi=\theta_{q}-\theta_{k}$. The radial component of the form factor
can then be written by analyzing the expansion of $\left\langle
c_{\mathbf{k}}\middle|c_{\mathbf{q}}\right\rangle\left\langle
v_{\mathbf{q}}\middle|v_{\mathbf{k}}\right\rangle$ as in terms of the $\ell$
factors in the complex exponential of the form $e^{i\ell\tau\varphi}$ in Eq.
(22). This then fixes the ranges on the sum over $\ell$ present in Eq. (22).
For a monolayer system, the form factor reads, in a somewhat abusive notation,

$\displaystyle\mathcal{A}_{\ell}\left(k,q\right)e^{i\ell\tau\varphi}$
$\displaystyle=\left\langle
c_{\mathbf{k}}\middle|c_{\mathbf{q}}\right\rangle\left\langle
v_{\mathbf{q}}\middle|v_{\mathbf{k}}\right\rangle,$ (23)

identical to what is present in Eq. (18). However, when one wishes to consider
a bilayer system where there is a distinction between interlayer and
intralayer phenomena, one must compute the form factor more carefully.
Recalling the eigenvectors from Eq. (7), separated into top- and bottom–layer
components as

$\displaystyle\left|\lambda\right\rangle$
$\displaystyle=\left[\begin{array}[]{c}\left|\lambda^{t\;}\right\rangle\\\
\left|\lambda^{b}\right\rangle\end{array}\right],$ (26)

we explicitly introduce a distinction between intralayer and interlayer
interactions by writing the form factor as

$\displaystyle\mathcal{A}_{\ell}\left(k,q\right)e^{i\ell\tau\varphi}=\left\langle
c_{\mathbf{k}}^{t}\middle|c_{\mathbf{q}}^{t}\right\rangle\left\langle
v_{\mathbf{q}}^{t}\middle|v_{\mathbf{k}}^{t}\right\rangle+\left\langle
c_{\mathbf{k}}^{b}\middle|c_{\mathbf{q}}^{b}\right\rangle\left\langle
v_{\mathbf{q}}^{b}\middle|v_{\mathbf{k}}^{b}\right\rangle$
$\displaystyle+e^{-d\left|\mathbf{k}-\mathbf{q}\right|}\left[\left\langle
c_{\mathbf{k}}^{t}\middle|c_{\mathbf{q}}^{t}\right\rangle\left\langle
v_{\mathbf{q}}^{b}\middle|v_{\mathbf{k}}^{b}\right\rangle+\left\langle
c_{\mathbf{k}}^{b}\middle|c_{\mathbf{q}}^{b}\right\rangle\left\langle
v_{\mathbf{q}}^{t}\middle|v_{\mathbf{k}}^{t}\right\rangle\right],$ (27)

where $d\approx 3.46\,\text{\AA}$ is the interlayer separation between the two
graphene layers[53, 54, 55, 56]. The extra
$e^{-d\left|\mathbf{k}-\mathbf{q}\right|}$ factor present in Eq. (27) in fact
originates from the Rytova–Keldysh potential, as written in Eq. (19),
corresponding to the charges being separated in the $z$-direction. It is,
however, simpler to write this factor in the form factor as to explicitly
couple only eigenvector components from opposing layers. As required, when the
interlayer separation vanishes there will no longer be a vertical separation
of the charges and the form factor exactly matches Eq. (23).

Finally, Eq. (22) is solved numerically via a simple numerical quadrature in a
tangent grid $q=\tan\left(x\frac{\pi}{2}\right)$, with $1500$ points
$x\in\left[0,1\right]$, following the procedure already outlined several times
in the literature, namely in Refs. [70, 71, 21, 62]. Additionally, due to the
introduction of the $e^{-d\left|\mathbf{k}-\mathbf{q}\right|}$ in Eq. (27),
the angular integral in Eq. (22) is done numerically in an uniform grid with
$1500$ points $\varphi\in\left[0,2\pi\right]$.

When discussing excitonic states, we will use a similar nomenclature as in the
2D Hydrogen atom[72] to distinguish the different angular momentum states
(_i.e._ , $s$, $p_{\pm}$, $d_{\pm}$ states). In Fig. 4, we plot the binding
energy of the first four states of the $s$, $p_{-}$ and $d_{-}$ series. As the
bias potential $V$ increases, the energy of the first $s$ and $p_{-}$ states
approach each other, eventually crossing in the same way as presented in
current _ab-initio_ computations[22]. In the model considered in this paper,
this crossing occurs for a bias potential $V_{s-p_{-}}\approx
180\,\mathrm{meV}$ and will lead to a swap of the first two peaks of the
linear conductivity, as seen in Ref. [22]. The second crossing, involving $s$
and $d_{-}$ states, will occur at $V_{s-d_{-}}\approx 245\,\mathrm{meV}$, and
will lead to the dominant resonance for the linear OOP response being at a
lower energy than the $1s$ resonance of the linear IP response.

Figure 4: Evolution of the binding energies of the first four $s$, $p_{-}$,
and $d_{-}$ -series excitons as a function of external bias $V$.

##  IV Excitonic Optical Response

Having discussed the method for solving the BSE in a BBG system, we now
proceed to the study of the excitonic optical conductivity, writing the linear
conductivity as[60, 50]

$\frac{\sigma_{\alpha\beta}\left(\omega\right)}{\sigma_{0}}=\frac{-i\hbar^{2}}{2\pi^{3}m_{0}^{2}}\sum_{n}\left[\frac{E_{n}X_{0n}^{\alpha}X_{n0}^{\beta}}{E_{n}-\hbar\omega}-\left(\omega\rightarrow-\omega\right)^{*}\right],$
(28)

where $\sigma_{0}=\frac{e^{2}}{4\hbar}$ is the conductivity of monolayer
graphene. The excitonic matrix elements are defined as[60, 50, 52]

$X_{0n}^{\alpha}=i\int\,\psi_{cv\mathbf{k}}^{\left(n\right)}\frac{P_{vc\mathbf{k}}^{\alpha}}{E_{cv\mathbf{k}}}d^{2}\mathbf{k}.$
(29)

For the nonlinear conductivity, we define
$\sigma_{2}=\frac{e^{3}a}{4E_{g}\hbar}$ and write the SHG
$\left(\omega_{p}=\omega_{q}\right)$ nonlinear conductivity[60, 50] as

$\frac{\sigma_{\alpha\beta\gamma}^{\mathrm{SHG}}\left(\omega\right)}{\sigma_{2}}=\frac{-iE_{g}\hbar^{2}}{2a\pi^{3}m_{0}^{2}}\sum_{n,m}\left[\frac{E_{n}X_{0n}^{\alpha}Q_{nm}^{\beta}X_{m0}^{\gamma}}{\left(E_{n}-2\hbar\omega\right)\left(E_{m}-\hbar\omega\right)}-\frac{E_{n}X_{n0}^{\alpha}Q_{mn}^{\beta}X_{0m}^{\gamma}}{\left(E_{n}+2\hbar\omega\right)\left(E_{m}+\hbar\omega\right)}-\frac{\left(E_{n}-E_{m}\right)X_{0n}^{\alpha}Q_{nm}^{\beta}X_{m0}^{\gamma}}{\left(E_{n}+\hbar\omega\right)\left(E_{m}-\hbar\omega\right)}\right],$
(30)

where $E_{g}$ is the bandgap of the material. As we are only considering
effects from the two $\eta=1$ bands, we will take the bandgap $E_{g}$ as
$E_{g}^{(11)}$, as discussed earlier. The two–state excitonic matrix elements
are defined as[60, 50, 52]

$Q_{nm}^{\alpha}=i\int\,\psi_{cv\mathbf{k}}^{\left(n\right)*}\left[\psi_{cv\mathbf{k}}^{\left(m\right)}\right]_{;k_{\alpha}}d^{2}\mathbf{k},$
(31)

where $\left[\psi_{cv\mathbf{k}}^{\left(m\right)}\right]_{;k_{\alpha}}$ is the
generalized derivative[51] in the $\alpha$–direction of the exciton wave
function for the state $m$ given in terms of the Berry connection
$\Omega_{ij\mathbf{k}}^{\alpha}$, defined as[60]

$\left[\psi_{cv\mathbf{k}}^{\left(m\right)}\right]_{;k_{\alpha}}=\frac{\partial\psi_{cv\mathbf{k}}^{\left(m\right)}}{\partial
k_{\alpha}}-i\left(\Omega_{cc\mathbf{k}}^{\alpha}-\Omega_{vv\mathbf{k}}^{\alpha}\right)\psi_{cv\mathbf{k}}^{\left(m\right)}.$
(32)

As the excitonic wave function will be independent of $k_{z}$, the
$\frac{\partial}{\partial k_{z}}\psi_{cv\mathbf{k}}^{\left(m\right)}$ term is
dropped, meaning that $Q_{nm}^{z}$ reads[52]

$\displaystyle Q_{nm}^{z}$
$\displaystyle=\int\,\psi_{cv\mathbf{k}}^{\left(n\right)*}\left(\Omega_{cc\mathbf{k}}^{z}-\Omega_{vv\mathbf{k}}^{z}\right)\psi_{cv\mathbf{k}}^{\left(m\right)}d^{2}\mathbf{k}.$
(33)

###  IV.1 In–Plane Linear Response

Figure 5: (Left) Real part of the linear IP (top) and OOP (bottom) optical
response for $\epsilon_{\mathrm{medium}}=6.9$, with
$\hbar\Gamma_{p_{-}}=\hbar\Gamma_{d_{-}}=1.3\,\mathrm{meV}$ and
$\hbar\Gamma_{s}=0.4\,\mathrm{meV}$. Orange curve corresponds to the excitonic
bound states, while blue line also includes continuum states. Vertical dashed
line represents the bandgap $E_{g}^{(11)}$ of the system, dot-dashed line
represents $2V=110\,\mathrm{meV}$. In the IP legend, $\gamma_{3}$ and
$\gamma_{0}$ identify which matrix elements allow for the transition in
question. (Right) Diagram of dominant excitonic selection rules in the
$\tau=1$ valley for linear IP (top) and OOP (bottom) optical response.
Transitions allowed with/without skew coupling are shown in dashed/solid
lines.

Recalling Eq. (29), the excitonic selection rules for IP linear response are
obtained directly by considering the angular integral

$\displaystyle\int_{0}^{2\pi}e^{i\ell_{n}\theta}P_{vc\mathbf{k}}^{x}d\theta,$
(34)

where $P_{vc\mathbf{k}}^{x}$ is as defined in Eq. (9). Direct inspection of
this integral tells us that the terms proportional to $\gamma_{0}$ will lead
to contributions from states with $\ell_{n}=\pm\tau$, while the presence of
skew coupling $\gamma_{3}$ allows transitions to states where $\ell_{n}=\pm
2\tau$. This means that the $X_{0n}^{x}$ matrix element can be written, in a
somewhat abusive but concise form, as

$\displaystyle X_{0n}^{x}$
$\displaystyle=X_{\ell_{n}+\tau=0}^{x}+X_{\ell_{n}-\tau=0}^{x}$
$\displaystyle\quad-i\frac{\gamma_{3}}{\gamma_{0}}\left[X_{\ell_{n}+2\tau=0}^{x}-X_{\ell_{n}-2\tau=0}^{x}\right],$
(35)

where the new indices restrict each term to the Kronecker $\delta$’s resulting
from the different angular integrals. Taking $\left|X_{0n}^{x}\right|^{2}$ as
in Eq. (28), the presence of the various Kronecker $\delta$’s originating from
each of the angular integrals means that non–vanishing contributions from
$\gamma_{3}$ to the linear response must be proportional to
$\left|\gamma_{3}\right|^{2}$. Therefore, the inclusion of an $i$ factor in
the Hamiltonian will not alter the results when compared to previous
studies[21]. It will, however, affect the nonlinear response as it will be
linear in $\gamma_{3}$ (in this simple approximation) and, therefore, the $i$
factor proves to be important.

Figure 6: Real part of the linear IP (top) and OOP (bottom) optical response
for $\epsilon=6.9$, with
$\hbar\Gamma_{p_{-}}=\hbar\Gamma_{d_{-}}=1.3\,\mathrm{meV}$ and
$\hbar\Gamma_{s}=0.4\,\mathrm{meV}$, for values of the bias potential $V$
between $20$ and $100\,\mathrm{meV}$. Vertical dashed lines represent the
bandgap $E_{g}^{(11)}$ of the system for each value of the bias potential $V$,
dot-dashed lines represent the gap at $k=0$, equal to $2V$.

Following from the results of Ref. [19], we will assume the broadening to be
dependent on the excitonic state $n$. Specifically, we set the broadening as
dependent on the angular momentum series, given by
$\hbar\Gamma_{p_{-}}=1.3\,\mathrm{meV}$ and
$\hbar\Gamma_{s}=0.4\,\mathrm{meV}$. With the broadening set, we plot the IP
linear response for $V=55\,\mathrm{meV}$ in the top panel of Fig. 5. Varying
the bias potential $V$, we can observe the tunability of the IP linear
response in the top panel of Fig. 6. Similarly to what was observed in Refs.
[21, 22], the dominant peak for the IP linear response quickly increases as
the bias potential increases. Additionally, the location of the peak
corresponding to the $1s$ resonance approaches that of the dominant $2p_{-}$
resonance, as indicated by the evolution of the binding energies in Fig. 4.

###  IV.2 In–Plane Nonlinear Response

Now for the nonlinear optical response, we recall Eqs. (31-32) and write the
integrand present in the two–state momentum matrix element as

$\displaystyle\psi_{cv\mathbf{k}}^{\left(n\right)*}\left[\frac{\partial\psi_{cv\mathbf{k}}^{\left(m\right)}}{\partial
k_{x}}-i\left(\Omega_{cc\mathbf{k}}^{x}-\Omega_{vv\mathbf{k}}^{x}\right)\psi_{cv\mathbf{k}}^{\left(m\right)}\right].$
(36)

Focusing on the angular integration of this expression, and recalling the
definition of the Berry connection present in Eq. (12), one obtains that only
transitions with $\ell_{m}-\ell_{n}=\pm 1$ are allowed. In a similar notation
to Eq. (35), this means that the Berry connection up to zeroth order in
$\frac{\gamma_{3}}{\gamma_{0}}$ can be written as

$Q_{nm}^{x}=Q_{\left|\ell_{m}-\ell_{n}\right|=1}^{x}.$ (37)

Multiplication of the angular momentum Kronecker $\delta$’s from Eqs. (35-37)
leads to four non–zero matrix elements which will be linear in
$\frac{\gamma_{3}}{\gamma_{0}}$. We denote these non–zero matrix elements by
the oscillator strength defined, for compactness, as

$\displaystyle\sigma_{\ell_{n};\ell_{m}}$ $\displaystyle\equiv
X_{\ell_{n}}^{x}Q_{\ell_{n},\ell_{m}}^{x}X_{\ell_{m}}^{x}.$ (38)

Explicitly, the transitions allowed by Eqs. (34-36) will be

$\displaystyle\sigma_{s;p_{-}}$ $\displaystyle\sigma_{p_{-};s}$
$\displaystyle\sigma_{f_{-};g_{-}}$ $\displaystyle\sigma_{g_{-};f_{-}},$ (39)

where the left column corresponds to those where $\ell_{m}-\ell_{n}=-1$ and
the right column to those where $\ell_{m}-\ell_{n}=1$.

Careful inspection of the matrix elements tells us that $\sigma_{s;p_{-}}$ and
$\sigma_{p_{-};s}$ will be the dominant contributions to the SHG nonlinear
response up to linear order in $\frac{\gamma_{3}}{\gamma_{0}}$. The real part
of the SHG nonlinear response is then plotted in the top panel of Fig. 7.
Similarly to previous results[52], we see that the dominant contribution will
originate from the excitonic states with the smallest angular momentum.

###  IV.3 Out–of–Plane Linear Response

Recalling Eq. (29), the dipole transition amplitude for the excitonic state
$n$ with angular momentum quantum number $\ell_{n}$ reads

$\displaystyle X_{0n}^{z}$
$\displaystyle=\int_{0}^{\infty}\frac{f_{cvk}^{(n)}}{E_{cvk}}kdk\int_{0}^{2\pi}e^{i\ell_{n}\theta}P_{vc\mathbf{k}}^{z}d\theta$
(40)

which immediately leads to the optical selection rule $\ell_{n}=0$ for the
linear response when one recalls the definition of $P_{vc\mathbf{k}}^{z}$ from
Eq. (15). Inspecting this equation also tells us that the OOP linear optical
response will have a quadratic dependence on the interlayer spacing,
originating directly from
$\left|\int_{0}^{2\pi}e^{i\ell_{n}\theta}P_{vc\mathbf{k}}^{z}\,d\theta\right|^{2}$.
Additionally, a more complex dependence on $d$ will be implicitly present in
the excitonic wave function $f_{cvk}^{(n)}$ due to the exponential
$e^{-d\left|\mathbf{k}-\mathbf{q}\right|}$ in the form factor present in Eq.
(27).

The real part of the OOP linear excitonic optical response is plotted in the
bottom panel of Fig. 5, together with the diagram of the allowed transitions
in the $\tau=1$ valley. Although it is not labelled in the figure due to its
much weaker amplitude, the $4d_{-}$ transition is still allowed. Contrary to
what was computed for the OOP linear response of a buckled monolayer[52], we
observe that the conductivity quickly stabilizes to a constant value for
$\hbar\omega>E_{g}^{(11)}$, while in a buckled monolayer an apparently linear
growth with $\hbar\omega$ was observed in the same regime. The observed
behaviour also matches with the results for the free–carrier conductivity,
present in the bottom panel of Fig. 2, in the region where the two–band
approximation is valid. Additionally, in the bottom panel of Fig. 6, we
present the OOP optical response for various values of the bias potential $V$
between $20$ and $100\,\mathrm{meV}$. In this figure, we can see the much
larger effect of the bias potential on the OOP optical response in comparison
to the IP response. This is most noticeable between $V=20\,\mathrm{meV}$ and
$55\,\mathrm{meV}$, where the maximum of the linear response for the IP and
OOP linear response increases by factors of approximately $3\times$ and
$5\times$, respectively.

###  IV.4 Out–of–Plane Nonlinear Response

Figure 7: (Left) Real part of the SHG optical response of BBG with diagonal IP
(top), diagonal OOP (middle) and non–diagonal OOP (bottom) conductivity for
$\epsilon_{\mathrm{medium}}=6.9$ and $\hbar\Gamma=1.5\,\mathrm{meV}$. Orange
curve corresponds to only excitonic bound states, while blue line also
includes continuum states. Vertical (dotted) dashed black lines represent
(half) the bandgap $E_{g}^{(11)}$ of the system. (Right) Diagram of dominant
excitonic selection rules in the $\tau=1$ valley for each component.
Transitions allowed with/without skew coupling are shown in dashed/solid
lines. Arrow direction and colour represent the specific resonance when
multiple contributions are present.

Focusing now on the OOP nonlinear regime, we recall the expression for
$Q_{nm}^{z}$ from Eq. (33) as

$\displaystyle
Q_{nm}^{z}=\int\psi_{cv\mathbf{k}}^{\left(n\right)*}\left(\Omega_{cc\mathbf{k}}^{z}-\Omega_{vv\mathbf{k}}^{z}\right)\psi_{cv\mathbf{k}}^{\left(m\right)}d^{2}\mathbf{k}.$
(41)

With this expression in mind, careful analysis of Eq. (16) leads immediately
to the selection rules

$\displaystyle\ell_{n}$ $\displaystyle=\ell_{m}.$ (42)

This selection rule is rather unrestrictive, meaning that it will be
$X_{0n}^{z}$ which will determine the allowed transitions.

Focusing first on the diagonal $\sigma_{zzz}$ response, the selection rules
follow immediately from analysis of $X_{0n}^{z}$ and $Q_{nm}^{z}$ as

$\displaystyle\ell_{n}$ $\displaystyle=\ell_{m}=0.$ (43)

Recalling the pseudo–spin factor, this means that the transitions present will
be associated with $d_{-}$–series states. Knowing the selection rules, we can
compute the SHG conductivity, plotted in the middle panel of Fig. 7. The
matrix elements $X_{0n}^{z}$ (around $1/5$ the magnitude of $X_{0n}^{x}$) and
$Q_{nm}^{z}$ (around $1/4$ the magnitude of the analogous matrix element in
$Q_{nm}^{x}$) are significantly smaller than those present in
$\sigma_{xxx}^{\mathrm{SHG}}$. This, in conjunction with the cubic
dependence[52] on the ratio $d/a\approx 0.7$, means that the peak amplitude of
$\sigma_{zzz}^{\mathrm{SHG}}$ is close to $800$ times smaller than that of
$\sigma_{xxx}^{\mathrm{SHG}}$. A big difference can also be observer in the
qualitative behaviour of $\sigma_{zzz}^{\mathrm{SHG}}$ versus its analogous
counterpart in a buckled monolayer[52]: while in the buckled monolayer the
$\sigma_{zzz}^{\mathrm{SHG}}$ response above the bandgap remains close to its
maximum, in BBG we observe that it quickly tends to a much smaller value,
similarly to what is observed in the $\sigma_{xxx}^{\mathrm{SHG}}$ response.

Secondly, we focus our attention on the off-diagonal $\sigma_{zxx}$ response.
While the presence of $X_{0n}^{z}$ means immediately that $\ell_{n}=0$, we
must now carefully analyze Eqs. (9-34), as well as Eqs. (12-36). From
$Q_{nm}^{x}$, we immediately get that

$\left|\ell_{m}-\ell_{n}\right|=1\Rightarrow\ell_{m}=\pm 1,$

both allowed by the selection rules for $X_{m0}^{x}$. Hence, we are restricted
to states

$\displaystyle\ell_{n}$ $\displaystyle=0,$ $\displaystyle\ell_{m}$
$\displaystyle=\pm 1.$

Again recalling the pseudo–spin factor, we can immediately associate the
selection rules with transitions between $d_{-}$ states and both $p_{-}$ and
$f_{-}$ states. Careful inspection of the oscillator strength for the two
distinct selection rules reveals that, as one would intuitively expect from
the lower angular momentum, resonances associated with the matrix elements
$\sigma_{d_{-};p_{-}}$ dominate. The off-diagonal
$\sigma_{zxx}^{\mathrm{SHG}}$ SHG response is then plotted in the bottom panel
of Fig. 7.

Figure 8: Real part of the SHG IP (top), diagonal OOP (middle) and nondiagonal
OOP (bottom) optical response for $\epsilon=6.9$, with broadening
$\hbar\Gamma=1.5\,\mathrm{meV}$, for values of the bias potential $V$ between
$20$ and $100\,\mathrm{meV}$. Vertical dashed lines represent the bandgap
$E_{g}^{(11)}$ of the system for each value of the bias potential $V$, dot-
dashed lines represent the gap at $k=0$, equal to $2V$. Leftmost vertical
lines correspond to response at $2\hbar\omega$, rightmost vertical lines
correspond to response at $\hbar\omega$.

Finally, in Fig. 8, we present the SHG IP (top), diagonal OOP (middle) and
nondiagonal OOP (bottom) optical response for three different values of the
bias potential $V$. In this figure, we keep a constant broadening
$\hbar\Gamma=1.5\,\mathrm{meV}$ as to more easily observe the effects of
tuning the bias potential. For the lowest value of the bias potential,
$V=20\,\mathrm{meV}$, the lowest bandgap $E_{g}^{(11)}=39.8\,\mathrm{meV}$ is
very close to $2V=40\,\mathrm{meV}$, we only keep the vertical lines
associated with $E_{g}^{(11)}$ and $E_{g}^{(11)}/2$ for improved readability.
Starting with the IP response $\sigma_{xxx}^{\mathrm{SHG}}$, we can observe
that, at $V=20\,\mathrm{meV}$, the proximity between the multiple resonances
leads to a sharp increase in intensity when compared with the response at
$V=55\,\mathrm{meV}$ and $100\,\mathrm{meV}$. For the diagonal OOP response
$\sigma_{zzz}^{\mathrm{SHG}}$, we observe that the magnitude of the excitonic
response quickly grows with the bias potential $V$. This follows what is
expected, as it is the OOP anisotropy due to the bias potential that leads to
a non–zero OOP response. Lastly, the nondiagonal OOP response
$\sigma_{zxx}^{\mathrm{SHG}}$ also presents an increasing magnitude with the
bias potential $V$, as expected from the growing OOP anisotropy. The nature of
the nondiagonal response, with its mixing of IP and OOP components, means that
the response at low bias potential (see $V=20\,\mathrm{meV}$ regime) does not
vanish as quickly as the $\sigma_{zzz}^{\mathrm{SHG}}$ component.
Additionally, direct comparison between top and bottom panels of Fig. 8 shows
that, as bias potential increases, the magnitude of the nondiagonal OOP
response approaches that of the IP response, with the maximum of
$\sigma_{zxx}^{\mathrm{SHG}}$ approximately $1/2$ that of
$\sigma_{xxx}^{\mathrm{SHG}}$ for $V=100\,\mathrm{meV}$.

##  V Summary

In this paper, we studied the excitonic linear and nonlinear optical response
of Bernal stacked BBG as a function of the gate voltage, both for in–plane
(IP) and out–of–plane (OOP) directions. Starting with the electronic structure
of BBG, we consider the influence of gate voltage on the band structure of the
system. Taking into account the three distinct bandgaps present, we introduce
a two–band approximation that greatly simplifies the numerical complexity of
the problem. Outlining the generic form of the momentum matrix elements and
Berry connections, we compute the linear IP and OOP free–carrier optical
response for both the two–band approximation and the full four–band
Hamiltonian. The large separation of the higher energy bands (for
$V=55\,\mathrm{meV}$, the lowest bandgap is $E_{g}^{(11)}=106\,\mathrm{meV}$
while the second lowest is $E_{g}^{(12)}=436\,\mathrm{meV}$) means that, for
frequency below $E_{g}^{(12)}$, the real part of the optical response is very
accurately described by the two–band approximation.

Under this two–band approximation, we then computed the excitonic states in
the system by numerical diagonalization of the Bethe–Salpeter equation,
studying the evolution of the excitonic binding energies as a function of
external bias potential.

Knowing the excitonic states, together with the momentum matrix elements and
Berry connections, allowed us to then explicitly discuss the excitonic
selection rules of the system for both the IP and OOP excitonic response. The
inclusion of skew coupling in the momentum matrix elements proves fundamental
in obtaining a non–zero IP nonlinear response, as well as introducing
$s-$series resonances in the IP linear response. For the OOP linear response
the OOP anisotropy proves to be sufficient for obtaining a non–zero linear and
nonlinear response.

Both linear and SHG nonlinear response are very sensitive to the tuning of the
bias voltage. Firstly, tuning the bias potential will directly affect the
bandgap of the system, leading to a tuning of the excitonic binding energies
and, therefore, the location of the resonant peaks in the optical response, as
well as on the interband momentum matrix elements. Additionally, as the OOP
anisotropy is controlled by the bias potential, its effects on the OOP
response are not as simple as those on the IP response.

## Acknowledgements

M.F.C.M.Q. acknowledges the International Iberian Nanotechnology Laboratory
(INL) and the Portuguese Foundation for Science and Technology (FCT) for the
Quantum Portugal Initiative (QPI) grant SFRH/BD/151114/2021. N.M.R.P.
acknowledges support by the Portuguese Foundation for Science and Technology
(FCT) in the framework of the Strategic Funding UIDB/04650/2020, COMPETE 2020,
PORTUGAL 2020, FEDER, and FCT through projects PTDC/FIS-MAC/2045/2021,
EXPL/FIS-MAC/0953/2021. N.M.R.P. also acknowledges the Independent Research
Fund Denmark (grant no. 2032-00045B) and the Danish National Research
Foundation (Project No. DNRF165).

## References

  * [1] K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, Y. Zhang, S. V. Dubonos, I. V. Grigorieva, and A. A. Firsov, “Electric Field Effect in Atomically Thin Carbon Films,” Science, vol. 306, pp. 666–669, 2004.
  * [2] A. Taghizadeh, K. S. Thygesen, and T. G. Pedersen, “Two-dimensional materials with giant optical nonlinearities near the theoretical upper limit,” ACS Nano, vol. 15, pp. 7155–7167, 2021.
  * [3] J. D. Caldwell, I. Aharonovich, G. Cassabois, J. H. Edgar, B. Gil, and D. N. Basov, “Photonics with hexagonal boron nitride,” Nat. Rev. Mater., vol. 4, pp. 552–567, 2019.
  * [4] G. Wang, A. Chernikov, M. M. Glazov, T. F. Heinz, X. Marie, T. Amand, and B. Urbaszek, “Colloquium: Excitons in atomically thin transition metal dichalcogenides,” Rev. Mod. Phys., vol. 90, p. 021001, 2018.
  * [5] H. Sahin, S. Cahangirov, M. Topsakal, E. Bekaroglu, E. Akturk, R. T. Senger, and S. Ciraci, “Monolayer honeycomb structures of group-IV elements and III-V binary compounds: First-principles calculations,” Phys. Rev. B, vol. 80, p. 155453, 2009.
  * [6] H. K. Avetissian and G. F. Mkrtchian, “Higher harmonic generation by massive carriers in buckled two-dimensional hexagonal nanostructures,” Phys. Rev. B, vol. 99, p. 085432, 2019.
  * [7] M.-Q. Le, H.-T. Nguyen, and T.-L. Bui, “Fracture of 28 buckled two-dimensional hexagonal sheets,” Mech. Adv. Mater. Struct., vol. 29, pp. 4993–5005, 2022\.
  * [8] R. Y. Kezerashvili and A. Spiridonova, “Effects of parallel electric and magnetic fields on Rydberg excitons in buckled two-dimensional materials,” Phys. Rev. B, vol. 103, p. 165410, 2021.
  * [9] A.-Y. Lu, H. Zhu, J. Xiao, C.-P. Chuu, Y. Han, M.-H. Chiu, C.-C. Cheng, C.-W. Yang, K.-H. Wei, Y. Yang, Y. Wang, D. Sokaras, D. Nordlund, P. Yang, D. A. Muller, M.-Y. Chou, X. Zhang, and L.-J. Li, “Janus monolayers of transition metal dichalcogenides,” Nat. Nanotechnol., vol. 12, pp. 744–749, 2017\.
  * [10] J. Zhang, S. Jia, I. Kholmanov, L. Dong, D. Er, W. Chen, H. Guo, Z. Jin, V. B. Shenoy, L. Shi, and J. Lou, “Janus monolayer transition-metal dichalcogenides,” ACS Nano, vol. 11, pp. 8192–8198, 2017.
  * [11] T. Zheng, Y.-C. Lin, Y. Yu, P. Valencia-Acuna, A. A. Puretzky, R. Torsi, C. Liu, I. N. Ivanov, G. Duscher, D. B. Geohegan, Z. Ni, K. Xiao, and H. Zhao, “Excitonic dynamics in janus MoSSe and WSSe monolayers,” Nano Lett., vol. 21, pp. 931–937, 2021.
  * [12] O. Dogadov, C. Trovatello, B. Yao, G. Soavi, and G. Cerullo, “Parametric nonlinear optics with layered materials and related heterostructures,” Laser Photonics Rev., vol. 16, p. 2100726, 2022.
  * [13] J. Shi, H. Xu, C. Heide, C. HuangFu, C. Xia, F. de Quesada, H. Shen, T. Zhang, L. Yu, A. Johnson, F. Liu, E. Shi, L. Jiao, T. Heinz, S. Ghimire, J. Li, J. Kong, Y. Guo, and A. M. Lindenberg, “Giant room-temperature nonlinearities in a monolayer janus topological semiconductor,” Nature Communications, vol. 14, no. 1, p. 4953, 2023.
  * [14] C.-H. Park and S. G. Louie, “Tunable excitons in biased bilayer graphene,” Nano Lett., vol. 10, pp. 426–431, 2010.
  * [15] Y. Gong, J. Lin, X. Wang, G. Shi, S. Lei, Z. Lin, X. Zou, G. Ye, R. Vajtai, B. I. Yakobson, H. Terrones, M. Terrones, B.-K. Tay, J. Lou, S. T. Pantelides, Z. Liu, W. Zhou, and P. M. Ajayan, “Vertical and in-plane heterostructures from WS2/MoS2 monolayers,” Nat. Mater., vol. 13, pp. 1135–1142, 2014.
  * [16] P. Rivera, K. L. Seyler, H. Yu, J. R. Schaibley, J. Yan, D. G. Mandrus, W. Yao, and X. Xu, “Valley-polarized exciton dynamics in a 2D semiconductor heterostructure,” Science, vol. 351, pp. 688–691, 2016.
  * [17] S. J. Brun and T. G. Pedersen, “Intense and tunable second-harmonic generation in biased bilayer graphene,” Phys. Rev. B, vol. 91, p. 205405, 2015.
  * [18] C.-H. Park and S. G. Louie, “Tunable Excitons in Biased Bilayer Graphene,” Nano Lett., vol. 10, pp. 426–431, 2010.
  * [19] L. Ju, L. Wang, T. Cao, T. Taniguchi, K. Watanabe, S. G. Louie, F. Rana, J. Park, J. Hone, F. Wang, and P. L. McEuen, “Tunable excitons in bilayer graphene,” Science, vol. 358, pp. 907–910, 2017.
  * [20] A. A. Avetisyan, A. P. Djotyan, and K. Moulopoulos, “Tunable excitons in bilayer graphene with opened energy gap,” Phys. At. Nucl., vol. 81, no. 6, pp. 799–803, 2018.
  * [21] J. C. G. Henriques, I. Epstein, and N. M. R. Peres, “Absorption and optical selection rules of tunable excitons in biased bilayer graphene,” Phys. Rev. B, vol. 105, p. 045411, 2022.
  * [22] M. O. Sauer and T. G. Pedersen, “Exciton absorption, band structure, and optical emission in biased bilayer graphene,” Phys. Rev. B, vol. 105, p. 115416, 2022.
  * [23] H. D. Scammell and O. P. Sushkov, “Dynamical screening and excitonic bound states in biased bilayer graphene,” Phys. Rev. B, vol. 107, p. 085104, 2023\.
  * [24] T. Ando and M. Koshino, “Optical absorption by interlayer density excitations in bilayer graphene,” J. Phys. Soc. Jpn., vol. 78, p. 104716, 2009.
  * [25] M.-E. Kleemann, R. Chikkaraddy, E. M. Alexeev, D. Kos, C. Carnegie, W. Deacon, A. C. de Pury, C. Große, B. de Nijs, J. Mertens, A. I. Tartakovskii, and J. J. Baumberg, “Strong-coupling of WSe2 in ultra-compact plasmonic nanocavities at room temperature,” Nat. Commun., vol. 8, p. 1296, 2017\.
  * [26] M. Stührenberg, B. Munkhbat, D. G. Baranov, J. Cuadra, A. B. Yankovich, T. J. Antosiewicz, E. Olsson, and T. Shegai, “Strong light–matter coupling between plasmons in individual gold bi-pyramids and excitons in mono- and multilayer WSe2,” Nano Lett., vol. 18, pp. 5938–5945, 2018.
  * [27] A. Dasgupta, J. Gao, and X. Yang, “Atomically thin nonlinear transition metal dichalcogenide holograms,” Nano Lett., vol. 19, pp. 6511–6516, 2019.
  * [28] G. Soavi, G. Wang, H. Rostami, D. G. Purdie, D. De Fazio, T. Ma, B. Luo, J. Wang, A. K. Ott, D. Yoon, S. A. Bourelle, J. E. Muench, I. Goykhman, S. Dal Conte, M. Celebrano, A. Tomadin, M. Polini, G. Cerullo, and A. C. Ferrari, “Broadband, electrically tunable third-harmonic generation in graphene,” Nat. Nanotechnol., vol. 13, pp. 583–588, 2018.
  * [29] S. Klimmer, O. Ghaebi, Z. Gan, A. George, A. Turchanin, G. Cerullo, and G. Soavi, “All-optical polarization and amplitude modulation of second-harmonic generation in atomically thin semiconductors,” Nat. Photon., vol. 15, pp. 837–842, 2021.
  * [30] Y. Zhang, Y. Wang, Y. Dai, X. Bai, X. Hu, L. Du, H. Hu, X. Yang, D. Li, Q. Dai, T. Hasan, and Z. Sun, “Chirality logic gates,” Sci. Adv., vol. 8, p. eabq8246, 2022.
  * [31] Y. Li, N. An, Z. Lu, Y. Wang, B. Chang, T. Tan, X. Guo, X. Xu, J. He, H. Xia, Z. Wu, Y. Su, Y. Liu, Y. Rao, G. Soavi, and B. Yao, “Nonlinear co-generation of graphene plasmons for optoelectronic logic operations,” Nat. Commun., vol. 13, p. 3138, 2022.
  * [32] Z. Zheng, K. Chang, and J. L. Cheng, “Gate voltage induced injection and shift currents in aa- and ab-stacked bilayer graphene,” Phys. Rev. B, vol. 108, p. 235401, 2023.
  * [33] T. Ohta, A. Bostwick, T. Seyller, K. Horn, and E. Rotenberg, “Controlling the electronic structure of bilayer graphene,” Science, vol. 313, no. 5789, pp. 951–954, 2006.
  * [34] F. Guinea, A. H. Castro Neto, and N. M. R. Peres, “Electronic states and landau levels in graphene stacks,” Phys. Rev. B, vol. 73, p. 245426, 2006\.
  * [35] T. Stauber, N. M. R. Peres, and A. K. Geim, “Optical conductivity of graphene in the visible region of the spectrum,” Phys. Rev. B, vol. 78, p. 085432, 2008.
  * [36] A. Rozhkov, A. Sboychakov, A. Rakhmanov, and F. Nori, “Electronic properties of graphene-based bilayer systems,” Phys. Rep., vol. 648, pp. 1–104, 2016\.
  * [37] C. Ortix, “Nonlinear hall effect with time-reversal symmetry: Theory and material realizations,” Adv. Quantum Technol., vol. 4, no. 9, p. 2100056, 2021.
  * [38] E. V. Castro, K. S. Novoselov, S. V. Morozov, N. M. R. Peres, J. M. B. L. dos Santos, J. Nilsson, F. Guinea, A. K. Geim, and A. H. C. Neto, “Biased bilayer graphene: Semiconductor with a gap tunable by the electric field effect,” Phys. Rev. Lett., vol. 99, p. 216802, 2007.
  * [39] J. M. B. Lopes dos Santos, N. M. R. Peres, and A. H. Castro Neto, “Graphene bilayer with a twist: Electronic structure,” Phys. Rev. Lett., vol. 99, p. 256802, 2007.
  * [40] Y. Cao, V. Fatemi, S. Fang, K. Watanabe, T. Taniguchi, E. Kaxiras, and P. Jarillo-Herrero, “Unconventional superconductivity in magic-angle graphene superlattices,” Nature, vol. 556, no. 7699, pp. 43–50, 2018.
  * [41] M. Martins Quintela, J. Guerra, and S. João, “Electronic properties of twisted bilayer graphene in the presence of a magnetic field,” EPJ Web Conf., vol. 233, p. 03004, 2020.
  * [42] E. Y. Andrei and A. H. MacDonald, “Graphene bilayers with a twist,” Nat. Mater., vol. 19, no. 12, pp. 1265–1275, 2020.
  * [43] L. Di Mauro Villari and A. Principi, “Optotwistronics of bilayer graphene,” Phys. Rev. B, vol. 106, p. 035401, 2022.
  * [44] H. Yu, X. Cui, X. Xu, and W. Yao, “Valley excitons in two-dimensional semiconductors,” Nat. Sci. Rev., vol. 2, pp. 57–70, 2015.
  * [45] E. McCann and M. Koshino, “The electronic properties of bilayer graphene,” Rep. Prog. Phys., vol. 76, no. 5, p. 056503, 2013.
  * [46] S. Di Sabatino, J. A. Berger, and P. Romaniello, “Optical spectra of 2D monolayers from time-dependent density functional theory,” Faraday Discuss., vol. 224, pp. 467–482, 2020.
  * [47] X. Zhang, W.-Y. Shan, and D. Xiao, “Optical Selection Rule of Excitons in Gapped Chiral Fermion Systems,” Phys. Rev. Lett., vol. 120, p. 077401, 2018.
  * [48] T. Cao, M. Wu, and S. G. Louie, “Unifying Optical Selection Rules for Excitons in Two Dimensions: Band Topology and Winding Numbers,” Phys. Rev. Lett., vol. 120, p. 087402, 2018.
  * [49] J. C. G. Henriques, B. Amorim, R. M. Ribeiro, and N. M. R. Peres, “Excitonic response of $\mathrm{AA^{\prime}}$ and AB stacked hBN bilayers,” Phys. Rev. B, vol. 105, p. 115421, 2022.
  * [50] A. Taghizadeh and T. G. Pedersen, “Nonlinear optical selection rules of excitons in monolayer transition metal dichalcogenides,” Phys. Rev. B, vol. 99, p. 235433, 2019.
  * [51] C. Aversa and J. E. Sipe, “Nonlinear optical susceptibilities of semiconductors: Results with a length-gauge analysis,” Phys. Rev. B, vol. 52, pp. 14636–14645, 1995.
  * [52] M. F. C. M. Quintela and T. G. Pedersen, “Anisotropic linear and nonlinear excitonic optical properties of buckled monolayer semiconductors,” Phys. Rev. B, vol. 107, p. 235416, 2023.
  * [53] K. Yoshizawa, T. Yumura, T. Yamabe, and S. Bandow, “The role of orbital interactions in determining the interlayer spacing in graphite slabs,” J. Am. Chem. Soc., vol. 122, no. 48, pp. 11871–11875, 2000.
  * [54] M. S. Alam, J. Lin, and M. Saito, “First-principles calculation of the interlayer distance of the two-layer graphene,” Jpn. J. Appl. Phys., vol. 50, no. 8R, p. 080213, 2011.
  * [55] B. Butz, C. Dolle, F. Niekiel, K. Weber, D. Waldmann, H. B. Weber, B. Meyer, and E. Spiecker, “Dislocations in bilayer graphene,” Nature, vol. 505, no. 7484, pp. 533–537, 2014.
  * [56] Y. Fukaya, Y. Zhao, H.-W. Kim, J. R. Ahn, H. Fukidome, and I. Matsuda, “Atomic arrangements of quasicrystal bilayer graphene: Interlayer distance expansion,” Phys. Rev. B, vol. 104, p. L180202, 2021.
  * [57] V. N. Genkin and P. M. Mednis, “Contribution to the theory of nonlinear effects in crystals with account taken of partially filled bands,” Sov. Phys. JETP, vol. 27, p. 609, 1968.
  * [58] B. Kirtman, F. L. Gu, and D. M. Bishop, “Extension of the Genkin and Mednis treatment for dynamic polarizabilities and hyperpolarizabilities of infinite periodic systems. I. Coupled perturbed Hartree–Fock theory,” J. Chem. Phys, vol. 113, pp. 1294–1309, 2000.
  * [59] V. A. Margulis, E. E. Muryumin, and E. A. Gaiduk, “Optical second-harmonic generation from two-dimensional hexagonal crystals with broken space inversion symmetry,” J. Phys. Condens. Matter, vol. 25, p. 195302, 2013\.
  * [60] T. G. Pedersen, “Intraband effects in excitonic second-harmonic generation,” Phys. Rev. B, vol. 92, p. 235432, 2015.
  * [61] M. F. C. Martins Quintela and N. M. R. Peres, “Tunable excitons in rhombohedral trilayer graphene,” Phys. Rev. B, vol. 105, p. 205417, 2022\.
  * [62] M. F. C. M. Quintela, J. C. G. Henriques, L. G. M. Tenório, and N. M. R. Peres, “Theoretical methods for excitonic physics in 2D materials,” Phys. Status Solidi B, vol. 259, p. 2200097, 2022.
  * [63] G. F. Glinskii and Z. Koinov, “Functional formulation of microscopic theory of exciton polaritons,” Theor. Math. Phys., vol. 70, no. 3, pp. 252–261, 1987\.
  * [64] S. K. Radha, W. R. L. Lambrecht, B. Cunningham, M. Grüning, D. Pashov, and M. van Schilfgaarde, “Optical response and band structure of LiCoO2 including electron-hole interaction effects,” Phys. Rev. B, vol. 104, p. 115120, 2021.
  * [65] S. N. Rytova, “The screened potential of a point charge in a thin film,” Mosc. Un. Phys. Bul., vol. 3, pp. 30–37, 1967.
  * [66] L. V. Keldysh, “Coulomb interaction in thin semiconductor and semimetal films,” Sov. J. Exp. and Theor. Phys. Lett., vol. 29, p. 658, 1979.
  * [67] P. Li and I. Appelbaum, “Excitons without effective mass: Biased bilayer graphene,” Phys. Rev. B, vol. 99, p. 035429, 2019.
  * [68] T. Tian, D. Scullion, D. Hughes, L. H. Li, C.-J. Shih, J. Coleman, M. Chhowalla, and E. J. G. Santos, “Electronic polarizability as the fundamental variable in the dielectric properties of two-dimensional materials,” Nano Lett., vol. 20, pp. 841–851, 2020.
  * [69] A. Laturia, M. L. Van de Put, and W. G. Vandenberghe, “Dielectric properties of hexagonal boron nitride and transition metal dichalcogenides: from monolayer to bulk,” npj 2D Mater. Appl., vol. 2, p. 6, 2018.
  * [70] C. Y.-P. Chao and S. L. Chuang, “Analytical and numerical solutions for a two-dimensional exciton in momentum space,” Phys. Rev. B, vol. 43, pp. 6530–6543, 1991.
  * [71] D. G. W. Parfitt and M. E. Portnoi, “The two-dimensional hydrogen atom revisited,” J. Math. Phys., vol. 43, pp. 4681–4691, 2002.
  * [72] X. L. Yang, S. H. Guo, F. T. Chan, K. W. Wong, and W. Y. Ching, “Analytic solution of a two-dimensional hydrogen atom. i. nonrelativistic theory,” Phys. Rev. A, vol. 43, pp. 1186–1196, 1991.


# Intelligent Reflecting Surface Aided Wireless Energy Transfer and Mobile
Edge Computing for Public Transport Vehicles

Shan Jaffry  School of Internet of Things,
Xi’an Jiaotong-Liverpool University.
<EMAIL_ADDRESS>
#  Intelligent Reflecting Surface Aided
Wireless Energy Transfer and Mobile Edge Computing for Public Transport
Vehicles
–
A Communication Eco-System For Transport Vehicles of 6G-era.

Shan Jaffry  School of Internet of Things,
Xi’an Jiaotong-Liverpool University.
<EMAIL_ADDRESS>
###### Abstract

In the forthcoming 6G era, Future Public transport vehicles (F-PTV), such as
buses, trains etc. will be designed to cater the communication needs of the
commuters that will carry numerous smart devices (smartphones, e-bands etc.).
These battery-powered devices need frequent recharging. Since recharging
facilities are not readily available while commuting, we envision F-PTVs that
will provide in-vehicle recharging via Wireless Energy Transfer (WET) through
in-vehicle Access Points. F-PTV will also be internally coated with
Intelligent Reflecting Surface (IRS) that reflect incident radio waves towards
the intended device to improve signal strength at the receiver, for both
information and energy transmissions. F-PTVs will also be equipped with Mobile
Edge Computing (MEC) servers to also serve multiple purposes, including
reduction in devices’ task completion latency and offering in-vehicle Cloud
services. MEC-server offloading will also relieve smart devices’ processors
off intensive tasks to preserve battery power. The challenges associated with
the IRS-aided integrated MEC-WET model within F-PTV are also discussed. We
present simulation to show the effectiveness of such a model for F-PTV.

###### Index Terms:

6G, Intelligent Reflective Surface, Wireless Energy Transfer, Wireless Power
Transfer, Mobile Edge Computing.

##  I Introduction

Modern applications running on smart devices and gadgets are increasingly
becoming sophisticated and consume excessive computation resources and battery
power. On the other hand, the available computational and battery resources
are limited for such devices, mainly due to their small form factor.
Replenishing devices’ battery levels on-the-go is becoming a serious
challenge. In particular, mobile users that are on the move, for example,
commuters inside public transport vehicles (PTVs), such as buses, trains,
subways etc., are more vulnerable to battery drainage as they have limited
recharging facilities available. These users collectively carry a large number
of devices and engage in online or offline activities during the commute. Such
activities (e.g. audio/video calling, gaming, web-surfing, etc.) drain
extensive amount of devices’ battery.

On the other hand, it has been proposed that future PTVs (F-PTVs) will have a
dedicated in-vehicle small cell to provide services to commuting users [1].
Primarily, Access Points within F-PTV (PTV-AP) are aimed to enhance in-vehicle
cellular user’s quality of service (QoS). The commuters within F-PTV will be
effectively decoupled from the out-of-vehicle macro-cell base station which
significantly reduce cellular network load. Furthermore, F-PTV will also
reduce the group handover into a collective hand-off by the vehicle [2].
Integrating in-vehicle wireless energy transfer (WET) methods with PTV-APs can
enable wireless charging of commuters’ energy-constrained devices and gadgets
while on the move.

However, due to hostility of air interface towards electromagnetic radio
signals, the magnitude of power delivered though WET remains effectively very
low. To make matter worse, signal quality degrades drastically in a crowded
public transport as signal passes through a human bodies. This body-shadowing
effect further diminishes the prospect of any energy gains with WET-enabled
PTV-AP. To bolster wireless channels within F-PTV, we further envision
installing the Intelligent Reflecting Surface (IRS) within the vehicle. The
IRS are passive reflecting elements that change the phase and amplitude of
wireless signals and can effectively undo the signal attenuation by reflecting
incident waves towards the receiver device.

Figure 1: A depiction of a communication-friendly futuristic Public Transport
Vehicle for 6G era(An idea).

On the other hand, to reduce increasing computing requirements from the smart
devices, researchers have focused on offloading services and tasks to the
near-by external servers which are relatively more powerful than smart devices
and gadgets. These powerful computing nodes, also known as MEC-servers, are
located at the network edge, for example with the base-stations and APs. These
MEC-servers can be collocated with PTV-AP to provide in-vehicle Cloud services
and augment computing requirements of commuting users inside F-PTV via task
offloading. The task offloading can be exploited for dual, albeit closely
linked, purposes. First, offloading reduce tasks’ processing time as MEC-
server’s Central Processing Units (CPUs) are more powerful than smart devices.
Second, the task offloading also preserves extensive amount of device’s power
that is otherwise consumed during the local computations.

This envisioned IRS-aided MEC-WET system creates a cellular-friendly eco-
system within F-PTV in the 6G era. In the rest of this article, we first
provide a brief overview of the envisioned combined IRS-aided integrated MEC-
WET system for F-PTVs. We then discuss potential challenges and possible
solutions associated with the envisioned F-PTV. Later, Simulation and
numerical results are provided to validate the effectiveness of the design.

##  II IRS-Aided MEC-WET within Public Transport Vehicle

In this section we will first discuss about the WET technology followed by IRS
and MEC system within F-PTV.

###  II-A Suitable WET Technology for PTV

WET is not an entirely new idea. In fact, it is more than century old concept
[3]. Famous inventor Nikola Tesla conducted wireless power transfer
experiments in the late 19th century [4]. However, modern improvements in
microelectronic circuit design have sparked a new interest in WET for smart
devices, gadgets, and sensor networks.

We broadly categorized WET technologies into two types. The first type, called
Inductive Coupling WET (IC-WET), exploits near-field electromagnetic energy
transmission phenomena. IC-WET is popular for charging smart gadgets, electric
vehicles, implanted medical devices etc. More recently, key smart device
manufacturers have built practical wireless charging accessories based on IC-
WET. However, IC-WET works only for very short distances and in most cases the
charging adapter and the rechargeable device needs to be within near-physical
contact. Due to this limitation, IC-WET is not a suitable fit for F-PTV use
case.

The second type of WET uses radio frequency (RF) signals to transmit energy
and offers flexibility of charging devices at a relatively larger distances
(few meters). This makes RF-WET a more suitable option for the use case within
F-PTV. However, electromagnetically radiated RF signals withers significantly
in-air before reaching the receiver. The highly attenuated signal received at
the harvesting device is significantly low, and even ineffective at times, to
provide any useful energy gains. The signal attenuation exacerbates within
F-PTV because of human-body shadowing due to presence of passengers.

###  II-B IRS Coating within Vehicle

The poor wireless propagation environment within the F-PTV can be alleviated
by internally coating passive IRS elements within the vehicle. An IRS can
virtually reconfigure RF propagation environment by controlling the phase and
amplitude of incident waves [5]. Interestingly, unlike legacy relaying
systems, IRS do not require external power source to boost signal strength.
Instead, IRS comprise of array of low-cost passive reflecting elements that
can be controlled in real-time to focus the direction of transmitted wave
towards the receiving device, much like a mirror. An in-vehicle IRS controller
is responsible for adjusting amplitude and phase of incident waves. The waves
can be added constructively to significantly improve the signal strength at
the desired receiver. The IRS-aided constructive interference can be achieved
both for information and energy harvesting signals. The same mechanism can be
used to induce destructive interference at unwanted or eavesdropping receivers
to nullify the undesirable interference effect.

Due to passive nature of IRS elements, additional power for signal decoding or
amplification is not required. Furthermore, unlike active relays IRS also do
not impose additional white noise as a consequence of decoding or
amplification [6]. The IRS’s RF reflecting elements can be conveniently coated
over the existing structures such as on the walls or windows of the vehicle
which significantly reduces the deployment complexity [5].

To fully exploit IRS benefits within F-PTV, controlling reflective elements is
a challenging task. The IRS controller must have real-time knowledge of
factors such as channel state information (CSI), number of wireless links,
optimal phase and amplitude etc [7]. However, even a fixed or sub-optimal
tuning of reflectors have shown to significantly improve the quality of
signal-to-noise-ratio (SNR) or harvested energy at the receiver [8, 6].

Despite IRS-aided improved channel conditions, effective energy harvesting
gain for the user devices remain limited due to low RF-to-DC power conversion
efficiency and other hardware constraints. Hence, next we discuss MEC
offloading to further improve energy preservation within F-PTV.

###  II-C Edge Computing for Energy Preservation

When it comes to WET, RF-to-DC power conversion efficiency is practically much
lower than theoretically anticipated [9]. This is in part true due to device’s
size limitation and hardware energy conversion constraints. Hence only limited
energy gain could be achieved by relying solely on WET, even after IRS
installation. On the other hand, a significant amount of devices’ energy
consumption is due to processing requirements of applications’ tasks. This is
true because a device’s CPU’s power consumption is proportional to the
operating frequency and the operational voltage of the CPU. By frequently
forcing device’s CPU into idle mode, significant energy savings can be
achieved. One way of relieving device’s CPU off continuous operations is by
offloading tasks to the MEC-servers collocated with PTV-AP. This MEC
offloading serves two purposes. Firstly, it allows faster completion of tasks
as CPU(s) at MEC-servers are more powerful than user devices. Secondly, it
preserves a significant amount of the device’s battery power, courtesy to task
offloading. Furthermore in [8] researchers have shown that IRS-aided wireless
environment also reduce system-level latency for MEC systems by improving SNR
at the end receivers and also overall improves energy efficiency.

The IRS-aided MEC-WET system within F-PTV does not come without challenges
which are discussed next.

##  III Design Challenges and Future Research Directions

In this section we discuss some of the challenges associated with envisioned
in-vehicle IRS-aided MEC-WET system.

###  III-A Crowded Channels within PTV

Due to relatively higher device-density (devices per square meter) within
public transport and installation of reflectors within F-PTV, in-vehicle radio
channels may get over-crowded in case excessive information and energy signals
are exchanged between PTV-AP and commuters. Even in the absence of online
activity, users may compete for channels to harvest RF energy or to upload and
download tasks to and from the MEC-servers, respectively. Given the
insufficient number of available channels for both information and energy
transmission, fair channel scheduling will be challenging for IRS aided
environment within PTV.

Alternatively, the mmWave band, which is suitable for indoor communication,
offers wider bandwidth and more radio channels to alleviate congestion
problems. However, due to large propagation losses, mmWave band may not be
suitable for energy signal transmission. One possible solution to resolve
above stated issues is to use a combination of mmWave and sub-6 GHz (e.g. 2.4
GHz and 3.5 GHz) bands within F-PTV such that information signals are carried
over mmWave frequencies and energy harvesting is carried over sub-6 GHz band.

###  III-B Limited MEC-WET Resources within PTV

As compared to Cloud storage, MEC servers have relatively limited storage and
computing resources. In particular, MEC servers within PTVs will have further
reduced capabilities due to physical size limitations. Hence, optimally
distributing the limited computing resources among several in-vehicle devices
would be another challenge.

Furthermore, energy signal transmission will base on time-sharing scheme. For
a given Transmission Time Internal (TTI) $\tau$, energy signal is transmitted
for the period $\tau_{v}=\gamma\tau$, where $\gamma\in[0,1]$. The remaining
time $\tau_{i}=(1-\gamma)\tau$, is dedicated for information transmission
($\tau_{i}>\tau_{v}$). Given that $\tau_{v}$ is the small proportion of
$\tau$, allocating slots for numerous devices need efficient time scheduling.
PTV-APs with Multiple-Input-Multiple-Output (MIMO) antenna may allow
simultaneous transmission of multiple energy signals via spatial diversity.
However, scheduling energy resources for large number of users is still a
challenge.

To address computing and energy resource management and scheduling, smart
optimization techniques are needed. One simplistic scheme could be to rely on
the tasks’ completion deadline ($Z_{n}$) and the battery level ($E_{n}$) for
each device. Based on this information, PTV-AP can allocate MEC and WET
resources for individual user. For example, as shown in Figure 2 if a device
has high battery level and lower deadline constraint, then that the device
task can be performed locally. On the contrary, if a device has low battery
levels then it is given higher priority for both MEC and WET resource
allocation.

Figure 2: A simplistic decision grid for scheduling computing and energy
resources scheduling among devices.

###  III-C Active and Passive Beamforming

The internal structure of most public transport is generally suitable for
coating IRS elements. For example, the roof or the upper sidewalls could be
covered with the reflecting elements. Given the passive nature of IRS elements
additional power is not required to improve signal quality as the IRS
controller adjust the phase and amplitude of incident waves striking the meta-
material surface. The installation of PTV-AP and IRS would follow Line-of-
Sight (LoS) link. Active beamforming at the PTV-AP and passive beamforming at
the IRS elements would guide the direction of the traveling waves. As
mentioned earlier, the presence of large number of passengers would increase
body-shadowing effect within PTV. In the downlink direction, the combined
active and passive beamforming can alleviate this effect. However, active
beamforming may not always be possible for the uplink transmission due to
device (such as smartphones, watches, e-bands etc.) limitations. Still,
researchers in [6, 10] have shown that in an IRS aided environment, even the
signals transmitted in random direction provide better quality than non-IRS
environment.

###  III-D Electromagnetic Radiation Concerns

Regular in-vehicle transmissions and additional IRS reflections inside public
transport increase electromagnetic radiation inside the vehicle which raises
health concerns for passengers. Hence, these concerns must be addressed by
suitable researchers in line of standardized guidelines outlined by the
governing bodies.

To demonstrate the feasibility of IRS-aided MEC-WET installation within PTV,
numerical results are presented in the next Section.

##  IV Simulation Results and Empirical Observation

The model diagram for the IRS-aided MEC-WET system within PTV is shown in
Figure 1. WE consider a PTV-AP that is capable of transmitting data and energy
signals. The MEC server is co-located with the PTV-AP and hence AP-to-server
latency is negligible. The IRS’ reflecting elements are coated on the
vehicle’s interior as shown in Figure 1. A single-antenna PTV-AP and a single
device are considered for the simulation to show the effectiveness of the
system. Figure 3 shows the geometrical representation of PTV-AP, IRS, and
device with default locations and distances.

A practical non-linear energy harvesting model [10] is considered for the
simulation in which the harvesting is successful if the received signal power
is above a threshold $P_{th}$, zero otherwise. The harvested power cannot
exceed an upper threshold limit of $P_{max}$. Parameters considered in [10]
are used in our simulation as well, unless otherwise stated.

Figure 3: Geometrical representation of elements within PTV used for
simulation.

The channel model consist of both small and large-scale fading. In particular,
exponentially distributed channels with zero mean and unit variance are
considered in both downlink and uplink direction. Channels are independent and
identically distributed. The IRS placement is considered such that PTV-AP and
IRS elements follows LoS link with pathloss parameter $\alpha_{o}=2.4$. The
pathloss parameter for PTV-to-User and IRS-to-User channels is
$\alpha_{n}=3.0$.

At the start of time-frame of length $\tau$, user device generates $l$ Kbits
of data per unit of time, where $l$ is distributed uniformly between 40 Kbits
and 50 Kbits. If the device wants to offload its task, it sends a message with
relevant information to the PTV-AP which makes decision regarding offloading
or local computing. In case of offloading, $l_{o}=\mu l$ Kbits of data is
offloaded to MEC server and $l_{c}=(1-\mu)l$ Kbits is computed locally by the
device, where $\mu\in[0,1]$. Device’s CPU takes $\psi=10^{3}$ cycles to
compute one bit of raw data and energy consumed by CPU to compute $l_{c}$ bits
is $\xi\psi^{3}l_{c}^{3}$ mJoule, where $\xi=1\times 10^{-28}$ is the
effective CPU capacitance [10]. After offloading $l_{o}$ Kbits, the device
saves $\xi\psi^{3}l_{o}^{3}$ mJoule energy. We considered PTV-AP transmit
power $P_{t}=1$ Watt. For the IRS, amplitude reflection coefficient is set to
1 [8] and phase shift coefficient ($\Phi$) is set to $\frac{\pi}{2}$. Each
simulation is iterated 1000 times and the averaged results are reported in the
following.

###  IV-A Numerical Results

First we discuss Figure 4 and Figure 5 in which device’s energy levels $E[t]$
are shown with respect to time $t$. Device’s new energy level in the
subsequent time steps is $E[t+1]=E[t]+E_{g}[t+1]-E_{c}[t+1]$, where
$E_{g}[t]=E_{v}[t]+E_{o}[t]$ is the instantaneous energy gain at time $t$ and
$E_{v}[t],E_{o}[t]$ denote the harvested energy and energy preserved due to
offloading, respectively. The energy consumption is calculated as
$E_{c}=E_{ckt}+E_{lc}+E_{tr}$, where $E_{ckt},E_{lc},E_{tr}$ denote device
circuit’s constant energy consumption, energy utilization for local computing,
and energy needed to offload task to MEC server. The mobile device’s energy
level is normalized to 1 mJoule at the start of the simulation (i.e. $t=0$).

Both figures (Fig. 4 and Fig. 5) show the consumption of device’s battery
energy levels under normal scenario. There is some energy gain due to WET,
however it is not very significant particularly due to poor channel condition
between PTV-AP and the device. Even though the IRS aided transmission
significantly improves signal strength at the receiving device, the benefits
reaped from the reflecting surface are still limited due to device-level
saturation [11]. This is because in a practical WET system, a device cannot
harvest energy beyond its saturation limit due to hardware constraints. To
fully exploit the IRS-aided MEC-WET, the maximum allowed harvesting power
threshold at the receiving device should be improved. This problem require
hardware-level solution from the device manufacturers and is out of scope of
theoretical research. As the maximum harvesting power threshold in Figure 4 is
$P_{max}=0.004$ Watts, WET does not offer any significant energy gains. On the
other hand, when $P_{max}$ is increased to 0.04 Watts, moderate WET-enabled
energy gain is observed in Figure 5, even without IRS.

Figure 4: Device’s energy levels with respect to time $P_{max}$= .004 Watts,
$\mu$= 0.75 , $\tau_{v}$= 0.25

Figure 4 and Figure 5 also show that the significant gain in energy occurs due
to MEC offloading. For example, for $\mu=0.75$ and $\tau_{0}.25$, negligible
energy losses are observed for the MEC-only case. This is due to the fact that
MEC server takes the toil of computation as result of task offloading which
allows energy preservation at user-device end. In case of IRS-aided MEC-WET
scenario, the device’s battery gains significantly more energy than that
consumed. In fact, device is able to actually store energy above the starting
level of 1 mJoule. It can be observed in Figure 5 that by increasing number of
reflecting elements from 64 to 200 the gain due to energy harvesting improved
considerably.

In Figure 6 we show the energy gain to consumption ratio for different amount
of data generated by user device. Considering $\mu=0.5$, we can observe that
the ratio for MEC-only scenario and remains below 1 for WET-only and IRS-aided
WET scenarios when data size is greater than 30 Kbits. For smaller data size
(e.g. 24 Kbits), energy gain in case of IRS-aided WET is twice that of energy
consumed by the device. On the other hand, the ratio is significantly higher
for MEC and IRS-aided MEC-WET scenario, particularly when data size are low.
The ratio decreases as the size of data increases because device need more
power for offloading and computing remaining data locally.

Figure 5: Device’s energy levels with respect to time $P_{max}$= .04 Watts,
$\mu$= 0.75 , $\tau_{v}$= 0.25

##  V Conclusion

In this paper we present an idea for the futuristic public transport vehicle
that is equipped with in-vehicle AP, Mobile Edge Computing and Wireless Energy
Transfer facilities and use Intelligent Reflecting Surface to improve in-
vehicle channel condition. Since preserving device battery levels for
commuting users is a serious challenge, future F-PTVs will replenish
commuters’ smart devices by employing IRS aided MEC-WET. The simulation
results for IRS-aided integrated MEC-WET system showed significantly increased
device battery’s energy levels. Some challenges associated with the proposed
F-PTV and their solutions are also presented in this paper. Such a model is
still open for research. As our future work, we are working on jointly
optimizing $\tau_{v},\mu$ and phase shift coefficient $\Phi$ within the
proposed framework, along with designing MEC-WET resource scheduling scheme
for multiple user scenario.

Figure 6: Ratio of energy gain-to-consumption ratio N=64, $P_{max}$ = 0.04
Watt, $N=64$, $\tau_{v}=0.25,\mu=0.5$.

##  VI Acknowledgment

This work is submitted to IEEE Vehicular Technology Magazine for publication.
It is supported by XJTLU’s Research Development Fund: RDF-20-01-15.

## References

  * [1] S. Jaffry, S. F. Hasan, and X. Gui, “Making a case for the moving small cells,” in _26th IEEE International Telecommunication Networks and Applications Conference_ , 2016.
  * [2] S. Jaffry, R. Hussain, X. Gui, and S. F. Hasan, “A comprehensive survey on moving networks,” _IEEE Communications Surveys & Tutorials_, 2020.
  * [3] N. Shinohara, “Power without wires,” _IEEE Microwave magazine_ , vol. 12, no. 7, pp. S64–S73, 2011.
  * [4] N. Tesla, “The transmission of electric energy without wires (the thirteenth anniversary number of the electrical world and engineer),” 1904.
  * [5] E. Basar, M. Di Renzo, J. De Rosny, M. Debbah, M.-S. Alouini, and R. Zhang, “Wireless communications through reconfigurable intelligent surfaces,” _IEEE Access_ , vol. 7, pp. 116 753–116 773, 2019.
  * [6] C. Pan, H. Ren, K. Wang, M. Elkashlan, A. Nallanathan, J. Wang, and L. Hanzo, “Intelligent reflecting surface aided mimo broadcasting for simultaneous wireless information and power transfer,” _IEEE JSAC_ , 2020.
  * [7] Q. Wu and R. Zhang, “Towards smart and reconfigurable environment: Intelligent reflecting surface aided wireless network,” _IEEE Comms. Magazine_ , vol. 58, no. 1, pp. 106–112, 2019.
  * [8] T. Bai, C. Pan, Y. Deng, M. Elkashlan, A. Nallanathan, and L. Hanzo, “Latency minimization for intelligent reflecting surface aided mobile edge computing,” _IEEE JSAC_ , vol. 38, no. 11, pp. 2666–2682, 2020.
  * [9] A. S. Almansouri, M. H. Ouda, and K. N. Salama, “A cmos rf-to-dc power converter with 86% efficiency and- 19.2-dbm sensitivity,” _IEEE Transactions on Microwave Theory and Techniques_ , vol. 66, no. 5, pp. 2409–2415, 2018.
  * [10] F. Zhou and R. Q. Hu, “Computation efficiency maximization in wireless-powered mobile edge computing networks,” _IEEE Trans. on Wireless Communications_ , vol. 19, no. 5.
  * [11] S. Wang, M. Xia, K. Huang, and Y.-C. Wu, “Wirelessly powered two-way communication with nonlinear energy harvesting model: Rate regions under fixed and mobile relay,” _IEEE Transactions on Wireless Communications_ , vol. 16, no. 12, pp. 8190–8204, 2017.