# Gate voltage induced injection and shift currents in AA- and AB-stacked
bilayer graphene

Ze Zheng  GPL Photonics Laboratory, State Key Laboratory of Luminescence and
Applications, Changchun Institute of Optics, Fine Mechanics and Physics,
Chinese Academy of Sciences, Changchun, Jilin, 130033 P. R. China  University
of Chinese Academy of Sciences, Beijing 100039, China.  Kainan Chang
<EMAIL_ADDRESS>GPL Photonics Laboratory, State Key Laboratory of
Luminescence and Applications, Changchun Institute of Optics, Fine Mechanics
and Physics, Chinese Academy of Sciences, Changchun, Jilin, 130033 P. R. China
University of Chinese Academy of Sciences, Beijing 100039, China.  Jin Luo
Cheng<EMAIL_ADDRESS>GPL Photonics Laboratory, State Key Laboratory of
Luminescence and Applications, Changchun Institute of Optics, Fine Mechanics
and Physics, Chinese Academy of Sciences, Changchun, Jilin, 130033 P. R. China
University of Chinese Academy of Sciences, Beijing 100039, China.

###### Abstract

Generating photogalvanic effects in centrosymmetric materials can provide new
opportunities for developing passive photodetectors and energy harvesting
devices. In this work, we investigate the photogalvanic effects in
centrosymmetric two-dimensional materials, AA- and AB-stacked bilayer
graphene, by applying an external gate voltage to break the symmetry. Using a
tight-binding model to describe the electronic states, the injection
coefficients for circular photogalvanic effects and shift conductivities for
linear photogalvanic effects are calculated for both materials with light
wavelengths ranging from THz to visible. We find that gate voltage induced
photogalvanic effects can be very significant for AB-stacked bilayer graphene,
with generating a maximal dc current in the order of mA for a 1 $\mu$m wide
sample illuminated by a light intensity of 0.1 GW/cm2, which is determined by
the optical transition around the band gap and van Hove singularity points.
Although such effects in AA-stacked bilayer graphene are about two orders of
magnitude smaller than those in AB-stacked bilayer graphene, the spectrum is
interestingly limited in a very narrow photon energy window, which is
associated with the interlayer coupling strength. A detailed analysis of the
light polarization dependence is also performed. The gate voltage and chemical
potential can be used to effectively control the photogalvanic effects.

##  I Introduction

Photogalvanic effects are nonlinear optical responses that generate direct
currents in homogeneous materials, and such a passive process is considered as
a direct and powerful photoelectric conversion method [1, 2, 3]. The widely
discussed photogalvanic effects can be induced by the one-color injection
current and shift current, which are second order nonlinear optical processes
occurring in noncentrosymmetric materials, or the two-color coherent current
injection processes, which are third (for “1+2” process) [4] or fifth (for
“2+3” process) [5] order nonlinear optical processes and are not sensitive to
the inversion symmetry of materials. According to the response to the light
polarization, second order photogalvanic effects are also phenomenologically
divided into circularly polarized photogalvanic effect and linearly polarized
photogalvanic effect, where the latter is light phase insensitive and can be
used for solar energy harvest without forming p-n junctions to surpass the
Shockley-Queisser limit [6, 7, 8]. One of the research topics in this field is
to find materials with significant photogalvanic effects at a specific
frequency range, and several studies have been conducted on various new
materials, including 2D materials [9, 10, 11, 12, 13], Dirac or Weyl
semimetals [1, 14, 15], ferroelectric materials [16, 17, 18, 19], and so on.

As the first two-dimension material, graphene is a potential candidate for
realizing new functionality in optoelectronic devices due to its superior
optical and electronic properties exceeding many traditional bulk materials.
However, because of its centrosymmetric crystal structure, one-color injection
and shift currents vanish in many few-layer graphene as well as their
nanostructures, while two-color coherent control has been well studied in both
theories [20, 21, 4, 22] and experiments [23, 24]. It is still meaningful to
generate one-color injection and shift currents in centrosymmetric graphene
based structure, in order to utilize its extraordinary physical properties.
The generation of second order response can be realized by forming an
asymmetric interface or edge [25], applying an external electric field [26],
forming surface curvature [27], considering the spatial variation of the light
field [28], and stacking graphene layers into asymmetric structure [29]. Wei
et al. [9] studied the gate field induced injection and shift currents in
zigzag graphene nanoribbons, and found that the subband and edge states
determine the generated currents with an effective modulation of their
amplitudes by the ribbon width and the static field strength. Xiong et al.
[30] investigated the light polarization dependence of in-plane shift current
in a AB-stacked bilayer graphene (AB-BG) with applying a gate voltage, and
their results clearly illustrated a sizeable photocurrent at a given light
frequency; however, neither the spectra of the shift conductivity nor the
injection current was present. By stacking two layers of monolayer graphene
with a relative rotation to form a twisted bilayer graphene, a large shift
current can be produced due to a huge density of states when the flat band is
formed at magic angles [13, 12, 31]. Surprisingly, whether the gate voltage
can generate photogalvanic effect in AA-stacked bilayer graphene (AA-BG) is
still not clear.

In this paper, we systematically study the spectra of the injection
coefficients and shift conductivities of AA-BG and AB-BG under applying a gate
voltage to break the inversion symmetry, as well as their dependence on the
gate voltage and chemical potential. Their electronic states are described by
widely adopted tight-binding model formed by the carbon $2p_{z}$ orbitals [26,
32], and the expressions for injection coefficient and shift conductivity are
employed from Ref. [33]. Our results confirm the feasibilities of generating
photogalvanic effects in AA-BG and AB-BG. Particularly, the response of AA-BG
distributes in a very narrow spectral region, while a maximal current in the
order of mA can be generated in AB-BG for a 1 $\mu$m wide sample at light
intensity of 0.1 GW/cm2.

This paper is organized as follows. In Sec. II we introduce the tight-binding
models for the AA-BG and AB-BG under applying a gate voltage, and give the
expressions for the injection coefficient and shift conductivity. In Sec. III
we present the spectra of injection coefficient and shift conductivity for AA-
BG and AB-BG, and discuss the effects of the gate voltage and chemical
potential. We conclude in Sec. IV.

##  II Models

###  II.1 Hamiltonian

Figure 1: Crystal structures and tight-binding hopping parameters for (a) AA-
BG and (b) AB-BG.

We consider the tight-binding Hamiltonian for the AA-BG and AB-BG, whose
crystal structures are illustrated in Fig. 1 (a) and (b), respectively. These
two structures have the same primitive lattice vectors
$\bm{a}_{1}=a_{0}\left(\frac{1}{2}\hat{\bm{x}}+\frac{\sqrt{3}}{2}\hat{\bm{y}}\right)$
and
$\bm{a}_{2}=a_{0}\left(-\frac{1}{2}\hat{\bm{x}}+\frac{\sqrt{3}}{2}\hat{\bm{y}}\right)$
with the lattice constant $a_{0}=2.46$ Å. The atomic positions in the unit
cell are taken as $\bm{\tau}_{A}=\bm{0}$,
$\bm{\tau}_{B}=(\bm{a}_{1}+\bm{a}_{2})/3$,
$\bm{\tau}_{A^{\prime}}=c\hat{\bm{z}}$, and
$\bm{\tau}_{B^{\prime}}=\bm{\tau}_{B}+c\hat{\bm{z}}$ for AA-BG, and
$\bm{\tau}_{A}=\bm{0}$, $\bm{\tau}_{B}=(\bm{a}_{1}+\bm{a}_{2})/3$,
$\bm{\tau}_{A^{\prime}}=\bm{\tau}_{B}+c\hat{\bm{z}}$, and
$\bm{\tau}_{B^{\prime}}=2\bm{\tau}_{B}+c\hat{\bm{z}}$ for AB-BG, where
$c=3.35$ Å is the interlayer distance. The primitive reciprocal lattice
vectors are
$\bm{b}_{1}=\frac{2\pi}{a_{0}}\left(\hat{\bm{x}}+\frac{1}{\sqrt{3}}\hat{\bm{y}}\right)$
and
$\bm{b}_{2}=\frac{2\pi}{a_{0}}\left(-\hat{\bm{x}}+\frac{1}{\sqrt{3}}\hat{\bm{y}}\right)$.
The electronic states are described by a tight-binding model employing carbon
$2p_{z}$ orbitals. The unperturbed Hamiltonian [32] for AA-BG is

$\displaystyle H^{\rm
AA}_{\bm{k}}=\left(\begin{array}[]{cccc}-\Delta&\gamma_{0}g_{\bm{k}}&\gamma_{1}&\gamma_{3}g_{\bm{k}}\\\
\gamma_{0}g^{*}_{\bm{k}}&-\Delta&\gamma_{3}g^{*}_{\bm{k}}&\gamma_{1}\\\
\gamma_{1}&\gamma_{3}g_{\bm{k}}&\Delta&\gamma_{0}g_{\bm{k}}\\\
\gamma_{3}g^{*}_{\bm{k}}&\gamma_{1}&\gamma_{0}g^{*}_{\bm{k}}&\Delta\\\
\end{array}\right)\,.$ (5)

Here $\bm{k}$ is the electron wavevector, and
$g_{\bm{k}}=1+e^{-i\bm{k}\cdot\bm{a}_{1}}+e^{-i\bm{k}\cdot\bm{a}_{2}}$. The
hopping parameters are illustrated in Fig. 1 (a) with $\gamma_{0}=2.569$ eV,
$\gamma_{1}=0.361$ eV, and $\gamma_{3}=-0.032$ eV. The on-site energies
$\pm\Delta$ are induced by a gate voltage. The Hamiltonian for AB-BG is given
from Ref. 26 as

$\displaystyle H^{\rm
AB}_{\bm{k}}=\left(\begin{array}[]{cccc}-\Delta-\frac{\Delta^{\prime}}{2}&\gamma_{0}^{\prime}g_{\bm{k}}&\gamma_{4}^{\prime}g_{\bm{k}}&\gamma_{3}^{\prime}g^{*}_{\bm{k}}\\\
\gamma_{0}^{\prime}g^{*}_{\bm{k}}&-\Delta+\frac{\Delta^{\prime}}{2}&\gamma_{1}^{\prime}&\gamma_{4}^{\prime}g_{\bm{k}}\\\
\gamma_{4}^{\prime}g^{*}_{\bm{k}}&\gamma_{1}^{\prime}&\Delta+\frac{\Delta^{\prime}}{2}&\gamma_{0}^{\prime}g_{\bm{k}}\\\
\gamma_{3}^{\prime}g_{\bm{k}}&\gamma_{4}^{\prime}g^{*}_{\bm{k}}&\gamma_{0}^{\prime}g^{*}_{\bm{k}}&\Delta-\frac{\Delta^{\prime}}{2}\\\
\end{array}\right)\,,$ (10)

where the hopping parameters (see Fig. 1 (b)) are $\gamma_{0}^{\prime}=-3.16$
eV, $\gamma_{1}^{\prime}=0.381$ eV, $\gamma_{3}^{\prime}=-0.38$ eV, and
$\gamma_{4}^{\prime}=0.14$ eV. The on-site potential difference
$\Delta^{\prime}=0.022$ eV is induced by the asymmetric environment of A, B
atoms in the crystal structure.

The eigenstates $C_{n\bm{k}}$ and eigenenergies $\epsilon_{n\bm{k}}$ at the
$n$th band are obtained by diagonalizing the Hamiltonian through

$\displaystyle H_{\bm{k}}C_{n\bm{k}}=\epsilon_{n\bm{k}}C_{n\bm{k}}\,.$ (11)

The calculation of the optical responses involves the position operator
$\widetilde{\bm{r}}_{\bm{k}}$ and velocity operator
$\widetilde{\bm{v}}_{\bm{k}}$, which are

$\displaystyle\widetilde{\bm{r}}_{\bm{k}}$
$\displaystyle=i\bm{\nabla}_{\bm{k}}+\left(\begin{array}[]{cccc}\bm{\tau}_{A}&0&0&0\\\
0&\bm{\tau}_{B}&0&0\\\ 0&0&\bm{\tau}_{A^{\prime}}&0\\\
0&0&0&\bm{\tau}_{B^{\prime}}\end{array}\right)\,,\quad\widetilde{\bm{v}}_{\bm{k}}=\frac{1}{i\hbar}[\widetilde{\bm{r}}_{\bm{k}},H_{\bm{k}}]\,,$
(16)

respectively. The matrix elements of the position operator give the Berry
connections $\bm{\xi}_{nm\bm{k}}$ by

$\displaystyle\bm{\xi}_{nm\bm{k}}=C^{\dagger}_{n\bm{k}}\widetilde{\bm{r}}_{\bm{k}}C_{m\bm{k}}\,,$
(17)

and those of the velocity operator are calculated as
$\bm{v}_{nm\bm{k}}=C^{\dagger}_{n\bm{k}}\widetilde{\bm{v}}_{\bm{k}}C_{m\bm{k}}$.
Due to the derivative with respect to the wavevector $\bm{k}$, a direct
calculation of $\bm{\xi}_{nm\bm{k}}$ from Eq. (17) requires that the
wavefunction $C_{n\bm{k}}$ is a smooth function of $\bm{k}$. However, this
becomes quite difficult in numerical calculation because of the phase
arbitrary for a numerical wavefunction. Practically, the off-diagonal terms of
$\bm{\xi}_{nm\bm{k}}$ can be also calculated from the velocity operator as

$\displaystyle\bm{r}_{nm\bm{k}}=\begin{cases}\bm{\xi}_{nm\bm{k}}=\frac{\bm{v}_{nm\bm{k}}}{i\omega_{nm\bm{k}}}&(n\neq
m)\\\ 0&(n=m)\end{cases}\,,$ (18)

with $\hbar\omega_{nm\bm{k}}=\epsilon_{n\bm{k}}-\epsilon_{m\bm{k}}$. The
diagonal terms $\xi^{a}_{nn\bm{k}}$ usually appear in the generalized
derivative of $(r^{c}_{\bm{k}})_{;nmk^{a}}=\frac{\partial
r_{nm\bm{k}}^{c}}{\partial
k^{a}}-i(\xi^{a}_{nn\bm{k}}-\xi^{a}_{mm\bm{k}})r^{c}_{nm\bm{k}}$, which is
calculated alternatively [9] by

$\displaystyle(r^{c}_{\bm{k}})_{;nmk^{a}}=$
$\displaystyle\frac{-ir^{c}_{nm\bm{k}}{\cal V}_{mn\bm{k}}^{a}+\hbar
M^{ca}_{nm\bm{k}}+i[r^{a}_{\bm{k}},v^{c}_{\bm{k}}]_{nm}}{i\omega_{nm\bm{k}}}\,,$
(19)

with ${\cal
V}_{mn\bm{k}}^{a}=v^{a}_{mm\bm{k}}-v^{a}_{nn\bm{k}}=\frac{\partial\omega_{mn\bm{k}}}{\partial
k^{a}}$ and

$\displaystyle
M^{ca}_{nm\bm{k}}=C^{\dagger}_{n\bm{k}}\frac{1}{i\hbar}[\widetilde{r}^{a}_{\bm{k}},\widetilde{v}^{c}_{\bm{k}}]C_{m\bm{k}}\,,$
(20)

where the Raman letters $a,c$ indicate the Cartesian directions $x,y,z$. Note
that the electron wavevector has only in-plane components $x,y$, the
derivative $\frac{\partial}{\partial k^{z}}$ thus gives zero and
$\left(r_{\bm{k}}^{a}\right)_{;nmk^{z}}=-i(\xi^{z}_{nn\bm{k}}-\xi^{z}_{mm\bm{k}})r^{a}_{nm\bm{k}}$.

###  II.2 Injection and shift currents

We focus on the injection and shift currents induced by a laser pulse centered
at frequency $\omega$, for which the electric field is
$\bm{E}(t)=\bm{E}_{0}(t)e^{-i\omega t}+c.c.$ and $\bm{E}_{0}(t)$ is a slow
varying envelop function. The response static currents can be written as

$\displaystyle\bm{J}_{0}(t)=\bm{J}_{\rm inj}(t)+\bm{J}_{\rm sh}(t)\,.$ (21)

Here the first term $\bm{J}_{\rm inj}(t)$ is a one-color injection current
satisfying

$\displaystyle\frac{dJ^{a}_{\rm
inj}(t)}{dt}=2i\eta^{abc}(\omega)E_{0}^{b}(t)[E_{0}^{c}(t)]^{*}\,,$ (22)

with the injection coefficient $\eta^{abc}(\omega)$ given by

$\displaystyle\eta^{abc}(\omega)=\frac{2e^{3}\pi}{\hbar^{2}}\int\frac{d\bm{k}}{4\pi^{2}}\sum_{nm}{\cal
V}_{mn\bm{k}}^{a}f_{nm\bm{k}}\text{Im}[r^{c}_{mn\bm{k}}r^{b}_{nm\bm{k}}]\delta(\omega_{mn\bm{k}}-\omega)\,.$
(23)

Here $f_{nm\bm{k}}=f_{n\bm{k}}-f_{m\bm{k}}$ is the population difference with
the Fermi-Dirac distribution
$f_{n\bm{k}}=[1-e^{(\epsilon_{n\bm{k}}-\mu)/k_{B}T}]^{-1}$ for given chemical
potential $\mu$ and temperature $T$. The second term $\bm{J}_{\rm sh}(t)$ in
Eq. (21) is a shift current written as

$\displaystyle J^{a}_{\rm sh}(t)$
$\displaystyle=2\sigma^{abc}(\omega)E_{0}^{b}(t)[E_{0}^{c}(t)]^{*}\,,$ (24)

with the shift conductivity $\sigma^{abc}(\omega)$ given by

$\displaystyle\sigma^{abc}(\omega)$ $\displaystyle=-\frac{i\pi
e^{3}}{\hbar^{2}}\int\frac{d\bm{k}}{4\pi^{2}}\sum_{nm}f_{nm\bm{k}}\left[r^{b}_{mn\bm{k}}\left(r_{\bm{k}}^{c}\right)_{;nmk^{a}}+r^{c}_{mn\bm{k}}\left(r_{\bm{k}}^{b}\right)_{;nmk^{a}}\right]\delta(\omega_{mn\bm{k}}-\omega)\,.$
(25)

Further discussion of photocurrents starts with a symmetry analysis on the
tensors of $\eta^{abc}(\omega)$ and $\sigma^{abc}(\omega)$. The presence of
time-reversal symmetry gives
$\bm{r}_{nm\bm{k}}=\bm{r}_{mn(-\bm{k})}=[\bm{r}_{nm(-\bm{k})}]^{*}$,
$\bm{v}_{nm\bm{k}}=-\bm{v}_{mn(-\bm{k})}=-[\bm{v}_{nm(-\bm{k})}]^{*}$,
$\epsilon_{n\bm{k}}=\epsilon_{n(-\bm{k})}$, and
$\left(r^{b}_{\bm{k}}\right)_{;nmk^{a}}=-\left(r^{b}_{-\bm{k}}\right)_{;mnk^{a}}=-[\left(r^{b}_{\bm{k}}\right)_{;nmk^{a}}]^{*}$.
Thus from Eqs. (23) and (25), we obtain $\eta^{abc}=[\eta^{abc}]^{*}$ and
$\sigma^{abc}=[\sigma^{abc}]^{*}$, which are both real numbers. At finite gate
voltage, the crystal point group of AB-BG is $C_{3v}$, whose symmetry is lower
than that of AA-BG with crystal point group $C_{6v}$. Thus we can check the
symmetry properties of AB-BG first, and then refine them to AA-BG. Combining
the point group and the time reversal symmetry, the nonzero tensor components
satisfy $\eta^{xzx}=\eta^{yzy}=\eta^{xxz}=\eta^{yyz}$,
$\sigma^{xzx}=\sigma^{yzy}=\sigma^{xxz}=\sigma^{yyz}$,
$\sigma^{zxx}=\sigma^{zyy}$, $\sigma^{zzz}$, and
$\sigma^{yyy}=-\sigma^{yxx}=-\sigma^{xxy}=-\sigma^{xyx}$. Then the injection
current becomes

$\displaystyle\frac{dJ_{\rm
inj}^{a}(t)}{dt}=4\eta^{xzx}(\omega)\text{Im}\\{E^{a}_{0}(t)[E^{z}_{0}(t)]^{\ast}\\}(1-\delta_{a,z})\,,$
(26)

and the shift current is

$\displaystyle J_{\rm sh}^{x}(t)$
$\displaystyle=4\sigma^{xzx}(\omega)\text{Re}\left\\{E^{z}_{0}(t)[E^{x}_{0}(t)]^{*}\right\\}-4\sigma^{yyy}(\omega)\text{Re}\left\\{E^{x}_{0}(t)[E^{y}_{0}(t)]^{*}\right\\}\,,$
(27a) $\displaystyle J_{\rm sh}^{y}(t)$
$\displaystyle=4\sigma^{xzx}(\omega)\text{Re}\left\\{E^{z}_{0}(t)[E^{y}_{0}(t)]^{*}\right\\}+2\sigma^{yyy}(\omega)\left[|E^{y}_{0}(t)|^{2}-|E^{x}_{0}(t)|^{2}\right]\,,$
(27b) $\displaystyle J_{\rm sh}^{z}(t)$
$\displaystyle=2\sigma^{zxx}(\omega)\left[|E^{x}_{0}(t)|^{2}+|E^{y}_{0}(t)|^{2}\right]+2\sigma^{zzz}(\omega)|E^{z}_{0}(t)|^{2}\,.$
(27c)

For AA-BG, the results are similar except that the $\sigma^{yyy}$ component
disappears due to the extra crystal symmetry.

The injection current in AA-BG or AB-BG requires an elliptically polarized
light incident obliquely, and its $z$-component vanishes due to the lack of
freely moving electrons along this quantum confined direction. The
$z$-component of shift current in AA-BG or AB-BG, induced by the charge shift
between the two layers under the light excitation, can be always generated.
Such shift current can lead to charge accumulation between these two layers,
which can further induce a gate voltage in this system, as discussed by Gao et
al. [34]. The in-plane components of the shift current in AA-BG can be
generated only for an elliptically polarized light incident obliquely, while
those in AB-BG have no such limit.

##  III Results

###  III.1 Analytical results for AA-BG

The Hamiltonian for the AA-BG can be analytically diagonalized. The
eigenstates are

$\displaystyle C_{n\bm{k}}$ $\displaystyle=\frac{\sqrt{1-\alpha_{n}{\cal
N}_{\beta_{n}\bm{k}}}}{2\sqrt{2}}\left(\begin{array}[]{c}-\hat{g}_{\bm{k}}\\\
-\beta_{n}\\\ \beta_{n}\hat{g}_{\bm{k}}\\\
1\end{array}\right)+\frac{\alpha_{n}\sqrt{1+\alpha_{n}{\cal
N}_{\beta_{n}\bm{k}}}}{2\sqrt{2}}\left(\begin{array}[]{c}\hat{g}_{\bm{k}}\\\
\beta_{n}\\\ \beta_{n}\hat{g}_{\bm{k}}\\\ 1\end{array}\right)\,,$ (36)

with $\hat{g}_{\bm{k}}={g_{\bm{k}}}/{|g_{\bm{k}}|}$ and

$\displaystyle{\cal
N}_{\beta_{n}\bm{k}}=\frac{\gamma_{3}|g_{\bm{k}}|+\beta_{n}\gamma_{1}}{\sqrt{\Delta^{2}+(\gamma_{3}|g_{\bm{k}}|+\beta_{n}\gamma_{1})^{2}}}\,.$
(37)

Here $n=1,2,3,4$ denotes the band index with $\alpha_{n}=-1,-1,+1,+1$ and
$\beta_{n}=-1,+1,-1,+1$, respectively. The associated eigenenergies are

$\displaystyle\epsilon_{n\bm{k}}=\beta_{n}\gamma_{0}|g_{\bm{k}}|+\alpha_{n}\sqrt{\Delta^{2}+(\gamma_{3}|g_{\bm{k}}|+\beta_{n}\gamma_{1})^{2}}\,.$
(38)

With the analytic wavefunctions in Eq. (36), Berry connections
$\bm{\xi}_{nm\bm{k}}$ can be calculated directly from Eq. (17), as listed in
Appendix A, where the relations between all components are also presented.
There exist selection rules for $r_{nm\bm{k}}^{z}$ as

$\displaystyle r^{z}_{13\bm{k}}$ $\displaystyle=r^{z}_{31\bm{k}}=\frac{c{\cal
N}_{-1\bm{k}}}{2}\,,\quad r^{z}_{24\bm{k}}=r^{z}_{42\bm{k}}=\frac{c{\cal
N}_{+1\bm{k}}}{2}\,.$ (39)

Therefore, $r_{nm\bm{k}}^{z}$ is nonzero only for the band pair $(n,m)=(1,3)$
or $(2,4)$. The injection coefficient becomes

$\displaystyle\eta^{xzx}(\omega)=$
$\displaystyle\frac{e^{3}}{2\pi\hbar^{2}}\int d\bm{k}\left\\{f_{13\bm{k}}{\cal
V}_{31\bm{k}}^{x}\text{Im}[r^{x}_{31\bm{k}}r^{z}_{13\bm{k}}]\delta(\omega_{31\bm{k}}-\omega)\right.$
$\displaystyle\left.\hskip 62.59596pt+f_{24\bm{k}}{\cal
V}_{42\bm{k}}^{x}\text{Im}[r^{x}_{42\bm{k}}r^{z}_{24\bm{k}}]\delta(\omega_{42\bm{k}}-\omega)\right\\}\,.$
(40)

The intraband Berry connections are obtained as

$\displaystyle\bm{\xi}_{nn\bm{k}}$
$\displaystyle=\frac{1}{2}\left[g_{\bm{k}}^{\ast}(i\bm{\nabla}_{\bm{k}})g_{\bm{k}}+\frac{a_{0}}{\sqrt{3}}\hat{\bm{y}}\right]+\frac{1}{2}c\hat{\bm{z}}\left(1+\alpha_{n}\sqrt{1-{\cal
N}_{{\beta}_{n}\bm{k}}^{2}}\right)\,,$ (41)

The matrix elements for $\xi_{nn\bm{k}}^{x/y}$ are independent of the band
index $n$, thus $(r_{\bm{k}}^{a})_{;nmk^{b}}=\frac{\partial
r_{nm\bm{k}}^{a}}{\partial k^{b}}$ for $b=x,y$ and
$(r_{\bm{k}}^{a})_{;nmk^{z}}=-i(\xi_{nn\bm{k}}^{z}-\xi_{mm\bm{k}}^{z})r_{nm\bm{k}}^{a}$.
The shift conductivities become

$\displaystyle\sigma^{xzx}(\omega)=$
$\displaystyle-i\frac{e^{3}}{4\pi\hbar^{2}}\int
d\bm{k}\left[f_{13\bm{k}}\left(r^{z}_{31\bm{k}}\frac{\partial
r^{x}_{13\bm{k}}}{\partial k_{x}}+r^{x}_{31\bm{k}}\frac{\partial
r^{z}_{13\bm{k}}}{\partial
k_{x}}\right)\delta(\omega_{31\bm{k}}-\omega)\right.$
$\displaystyle\left.\hskip
79.6678pt+f_{24\bm{k}}\left(r^{z}_{42\bm{k}}\frac{\partial
r^{x}_{24\bm{k}}}{\partial k_{x}}+r^{x}_{42\bm{k}}\frac{\partial
r^{z}_{24\bm{k}}}{\partial
k_{x}}\right)\delta(\omega_{42\bm{k}}-\omega)\right]\,,$ (42a)
$\displaystyle\sigma^{zzz}(\omega)=$
$\displaystyle\frac{e^{3}}{2\pi\hbar^{2}}\int
d\bm{k}\left[f_{12\bm{k}}|r^{z}_{31\bm{k}}|^{2}(\xi_{33\bm{k}}^{z}-\xi_{11\bm{k}}^{z})\delta(\omega_{31\bm{k}}-\omega)\right.$
$\displaystyle\left.\hskip
59.75095pt+f_{24\bm{k}}|r^{z}_{42\bm{k}}|^{2}(\xi_{44\bm{k}}^{z}-\xi_{22\bm{k}}^{z})\delta(\omega_{42\bm{k}}-\omega)\right]\,,$
(42b) $\displaystyle\sigma^{zxx}(\omega)=$
$\displaystyle\frac{e^{3}}{2\pi\hbar^{2}}\int
d\bm{k}\sum_{nm}f_{nm\bm{k}}|r^{x}_{mn\bm{k}}|^{2}(\xi_{mm\bm{k}}^{z}-\xi_{nn\bm{k}}^{z})\delta(\omega_{mn\bm{k}}-\omega)\,,$
(42c)

It can be seen that the coefficients $\eta^{xzx}$, $\sigma^{xzx}$, and
$\sigma^{zzz}$ are induced by the transitions only from the band 1 to 3 or
from the band 2 to 4, while $\sigma^{zxx}$ has no such limit. These
coefficients can be further simplified with the analytical expressions of all
these quantities, which can be obtained under the linear dispersion
approximation around the Dirac points, as shown in Appendix B.

Figure 2: Band structure (a) and JDOS (b) for AA-BG at $\Delta=0$ (dashed
curves) and $\Delta=0.4$ eV (solid curves).

Figure 2 (a) shows the band structure of AA-BG for $\Delta=0$ and $0.4$ eV.
With applying a gate voltage, the interlayer coupling shifts the energies of
the Dirac cones of each layer, while the electronic states at zero energy are
still degenerate. The bands 1 and 3 (or 2 and 4) are approximately parallel to
each other, and their energy differences are in the range of
$2\sqrt{\Delta^{2}+(\gamma_{1}+3\gamma_{3})^{2}}\leq\hbar\omega_{42\bm{k}}\leq
2\sqrt{\Delta^{2}+\gamma_{1}^{2}}\leq\hbar\omega_{31\bm{k}}\leq
2\sqrt{\Delta^{2}+(\gamma_{1}-3\gamma_{3})^{2}}$ due to $0\leq|g_{\bm{k}}|\leq
3$, where the middle value is obtained at the Dirac points and the other two
values are obtained at the M points. Figure 2 (b) gives the joint density of
states (JDOS) ${\cal J}_{31}(\omega)$ and ${\cal J}_{42}(\omega)$ for related
two pairs of bands, which are defined as

$\displaystyle{\cal J}_{nm}(\omega)$ $\displaystyle=\int
d\bm{k}\delta(\hbar\omega_{nm\bm{k}}-\hbar\omega)\,.$ (43)

These two JDOS are strongly localized in energy, regardless of whether there
is the gate voltage. For $\Delta=0.4$ eV, ${\cal D}_{42}(\omega)$ is nonzero
in the energy range of $[0.95,1.08]$ eV and ${\cal D}_{31}(\omega)$ is nonzero
in the energy range of $[1.08,1.21]$ eV.

###  III.2 Band structure of AB-BG

Figure 3: (a) Band structure and (b) JDOS for AB-BG at $\Delta=0$ (dashed
curves) and $\Delta=0.4$ eV (solid curves). The energetic locations of the
band gap ($E_{g}$) and the maximal values of ${\cal J}_{32}$, ${\cal J}_{42}$,
and ${\cal J}_{31}$ ($E_{1}$, $E_{2}$, $E_{3}$) are indicated. The inset:
$\bm{k}$-resolved energy difference $\hbar\omega_{32\bm{k}}$ for $\Delta=0.4$
eV. (c) $\Delta$ dependence of the $\bm{k}$ location for the minimum of
$\hbar\omega_{nm\bm{k}}$ along the K-M and K-$\Gamma$ directions. (d) $\Delta$
dependence of $E_{g}$, $E_{1}$, $E_{2}$, and $E_{3}$.

The Hamiltonian in Eq. (10) for AB-BG can be also analytically diagonalized,
as shown in Appendix C, but the expressions for the eigenenergies are too
complicated to provide meaningful physical insight, thus we discuss the band
structure based on numerical calculation. This work focuses on the electronic
transitions around the Dirac points, for convenience, the wavevectors are
expressed as
$\bm{k}=\bar{k}\frac{2\pi}{a_{0}}(\hat{\bm{x}}\cos\theta+\hat{\bm{y}}\sin\theta)+\bm{K}$
with $\theta=2n\pi/3$ along the K-M directions, and $\theta=(2n+1)\pi/3$ along
the K-$\Gamma$ directions. Figure 3 (a) gives the band structure for AB-BG at
gate voltages $\Delta=0$ and 0.4 eV. At $\Delta=0$, in each Dirac cone, the
two middle bands are degenerate at the Dirac points with $\bar{k}=0$ and other
three $\bm{k}$ points on the K-M paths with
$\bar{k}=-\frac{\gamma_{1}^{\prime}\gamma_{3}^{\prime}}{\sqrt{3}\pi{\gamma_{0}^{\prime}}^{2}}\sim
0.003$ (see details in Appendix C). Meanwhile, the energy differences,
$\hbar\omega_{31\bm{k}}$ and $\hbar\omega_{42\bm{k}}$, have minima at the
Dirac points. For nonzero gate voltage, the degeneracy at these points is
lifted. The eigenenergies at the Dirac points are
$\pm\Delta-\frac{\Delta^{\prime}}{2}$,
$\pm\sqrt{\Delta^{2}+\gamma_{1}^{2}}+\frac{\Delta^{\prime}}{2}$, and the
middle two bands around the Dirac points have the Mexican hat shape [35]. At
$\Delta=0.4$ eV, the energy difference $\hbar\omega_{32\bm{k}}$ shows a
minimum with increasing $\bar{k}$ for each $\theta$, as shown in the
$\bm{k}$-resolved energy difference in the inset, where the three-fold
rotational symmetry can be clearly seen around this Dirac point. Along the K-M
directions, the minima of $\hbar\omega_{32\bm{k}}$ appear around
$\bar{k}=0.027$ to give the band gap of $E_{g}=0.28$ eV; and along the
K-$\Gamma$ directions, the minima appear around $\bar{k}=0.023$, which have an
energy $E_{1}=0.4$ eV higher than the band gap and give a van Hove singularity
(VHS). Similar results can be found for $\hbar\omega_{42\bm{k}}$, and another
VHS appears with energy $E_{2}=0.97$ eV; however, $\hbar\omega_{31\bm{k}}$
shows a minimum at the Dirac points but no VHS appears. Figure 3 (b) gives
JDOS of ${\cal J}_{31}(\omega)$, ${\cal J}_{32}(\omega)$, ${\cal
J}_{41}(\omega)$, and ${\cal J}_{42}(\omega)$ at $\Delta=0$ and $0.4$ eV. The
gate voltage changes these JDOS significantly around the band edge. ${\cal
J}_{32}(\omega)$ and ${\cal J}_{42}(\omega)$ have divergences at the VHS
points with energies $E_{1}$ and $E_{2}$, respectively; and ${\cal
J}_{31}(\omega)$ has a peak located at $E_{3}\sim 0.97$ eV around the band
edge, which is induced by the nearly parallel bands (1, 3) around the Dirac
points.

The VHS points do not appear for all gate voltages. Figure 3 (c) exhibits
$\Delta$ dependence of the $\bar{k}$ value for the minimal energy of
$\hbar\omega_{32\bm{k}}$ and $\hbar\omega_{42\bm{k}}$ for $\theta$ along the
K-M and K-$\Gamma$ directions, respectively. Along the K-M directions,
$\hbar\omega_{32\bm{k}}$ has a minimum value at nonzero $\bar{k}$ for all
$\Delta$, which gives the band gap $E_{g}$ of the system; while along the
K-$\Gamma$ directions, the minimum energy $E_{1}$ moves to a nonzero $\bar{k}$
only for $\Delta\geq 0.023$ eV, where VHS appears as well. Note that the JDOS
${\cal J}_{32\bm{k}}$ shows a maximum at the band edge when there is no VHS
for $\Delta<0.023$ eV. However, the minima of $\hbar\omega_{42\bm{k}}$ along
the K-M and K-$\Gamma$ directions locate not at the Dirac points only for
$\Delta\geq 0.174$ eV, where VHS appears as well. For $\Delta<0.174$ eV,
${\cal J}_{42}(\omega)$ also shows a maximum at the band edge between the
bands 4 and 2, where this energy is still noted as $E_{2}$; the maximum of
${\cal J}_{31}(\omega)$ also locates at the band edge between bands 3 and 1,
where this energy is still noted as $E_{3}$. The gate voltage dependences of
these energies $E_{g}$, $E_{1}$, $E_{2}$, and $E_{3}$ are shown in Fig. 3 (d).

###  III.3 Injection coefficients and shift conductivities at $\Delta=0.4$ eV

In this section we present the numerical results for injection coefficient
$\eta^{xzx}(\omega)$ and shift conductivities $\sigma^{yyy}(\omega)$,
$\sigma^{xzx}(\omega)$, $\sigma^{zxx}(\omega)$, and $\sigma^{zzz}(\omega)$.
The parameters are chosen as $T=300$ K, $\mu=$ 0, $\Delta=0.4$ eV. During the
numerical calculation, the Brillouin zone is divided into a $3000\times 3000$
homogeneous grid. The $\delta$ functions in Eqs. (23) and (25) are
approximated by a Gaussian function as
$\delta(\omega)=\frac{\hbar}{\sqrt{\pi}\Gamma}e^{-(\hbar\omega)^{2}/\Gamma^{2}}$
with the Gaussian broadening $\Gamma=10$ meV.

Figure 4: Injection coefficients and shift conductivities for AA-BG and AB-BG
at $\Delta=0.4$ eV. (a) $\eta^{xzx}$ for AA-BG and AB-BG, (b) $\sigma^{yyy}$
for AB-BG, (c) $\sigma^{xzx}$, $\sigma^{zxx}$, and $\sigma^{zzz}$ for AA-BG,
and (d) $\sigma^{xzx}$, $\sigma^{zxx}$, and $\sigma^{zzz}$ for AB-BG.

Figure 4 (a) shows the injection coefficient spectra for AA-BG and AB-BG. For
the injection in AA-BG, the spectrum is just a peak located in a very narrow
energy range $1.069$ eV$<\hbar\omega<1.087$ eV with an absolute value about
$0.067$ $\rm A\cdot s^{-1}\cdot m/V^{2}$. From the analytic results shown in
Eq. (64), the spectra include two contributions at different photon energy
regions: one is from the optical transition between the bands (1, 3) for
photon energy $\hbar\omega>2\sqrt{\Delta^{2}+\gamma_{1}^{2}}$ or $1.078$
eV$<\hbar\omega<$1.087 eV, and the other is between the bands (2, 4) for
$\hbar\omega<2\sqrt{\Delta^{2}+\gamma_{1}^{2}}$ or $1.069$
eV$<\hbar\omega<$1.078 eV; both magnitudes are nearly proportional to
$\hbar\omega-2\sqrt{\Delta^{2}+\gamma_{1}^{2}}$. These two contributions merge
as a single peak just because the $\delta$ function is numerically broadened
with $\Gamma=10$ meV, which is even larger than each energy region. The
injection coefficient $\eta^{xzx}$ in AB-BG starts with photon energy higher
than the gap, i.e., $\hbar\omega>0.28$ eV, and reaches its maximum value of
$25~{}\rm A\cdot s^{-1}\cdot m/V^{2}$ in amplitude at $\hbar\omega=0.45$ eV,
which is slightly larger than the first VHS energy of JDOS $E_{1}$; the energy
difference arises from the zero electron velocity at this VHS. Considering the
thickness of a bilayer graphene as $2c=6.7$ Å, the effective bulk injection
coefficient is $3.7\times 10^{10}~{}\rm\mu A\cdot s^{-1}V^{-2}$, which is
nearly 50 times larger that that in bulk GaAs [36]. After this peak, the
amplitude of injection coefficient decreases as the photon energy increases,
except for a small peak located around the JDOS peak at higher energy $E_{2}$
or $E_{3}$. It can be seen that the injection coefficient for AB-BG is about
two orders of magnitude larger than that for AA-BG. To have a direct
impression on these values, we give an estimation on how large the injection
current can be in AB-BG. Based on the Eq. (26), when the laser is a
$45^{\circ}$ obliquely incident $p$-polarized light with photon energy of 0.45
eV, light intensity of $I=0.1~{}\rm GW/cm^{2}$, and pulse duration of $\tau=1$
ps, the generated injection current is
$2\eta^{xzx}\frac{I}{2c\epsilon_{0}}W\tau\sim 9~{}\rm mA$ for an electrode
with a width $W=1~{}\rm\mu m$.

Then we turn to the shift conductivities, as shown in Figs. 4 (b–d). Figure 4
(c) gives the shift conductivity for AA-BG. It can be seen that the component
$\sigma^{zzz}$ is about one order of magnitude larger than $\sigma^{xzx}$, or
is at least two order of magnitude larger than $\sigma^{zxx}$. Both
$\sigma^{zzz}$ and $\sigma^{xzx}$ have nonzero values only in the very narrow
energy regions, similar to the injection coefficient. These results are
consistent with the analytic results shown in Eqs. (65–66). Interestingly,
$\sigma^{xzx}$ includes the contributions from the band 1 to 3 and from the
band 2 to 4 but with opposite signs. For AB-BG shown in Figs. 4 (b) and (d),
all nonzero components start from the band edge $\hbar\omega\geq E_{g}$.
Different from the injection coefficients, the shift conductivities at the
band edge are nonzero, and show prominent peaks. Especially, $\sigma^{yyy}$
shows a large value about $6\times 10^{-13}~{}\rm A\cdot m/V^{2}$ at the band
edge and it drops quickly with increasing the photon energy. The effective
bulk shift conductivity is $896~{}\rm\mu A/V^{2}$, which is several times
larger than in GeSe ($200~{}\rm\mu A/V^{2}$) [36]. Besides, the component
$\sigma^{zzz}$ is at least one order of magnitude smaller than other nonzero
components, totally different from the case of AA-BG, where it is the largest
one. The spectra of $\sigma^{xzx}$ and $\sigma^{zxx}$ have similar amplitude
around a few $10^{-14}~{}\rm A\cdot m/V^{2}$, which is a few tens of times
smaller than the peak of $\sigma^{yyy}$; they also show some fine structures
around those characteristic energies $E_{1}$, $E_{2}$, and $E_{3}$. We repeat
the above estimation for the shift current using the same parameters but
$\hbar\omega=0.3$ eV, and then obtain the generated shift current of
$2\sigma^{yyy}\frac{I}{2c\epsilon_{0}}W\sim 0.23$ mA.

###  III.4 Effects of Gate voltage

Figure 5: Gate voltage dependence of injection coefficients and shift
conductivities at zero chemical potential. (a) $\eta^{xzx}(\omega)$ and (b)
$\sigma^{zzz}(\omega)$ for AA-BG, (c) $\eta^{xzx}(\omega)$, (d)
$\sigma^{xzx}(\omega)$, (e) $\sigma^{yyy}(\omega)$, and (f)
$\sigma^{zxx}(\omega)$ for AB-BG. In (c–g), dashed curves indicate the
characteristic energies $E_{g}$, $E_{1}$, $E_{2}$, and $E_{3}$ for AB-BG.

Figure 5 gives the gate voltage dependence of the injection coefficients and
shift conductivities for AA-BG and AB-BG at zero chemical potential. Note that
the negative gate voltage leads to opposite coefficients, which are consistent
with the results by Xiong et al. [30], thus only positive gate voltages are
shown here.

Figures 5 (a) and (b) show the spectra of $\eta^{xzx}$ and $\sigma^{zzz}$ for
AA-BG, respectively. As indicated in previous section, both spectra for
different gate voltages are nonzero in a very narrow photon energy region.
With the increase of the gate voltage, the region moves to larger energy and
the values of both spectra increase, which are indicated by $\propto\Delta$ in
Eqs. (68) and (70). Figure 5 (c) gives the injection coefficient $\eta^{xzx}$
for AB-BG. At each gate voltage, the injection coefficient shows two peaks
located at photon energies slightly larger than $E_{1}$ and $E_{2}$, which
have been discussed in previous section. As the gate voltage $\Delta$ varies,
the peak amplitude reaches a maximum at $\Delta\sim 0.2$ eV. The shift
conductivities $\sigma^{xzx},\sigma^{yyy}$ and $\sigma^{zxx}$ for AB-BG are
plotted in Figs. 5 (d–f). They show some similar characteristics: (1) The
spectra are located at about the band gap similar to the case of $\Delta=0.4$
eV, and their amplitudes increase with the decrease of $\Delta$;
$\sigma^{xzx}$ and $\sigma^{zxx}$ increase much faster than $\sigma^{yyy}$.
(2) There exist sign changes of shift conductivities.

###  III.5 Effects of Chemical potential

Figure 6: Chemical potential dependence of injection coefficients and and
shift conductivities at $\Delta=0.4$ eV. (a) $\eta^{xzx}(\omega)$ and (b)
$\sigma^{zzz}(\omega)$ for AA-BG, (c) $\eta^{xzx}(\omega)$, (d)
$\sigma^{xzx}(\omega)$, (e) $\sigma^{yyy}(\omega)$, and (f)
$\sigma^{zxx}(\omega)$ for AB-BG. The dashed lines in (c–f) indicate the
position of the conduction and valence band edges for AB-BG.

The chemical potential $\mu$ dependence of injection coefficients and shift
conductivities at $\Delta=0.4$ eV are depicted in Fig. 6 with the same layout
as Fig. 5. For AA-BG in Figs. 6 (a) and (b), they show very similar asymmetric
dependence on the chemical potential: with the increase of the chemical
potential, the values of all coefficients increase and the locations shift to
higher or lower photon energies depending on the sign of the chemical
potential. For positive chemical potential, the transitions between bands (1,
3) are suppressed according to the Pauli blocking effects, while new extra
transitions between bands (2, 4) appear due to the additional free electrons
in the band 2. The extra transitions require lower photon energy and red shift
the spectra, and they also correspond to larger JDOS, leading to larger
coefficients. Similar results can be analyzed for negative chemical potential,
but with switching the band pairs (1, 2) and (3, 4).

In AB-BG, the chemical potential $\mu$ has different effects, as shown in
Figs. 6 (c–f). Due to the existence of the band gap, the spectra are hardly
changed when the chemical potential lies in the gap. When $\mu$ is above the
conduction band edge or below the valence band edge, the main peak of
$\eta^{xzx}$ around 0.5 eV is reduced gradually due to the Pauli blocking, and
there appear new transitions between the bands (1, 2) or (3, 4) to give
additional injections with opposite signs. Similar results are obtained for
the shift conductivities.

##  IV Conclusion

In this paper we have studied the gate voltage induced injection current and
shift current in AA- and AB-stacked bilayer graphene. The gate voltage plays a
crucial role in breaking the inversion symmetry of bilayer graphene to induce
photogalvanic effects, and at the same time it effectively changes the band
structure for AB-BG with opening gaps located in the K-M directions and
inducing additional VHS located in the K-$\Gamma$ directions. In AA-BG, the
injection and shift currents are mainly induced by optical transitions between
two pairs of nearly parallel bands; the coefficient spectra locate in a very
narrow photon energy region of about 20 meV. In AB-BG, the optical transition
can occur between any possible band pairs, and the structure of spectra are
strongly determined by the band gap and the VHS energies. For both structures,
the injection and shift currents can be generated by the existence of an
oblique $p$-polarized light component, while the in-plane shift currents in
AB-BG can also be generated by normal incident lights. The out-of-plane shift
current finally results in a static electric polarization between layers. The
stacking order has significant effects on both currents. The injection
coefficient for AA-BG is about two orders of magnitude smaller than that for
AB-BG, while the shift conductivities are mostly in the same order of
magnitude. All these coefficients can be effectively modulated by the gate
voltage and the chemical potential. Our results suggest that gate voltage
controlled bilayer graphene can be used to realize tunable optoelectronic
detectors working in the mid-infrared.

###### Acknowledgements.

This work has been supported by National Natural Science Foundation of China
Grant No. 12034003, 12004379, and 62250065. J.L.C. acknowledges the support
from Talent Program of CIOMP.

##  Appendix A Berry connections of AA-BG

The general expression for the Berry connection of AA-BG is

$\displaystyle\bm{\xi}_{nm\bm{k}}$
$\displaystyle=\left(\sqrt{1-\alpha_{n}{\cal
N}_{\beta_{n}\bm{k}}}\sqrt{1-\alpha_{m}{\cal
N}_{\beta_{m}\bm{k}}}+\alpha_{n}\alpha_{m}\sqrt{1+\alpha_{n}{\cal
N}_{\beta_{n}\bm{k}}}\sqrt{1+\alpha_{m}{\cal
N}_{\beta_{m}\bm{k}}}\right)\times$
$\displaystyle\times\frac{1+\beta_{n}\beta_{m}}{8}\left[\hat{g}_{\bm{k}}^{\ast}(i\bm{\nabla}_{\bm{k}}\hat{g}_{\bm{k}})+\hat{\bm{y}}d\right]$
$\displaystyle+\left(\alpha_{m}\sqrt{1-\alpha_{n}{\cal
N}_{\beta_{n}\bm{k}}}\sqrt{1+\alpha_{m}{\cal
N}_{\beta_{m}\bm{k}}}+\alpha_{n}\sqrt{1+\alpha_{n}{\cal
N}_{\beta_{n}\bm{k}}}\sqrt{1-\alpha_{m}{\cal N}_{\beta_{m}\bm{k}}}\right)$
$\displaystyle\times\frac{\beta_{n}\beta_{m}-1}{8}\left[\hat{g}_{\bm{k}}^{\ast}(i\bm{\nabla}_{\bm{k}}\hat{g}_{\bm{k}})-\hat{\bm{y}}d\right]$
$\displaystyle+\frac{i\delta_{\beta_{n}\beta_{m}}}{2}\left(\sqrt{1-\alpha_{n}{\cal
N}_{\beta_{n}\bm{k}}}\bm{\nabla}_{\bm{k}}\sqrt{1-\alpha_{m}{\cal
N}_{\beta_{m}\bm{k}}}+\alpha_{n}\alpha_{m}\sqrt{1+\alpha_{n}{\cal
N}_{\beta_{n}\bm{k}}}\bm{\nabla}_{\bm{k}}\sqrt{1+\alpha_{m}{\cal
N}_{\beta_{m}\bm{k}}}\right)$ $\displaystyle+\left(\sqrt{1-\alpha_{n}{\cal
N}_{\beta_{n}\bm{k}}}+\alpha_{n}\sqrt{1+\alpha_{n}{\cal
N}_{\beta_{n}\bm{k}}}\right)\left(\sqrt{1-\alpha_{m}{\cal
N}_{\beta_{m}\bm{k}}}+\alpha_{m}\sqrt{1+\alpha_{m}{\cal
N}_{\beta_{m}\bm{k}}}\right)$
$\displaystyle\times\frac{1+\beta_{n}\beta_{m}}{8}c\hat{\bm{z}}$ (44)

with $d=\sqrt{3}/3a_{0}$. Here we give the $x$-component between different
bands as

$\displaystyle r^{x}_{13\bm{k}}$
$\displaystyle=-r^{x}_{31\bm{k}}=-\frac{i}{2}\frac{\frac{\partial{\cal
N}_{-1\bm{k}}}{\partial k_{x}}}{\sqrt{1-{\cal
N}_{-1\bm{k}}^{2}}}=-\frac{i}{2}\frac{\gamma_{3}}{|\Delta|}(1-{\cal
N}_{-1\bm{k}}^{2})\frac{\partial|g_{\bm{k}}|}{\partial k_{x}}\,,$ (45a)
$\displaystyle r^{x}_{24\bm{k}}$
$\displaystyle=-r^{x}_{42\bm{k}}=-\frac{i}{2}\frac{\frac{\partial{\cal
N}_{+1\bm{k}}}{\partial k_{x}}}{\sqrt{1-{\cal
N}_{+1\bm{k}}^{2}}}=-\frac{i}{2}\frac{\gamma_{3}}{|\Delta|}(1-{\cal
N}_{+1\bm{k}}^{2})\frac{\partial|g_{\bm{k}}|}{\partial k_{x}}\,,$ (45b)
$\displaystyle r^{x}_{12\bm{k}}$
$\displaystyle=r^{x}_{21\bm{k}}=-r^{x}_{34\bm{k}}=-r^{x}_{43\bm{k}}$
$\displaystyle=\frac{1}{4}\left[\sqrt{1+{\cal N}_{-1\bm{k}}}\sqrt{1-{\cal
N}_{+1\bm{k}}}+\sqrt{1-{\cal N}_{-1\bm{k}}}\sqrt{1+{\cal
N}_{+1\bm{k}}}\right]\left[\hat{g}_{\bm{k}}^{\ast}\left(i\frac{\partial\hat{g}_{\bm{k}}}{\partial
k_{x}}\right)\right]\,,$ (45c) $\displaystyle r^{x}_{32\bm{k}}$
$\displaystyle=-r^{x}_{23\bm{k}}=r^{x}_{14\bm{k}}=-r^{x}_{41\bm{k}}$
$\displaystyle=\frac{1}{4}\left[\sqrt{1+{\cal N}_{-1\bm{k}}}\sqrt{1-{\cal
N}_{+1\bm{k}}}-\sqrt{1-{\cal N}_{-1\bm{k}}}\sqrt{1+{\cal
N}_{+1\bm{k}}}\right]\left[\hat{g}_{\bm{k}}^{\ast}\left(i\frac{\partial\hat{g}_{\bm{k}}}{\partial
k_{x}}\right)\right]\,.$ (45d)

Combining with other quantities in Eqs. (39) and (41), the injection
coefficients and the shift conductivities can be evaluated. For the latter
use, we also need

$\displaystyle{\cal V}_{21\bm{k}}$
$\displaystyle=\frac{2\gamma_{3}}{\hbar}{\cal
N}_{-1\bm{k}}\frac{\partial|g_{\bm{k}}|}{\partial k_{x}}\,,$ (46)
$\displaystyle{\cal V}_{43\bm{k}}$
$\displaystyle=\frac{2\gamma_{3}}{\hbar}{\cal
N}_{+1\bm{k}}\frac{\partial|g_{\bm{k}}|}{\partial k_{x}}\,.$ (47)

##  Appendix B Analytical expressions of $\eta^{xzx}$, $\sigma^{xzx}$, and
$\sigma^{zzz}$ in AA-BG under the linear dispersion approximation

Here we give the analytic results for $\eta^{xzx}$ in Eq. (40), $\sigma^{xzx}$
in Eq. (42a) and $\sigma^{zzz}$ in Eq. (42b) under the linear dispersion
approximation around the Dirac points. The term of $\sigma^{zxx}$ is not
discussed due to its very small magnitude, as shown in Fig. 4 (c).

The integrands of $\eta^{xzx}$, $\sigma^{xzx}$, and $\sigma^{zzz}$ are
functions of $|g_{\bm{k}}|$, $\frac{\partial|g_{\bm{k}}|}{\partial k_{x}}$,
and $\frac{\partial^{2}|g_{\bm{k}}|}{\partial k_{x}^{2}}$, where all terms
involving $|g_{\bm{k}}|$ can be simplified by using the properties of the
$\delta$ function. The function $\delta(\hbar\omega_{nm\bm{k}}-\hbar\omega)$
is nonzero only for $|g_{\bm{k}}|=G_{nm}$ with

$\displaystyle\gamma_{3}G_{31}$
$\displaystyle=\gamma_{1}-\sqrt{\left(\frac{\hbar\omega}{2}\right)^{2}-\Delta^{2}}\,,\text{
for }\hbar\omega\geq 2\sqrt{\Delta^{2}+\gamma_{1}^{2}}\,,$ (48)
$\displaystyle\gamma_{3}G_{42}$
$\displaystyle=\sqrt{\left(\frac{\hbar\omega}{2}\right)^{2}-\Delta^{2}}-\gamma_{1}\,,\text{
for }\hbar\omega\leq 2\sqrt{\Delta^{2}+\gamma_{1}^{2}}\,.$ (49)

Further we get

$\displaystyle\left.\left({\cal
N}_{-1\bm{k}}\right)\right|_{|g_{\bm{k}}|=G_{31}}=-\left.\left({\cal
N}_{+1\bm{k}}\right)\right|_{|g_{\bm{k}}|=G_{42}}=-\sqrt{1-\left(\frac{2\Delta}{\hbar\omega}\right)^{2}}\,.$
(50)

  1. 1.

By substituting the expressions of ${\cal V}_{nm\bm{k}}^{x}$,
$r_{31\bm{k}}^{x}$, $r_{13\bm{k}}^{z}$, $r_{42\bm{k}}^{x}$, and
$r_{24\bm{k}}^{z}$, $\eta^{xzx}$ becomes

$\displaystyle\eta^{xzx}=$ $\displaystyle\frac{e^{3}}{2\pi\hbar^{2}}\int
d\bm{k}\left(\frac{c\gamma_{3}^{2}}{2\hbar|\Delta|}\right)\left\\{f_{12\bm{k}}{\cal
N}_{-1\bm{k}}^{2}(1-{\cal
N}_{-1\bm{k}}^{2})\left(\frac{\partial|g_{\bm{k}}|}{\partial
k_{x}}\right)^{2}\delta(\omega_{31\bm{k}}-\omega)\right.$
$\displaystyle\left.+f_{34\bm{k}}{\cal N}_{+1\bm{k}}^{2}(1-{\cal
N}_{+1\bm{k}}^{2})\left(\frac{\partial|g_{\bm{k}}|}{\partial
k_{x}}\right)^{2}\delta(\omega_{42\bm{k}}-\omega)\right\\}$ $\displaystyle=$
$\displaystyle\frac{e^{3}c|\Delta|}{\pi\hbar^{2}(\hbar\omega)^{2}}\left[1-\left(\frac{2\Delta}{\hbar\omega}\right)^{2}\right]\left\\{f_{13\bm{k}}|_{|g_{\bm{k}}|=G_{31}}{\cal
F}_{31}(\omega)+f_{24\bm{k}}|_{|g_{\bm{k}}|=G_{42}}{\cal
F}_{42}(\omega)\right\\}\,,$ (51)

with

$\displaystyle{\cal F}_{nm}(\omega)$ $\displaystyle=\int
d\bm{k}\left(\gamma_{3}\frac{\partial|g_{\bm{k}}|}{\partial
k_{x}}\right)^{2}\delta(\hbar\omega_{nm\bm{k}}-\hbar\omega)\,.$ (52)

  2. 2.

To get the result for $\sigma^{xzx}$, we use

$\displaystyle\frac{\partial{\cal N}_{-1\bm{k}}}{\partial k_{x}}=(1-{\cal
N}_{-1\bm{k}}^{2})^{3/2}\frac{\gamma_{3}}{|\Delta|}\frac{\partial|g_{\bm{k}}|}{\partial
k_{x}}$ (53)

to get

$\displaystyle r_{31\bm{k}}^{z}\frac{\partial r_{13\bm{k}}^{x}}{\partial
k_{x}}+r_{31\bm{k}}^{x}\frac{\partial r_{13\bm{k}}^{z}}{\partial k_{x}}=$
$\displaystyle\frac{ic}{4}(1+{\cal N}_{-1\bm{k}}^{2})(1-{\cal
N}_{-1\bm{k}}^{2})^{3/2}\left(\frac{\gamma_{3}}{|\Delta|}\frac{\partial|g_{\bm{k}}|}{\partial
k_{x}}\right)^{2}$ $\displaystyle-\frac{ic}{4}{\cal N}_{-1\bm{k}}(1-{\cal
N}_{-1\bm{k}}^{2})\frac{\gamma_{3}}{|\Delta|}\frac{\partial^{2}|g_{\bm{k}}|}{\partial
k_{x}^{2}}\,.$ (54)

Similar expressions can be obtained for terms involving $\bm{r}_{32\bm{k}}$.
Then we get

$\displaystyle\sigma^{xzx}=$
$\displaystyle\frac{e^{3}c}{4\pi\hbar(\hbar\omega)^{2}}\left\\{\left[2-\left(\frac{2\Delta}{\hbar\omega}\right)^{2}\right]\frac{2|\Delta|}{\hbar\omega}\left[f_{13\bm{k}}|_{|g_{\bm{k}}|=G_{31}}{\cal
F}_{31}(\omega)+f_{24\bm{k}}|_{|g_{\bm{k}}|=G_{42}}{\cal
F}_{42}(\omega)\right]\right.$
$\displaystyle\left.-|\Delta|\sqrt{1-\left(\frac{2\Delta}{\hbar\omega}\right)^{2}}\left[f_{13\bm{k}}|_{|g_{\bm{k}}|=G_{31}}{\cal
Q}_{31}(\omega)-f_{24\bm{k}}|_{|g_{\bm{k}}|=G_{42}}{\cal
Q}_{42}(\omega)\right]\right\\}\,.$ (55)

with

$\displaystyle{\cal Q}_{nm}(\omega)$ $\displaystyle=\int
d\bm{k}\gamma_{3}\frac{\partial^{2}|g_{\bm{k}}|}{\partial
k^{2}_{x}}\delta(\hbar\omega_{nm\bm{k}}-\hbar\omega)\,.$ (56)

  3. 3.

The term of $\sigma^{zzz}(\omega)$ becomes

$\displaystyle\sigma^{zzz}(\omega)=$
$\displaystyle\frac{e^{3}}{2\pi\hbar^{2}}\int
d\bm{k}\left\\{f_{13\bm{k}}\frac{c^{2}}{4}{\cal
N}_{-1\bm{k}}^{2}c\sqrt{1-{\cal
N}^{2}_{-1\bm{k}}}\delta(\omega_{31\bm{k}}-\omega)\right.$
$\displaystyle\left.+f_{24\bm{k}}\frac{c^{2}}{4}{\cal
N}_{+1\bm{k}}^{2}c\sqrt{1-{\cal
N}^{2}_{+1\bm{k}}}\delta(\omega_{42\bm{k}}-\omega)\right\\}$ $\displaystyle=$
$\displaystyle\frac{e^{3}c^{3}}{4\pi\hbar}\frac{|\Delta|}{\hbar\omega}\left[1-\left(\frac{2\Delta}{\hbar\omega}\right)^{2}\right]\left[f_{13\bm{k}}|_{|g_{\bm{k}}|=G_{31}}{\cal
J}_{31}(\omega)+f_{24\bm{k}}|_{|g_{\bm{k}}|=G_{42}}{\cal
J}_{42}(\omega)\right]\,.$ (57)

with

$\displaystyle{\cal J}_{nm}(\omega)$ $\displaystyle=\int
d\bm{k}\delta(\hbar\omega_{nm\bm{k}}-\hbar\omega)\,.$ (58)

When the optical transition occurs just around the Dirac points $\bm{K}$, we
can approximate $|g_{\bm{k}+\bm{K}}|=\sqrt{3}a_{0}k/2$, then the $\delta$
functions can be worked out as

$\displaystyle\delta(2\sqrt{\Delta^{2}+(\gamma_{3}|g_{\bm{k}}|-\gamma_{1})^{2}}-\hbar\omega)$
$\displaystyle=\frac{\delta\left(k-{2G_{31}}/(\sqrt{3}a_{0})\right)}{\sqrt{3}a_{0}|\gamma_{3}|\sqrt{1-\left(\frac{2\Delta}{\hbar\omega}\right)^{2}}}\theta(\hbar\omega-2\sqrt{\Delta^{2}+\gamma_{1}^{2}})\,,$
(59)
$\displaystyle\delta(2\sqrt{\Delta^{2}+(\gamma_{3}|g_{\bm{k}}|+\gamma_{1})^{2}}-\hbar\omega)$
$\displaystyle=\frac{\delta\left(k-{2G_{42}}/(\sqrt{3}a_{0})\right)}{\sqrt{3}a_{0}|\gamma_{3}|\sqrt{1-\left(\frac{2\Delta}{\hbar\omega}\right)^{2}}}\theta(2\sqrt{\Delta^{2}+\gamma_{1}^{2}}-\hbar\omega)\,.$
(60)

Then we get

$\displaystyle\begin{pmatrix}{\cal J}_{31}(\omega)\\\ {\cal
J}_{42}(\omega)\end{pmatrix}$
$\displaystyle=\frac{8\pi}{3a_{0}^{2}\gamma_{3}^{2}{\sqrt{1-\left(\frac{2\Delta}{\hbar\omega}\right)^{2}}}}\left|\gamma_{1}-\sqrt{\left(\frac{\hbar\omega}{2}\right)^{2}-\Delta^{2}}\right|\begin{pmatrix}\theta(\hbar\omega-2\sqrt{\Delta^{2}+\gamma_{1}^{2}})\\\
\theta(2\sqrt{\Delta^{2}+\gamma_{1}^{2}}-\hbar\omega)\end{pmatrix}\,,$ (61)
$\displaystyle\begin{pmatrix}{\cal F}_{31}(\omega)\\\ {\cal
F}_{42}(\omega)\end{pmatrix}$
$\displaystyle=\frac{3a_{0}^{2}\gamma_{3}^{2}}{8}\begin{pmatrix}{\cal
J}_{31}(\omega)\\\ {\cal J}_{42}(\omega)\end{pmatrix}\,,$ (62)
$\displaystyle\begin{pmatrix}{\cal Q}_{31}(\omega)\\\ {\cal
Q}_{42}(\omega)\end{pmatrix}$
$\displaystyle=-\frac{\pi}{{\sqrt{1-\left(\frac{2\Delta}{\hbar\omega}\right)^{2}}}}\begin{pmatrix}\theta(\hbar\omega-2\sqrt{\Delta^{2}+\gamma_{1}^{2}})\\\
\theta(2\sqrt{\Delta^{2}+\gamma_{1}^{2}}-\hbar\omega)\end{pmatrix}\,,$ (63)

where two Dirac points have been counted in the integration. In such
approximation, the expressions for $\eta^{xzx}$, $\sigma^{xzx}$, and
$\sigma^{zzz}$ are expressed as

$\displaystyle\eta^{xzx}(\omega)=$
$\displaystyle\frac{e^{3}c|\Delta|\sqrt{1-\left(\frac{2\Delta}{\hbar\omega}\right)^{2}}}{\hbar^{2}(\hbar\omega)^{2}}\left|\gamma_{1}-\sqrt{\left(\frac{\hbar\omega}{2}\right)^{2}-\Delta^{2}}\right|\left({\cal
M}_{31}(\omega)+{\cal M}_{42}(\omega)\right)\,,$ (64)
$\displaystyle\sigma^{xzx}(\omega)=$
$\displaystyle\frac{e^{3}c|\Delta|\left(\hbar^{2}\omega^{2}-2\Delta^{2}\right)}{2\hbar(\hbar\omega)^{4}\sqrt{1-\left(\frac{2\Delta}{\hbar\omega}\right)^{2}}}\left|\sqrt{1-\left(\frac{2\Delta}{\hbar\omega}\right)^{2}}-\frac{2\gamma_{1}}{\hbar\omega}\right|\left({\cal
M}_{31}(\omega)+{\cal M}_{42}(\omega)\right)$
$\displaystyle-\frac{ce^{3}|\Delta|}{4\hbar(\hbar\omega)^{2}}\left({\cal
M}_{31}(\omega)-{\cal M}_{42}(\omega)\right)\,,$ (65)
$\displaystyle\sigma^{zzz}(\omega)=$
$\displaystyle\frac{e^{3}c^{3}|\Delta|\sqrt{1-\left(\frac{2\Delta}{\hbar\omega}\right)^{2}}}{3\hbar(a_{0}\gamma_{3})^{2}}\left|\sqrt{1-\left(\frac{2\Delta}{\hbar\omega}\right)^{2}}-\frac{2\gamma_{1}}{\hbar\omega}\right|\left({\cal
M}_{31}(\omega)+{\cal M}_{42}(\omega)\right)\,,$ (66)

respectively, with

$\displaystyle\begin{pmatrix}{\cal M}_{31}(\omega)\\\ {\cal
M}_{42}(\omega)\end{pmatrix}=\begin{pmatrix}f_{13\bm{k}}|_{|g_{\bm{k}}|=G_{31}}\theta(\hbar\omega-2\sqrt{\Delta^{2}+\gamma_{1}^{2}})\\\
f_{24\bm{k}}|_{|g_{\bm{k}}|=G_{42}}\theta(2\sqrt{\Delta^{2}+\gamma_{1}^{2}}-\hbar\omega)\end{pmatrix}\,.$
(67)

Through the Taylor expansion, the above expressions around frequency
$2\sqrt{\Delta^{2}+\gamma_{1}^{2}}$ can be approximated as

$\displaystyle\eta^{xzx}(\omega)\approx$
$\displaystyle\frac{ce^{3}|\Delta|\left|2\sqrt{\gamma_{1}^{2}+\Delta^{2}}-\hbar\omega\right|}{8\hbar^{2}(\gamma_{1}^{2}+\Delta^{2})}\left({\cal
M}_{31}(\omega)+{\cal M}_{42}(\omega)\right)\,,$ (68)
$\displaystyle\sigma^{xzx}(\omega)\approx$
$\displaystyle\frac{ce^{3}|\Delta|(2\gamma_{1}^{2}+\Delta^{2})\left|2\sqrt{\gamma_{1}^{2}+\Delta^{2}}-\hbar\omega\right|}{32\hbar\gamma_{1}^{2}\sqrt{\gamma_{1}^{2}+\Delta^{2}}^{3}}\left({\cal
M}_{31}(\omega)+{\cal M}_{42}(\omega)\right)$
$\displaystyle-\frac{ce^{3}|\Delta|}{16\hbar(\gamma_{1}^{2}+\Delta^{2})}\left({\cal
M}_{31}(\omega)-{\cal M}_{42}(\omega)\right)\,,$ (69)
$\displaystyle\sigma^{zzz}(\omega)\approx$
$\displaystyle\frac{ce^{3}|\Delta|\left|2\sqrt{\gamma_{1}^{2}+\Delta^{2}}-\hbar\omega\right|}{6\hbar
a_{0}^{2}\gamma_{3}^{2}(\gamma_{1}^{2}+\Delta^{2})}\left({\cal
M}_{31}(\omega)+{\cal M}_{42}(\omega)\right)\,.$ (70)

##  Appendix C Eigenenergies of AB-BG

The eigenenergies $\epsilon$ satisfy the equation

$\displaystyle|H^{\text{AB}}_{\bm{k}}-\epsilon|=0\,,$ (71)

or

$\displaystyle\epsilon^{4}+x_{2}\epsilon^{2}+x_{1}\epsilon+x_{0}=0\,,$ (72)

with

$\displaystyle x_{2}=$
$\displaystyle-{\gamma_{1}^{\prime}}^{2}-\left(2{\gamma_{0}^{\prime}}^{2}+{\gamma_{3}^{\prime}}^{2}+2{\gamma_{4}^{\prime}}^{2}\right)|g_{\bm{k}}|^{2}-2\left[\Delta^{2}+\left(\frac{\Delta^{\prime}}{2}\right)^{2}\right]\,,$
(73) $\displaystyle x_{1}=$
$\displaystyle-4{\gamma_{0}^{\prime}}{\gamma_{4}^{\prime}}\left({\gamma_{1}^{\prime}}|g_{\bm{k}}|^{2}+{\gamma_{3}^{\prime}}\text{Re}\left[g_{\bm{k}}^{3}\right]\right)+\Delta^{\prime}\left({\gamma_{3}^{\prime}}^{2}|g_{\bm{k}}|^{2}-{\gamma_{1}^{\prime}}^{2}\right)\,,$
(74) $\displaystyle x_{0}=$
$\displaystyle\left({\gamma_{0}^{\prime}}^{2}-{\gamma_{4}^{\prime}}^{2}\right)^{2}|g_{\bm{k}}|^{4}-2{\gamma_{3}^{\prime}}\left[{\gamma_{1}^{\prime}}\left({\gamma_{0}^{\prime}}^{2}+{\gamma_{4}^{\prime}}^{2}\right)-{\gamma_{0}^{\prime}}{\gamma_{4}^{\prime}}\Delta^{\prime}\right]\text{Re}[g_{\bm{k}}^{3}]$
$\displaystyle+\left\\{{\gamma_{3}^{\prime}}^{2}\left[{\gamma_{1}^{\prime}}^{2}+\Delta^{2}-\left(\frac{\Delta^{\prime}}{2}\right)^{2}\right]-\left(2{\gamma_{0}^{\prime}}^{2}-{\gamma_{3}^{\prime}}^{2}\right)\left[\Delta^{2}-\left(\frac{\Delta^{\prime}}{2}\right)^{2}\right]-2{\gamma_{0}^{\prime}}{\gamma_{1}^{\prime}}{\gamma_{4}^{\prime}}\Delta^{\prime}\right\\}|g_{\bm{k}}|^{2}$
$\displaystyle+\left[\Delta^{2}-\left(\frac{\Delta^{\prime}}{2}\right)^{2}\right]\left[{\gamma_{1}^{\prime}}^{2}+\Delta^{2}-\left(\frac{\Delta^{\prime}}{2}\right)^{2}\right]\,.$
(75)

Then the analytic expressions of the eigenenergies are

$\displaystyle\epsilon_{n\bm{k}}=\frac{1}{2}\left[\alpha_{n}\sqrt{-2x_{2}-\beta_{n}\frac{2x_{1}}{\sqrt{y}}-y}+\beta_{n}\sqrt{y}\right]\,,\quad\text{
for }n=1,2,3,4\,.$ (76)

with

$\displaystyle y$
$\displaystyle=\frac{1}{6}\left[4^{\frac{1}{3}}\left(y_{1}+\sqrt{y_{1}^{2}-4y_{2}^{3}}\right)^{\frac{1}{3}}+\frac{4^{\frac{2}{3}}y_{2}}{\left(y_{1}+\sqrt{y_{1}^{2}-4y_{2}^{3}}\right)^{\frac{1}{3}}}-4x_{2}\right]\,,$
(77) $\displaystyle y_{1}$
$\displaystyle=2x_{2}^{3}+27x_{1}^{2}-72x_{2}x_{0}\,,$ (78) $\displaystyle
y_{2}$ $\displaystyle=x_{2}^{2}+12x_{0}\,.$ (79)

At the Dirac points with $g_{\bm{k}}=0$, the four eigenenergies are
$\pm\Delta-\frac{\Delta^{\prime}}{2}$,
$\pm\sqrt{\Delta^{2}+\gamma_{1}^{2}}+\frac{\Delta^{\prime}}{2}$.

In general the electron-hole symmetry for AB-BG is broken due to the nonzero
of $\gamma_{4}^{\prime}$ and $\Delta^{\prime}$. However, we find that
$\gamma_{4}^{\prime}$ and $\Delta^{\prime}$ have negligble effects on the
optical transition between the bands (2, 3). With setting
$\gamma_{4}^{\prime}=0$ and $\Delta^{\prime}=0$, the eigenvalues become

$\displaystyle\epsilon_{n\bm{k}}$
$\displaystyle=\alpha_{n}\frac{1}{\sqrt{2}}\sqrt{z_{1}+\alpha_{n}\beta_{n}\sqrt{z_{2}}}\,,$
(80)

with

$\displaystyle z_{1}$
$\displaystyle={\gamma_{1}^{\prime}}^{2}+2\Delta^{2}+\left(2{\gamma_{0}^{\prime}}^{2}+{\gamma_{3}^{\prime}}^{2}\right)|g_{\bm{k}}|^{2}\,,$
(81) $\displaystyle z_{2}$
$\displaystyle=4{\gamma_{0}^{\prime}}^{2}\left[{\gamma_{3}^{\prime}}^{2}|g_{\bm{k}}|^{4}+2{\gamma_{1}^{\prime}}{\gamma_{3}^{\prime}}\text{Re}[g_{\bm{k}}^{3}]+({\gamma_{1}^{\prime}}^{2}+4\Delta^{2})|g_{\bm{k}}|^{2}\right]+\left({\gamma_{3}^{\prime}}^{2}|g_{\bm{k}}|^{2}-{\gamma_{1}^{\prime}}^{2}\right)^{2}\,.$
(82)

Obviously, the electronic states become electron-hole symmetric. Using Eq.
(80), we can have analytic discussion on the band gap $E_{g}$ and the VHS for
${\cal J}_{32}$. Around the Dirac point $\bm{K}$, the approximation
$g_{\bm{k}+\bm{K}}=-re^{i\theta}$ can be adopted for
$\bm{k}=\frac{2r}{\sqrt{3}a_{0}}(\cos\theta\hat{\bm{x}}+\sin\theta\hat{\bm{y}})$.
For zero $\Delta$, the zero energy of $\epsilon_{3\bm{k}}$ can be directly
found from Eq. (80) at $r=0$ or
$r=r_{0}=-\frac{\gamma_{1}^{\prime}\gamma_{3}^{\prime}}{{\gamma_{0}^{\prime}}^{2}}$
and $\theta=(2n+1)\pi/3$. Therefore, there exist in total four degenerate zero
energy points in one Dirac cone at $\Delta=0$; one is at this Dirac point, and
the other three locate along the K-M directions. Furthermore, for small $r$,
$\epsilon_{3\bm{k}}$ can be approximated by

$\displaystyle\epsilon_{3\bm{k}}^{2}=\Delta^{2}+c_{2}r^{2}+c_{3}\cos(3\theta)r^{3}+c_{4}r^{4}\,,$
(83)

with

$\displaystyle c_{2}$
$\displaystyle={\gamma_{3}^{\prime}}^{2}-\frac{4{\gamma_{0}^{\prime}}^{2}\Delta^{2}}{{\gamma_{1}^{\prime}}^{2}}\,,$
(84) $\displaystyle c_{3}$
$\displaystyle=-\frac{2{\gamma_{0}^{\prime}}^{2}{\gamma_{3}^{\prime}}}{\gamma_{1}^{\prime}}\,,$
(85) $\displaystyle c_{4}$
$\displaystyle=\frac{{\gamma_{0}^{\prime}}^{2}}{{\gamma_{1}^{\prime}}^{2}}\left[{\gamma_{0}^{\prime}}^{2}-2{\gamma_{3}^{\prime}}^{2}+\frac{4\Delta^{2}(2{\gamma_{0}^{\prime}}^{2}-{\gamma_{3}^{\prime}}^{2})}{{\gamma_{1}^{\prime}}^{2}}+\frac{16{\gamma_{0}^{\prime}}^{2}\Delta^{4}}{{\gamma_{1}^{\prime}}^{4}}\right]\,.$
(86)

From Eq. (83) the band structure around the Dirac points has following
features:

  1. 1.

For nonzero $\Delta$, the energy $\epsilon_{3\bm{k}}$ at the Dirac point
$\bm{K}$ is an extreme, and it is a local minimum (maximum) as $c_{2}>0$
($c_{2}<0$), which corresponds to $|\Delta|<\Delta_{c}$
($|\Delta|>\Delta_{c}$) with
$\Delta_{c}=|\gamma_{3}^{\prime}\gamma_{1}^{\prime}/(2\gamma_{0}^{\prime})|=0.0229$
eV.

  2. 2.

We first look at the case $|\Delta|>\Delta_{c}$ ($c_{2}<0$). For a fixed
$\theta$, $\epsilon_{3\bm{k}}$ around the Dirac point $\bm{K}$ has one more
local minimum located at $r=r_{e}(\cos 3\theta)$ with

$\displaystyle r_{e}(\cos 3\theta)=\frac{-3c_{3}\cos
3\theta+\sqrt{9c_{3}^{2}\cos^{2}3\theta-32c_{2}c_{4}}}{8c_{4}}\,.$ (87)

When $r$ is fixed and $\theta$ varies, $\epsilon_{3\bm{k}}$ has local maxima
as $\cos 3\theta=1$ and local minima as $\cos 3\theta=-1$. When both $r$ and
$\theta$ are considered, there exists a minimum at $r=r_{e}(-1)$ and
$\theta=(2n+1)\pi/3$ (along the K-$\Gamma$ directions for integer $n$), and a
VHS point at $r=r_{e}(1)$ and $\theta=2n\pi/3$ (along the K-M directions).

  3. 3.

For the case $|\Delta|<\Delta_{c}$ ($c_{2}>0$), $\epsilon_{3\bm{k}}$ has no
VHS point around the Dirac points but the minimum along K-$\Gamma$ directions
still exists.

  4. 4.

Similar analysis can be applied to study the JDOS ${\cal J}_{42}={\cal
J}_{31}$. After ignoring $\gamma_{4}^{\prime}$ and $\Delta^{\prime}$,
$\epsilon_{4\bm{k}}-\epsilon_{2\bm{k}}$ has a local minimum at the $\bm{K}$
point, and there is no VHS in ${\cal J}_{42}$. Therefore,
$\gamma_{4}^{\prime}$ and $\Delta^{\prime}$ play a key role in forming a VHS
in ${\cal J}_{42}$.

## References

  * Osterhoudt  _et al._ [2019] G. B. Osterhoudt, L. K. Diebel, M. J. Gray, X. Yang, J. Stanco, X. Huang, B. Shen, N. Ni, P. J. W. Moll, Y. Ran, and K. S. Burch, Colossal mid-infrared bulk photovoltaic effect in a type-I Weyl semimetal, Nature Materials 18, 471 (2019).
  * Sauer  _et al._ [2023] M. O. Sauer, A. Taghizadeh, U. Petralanda, M. Ovesen, K. S. Thygesen, T. Olsen, H. Cornean, and T. G. Pedersen, Shift current photovoltaic efficiency of 2D materials, npj Computational Mathematics 9, 35 (2023).
  * Spanier  _et al._ [2016] J. E. Spanier, V. M. Fridkin, A. M. Rappe, A. R. Akbashev, A. Polemi, Y. Qi, Z. Gu, S. M. Young, C. J. Hawley, D. Imbrenda, G. Xiao, A. L. Bennett-Jackson, and C. L. Johnson, Power conversion efficiency exceeding the Shockley-Queisser limit in a ferroelectric insulator, Nature Photonics 10, 611 (2016).
  * Salazar  _et al._ [2016] C. Salazar, J. L. Cheng, and J. E. Sipe, Coherent control of current injection in zigzag graphene nanoribbons, Physical Review B 93, 075442 (2016).
  * Wang  _et al._ [2019] K. Wang, R. A. Muniz, J. E. Sipe, and S. T. Cundiff, Quantum interference control of photocurrents in semiconductors by nonlinear optical absorption processes, Physical Review Letters 123, 067402 (2019).
  * Shockley and Queisser [1961] W. Shockley and H. J. Queisser, Detailed balance limit of efficiency of p‐n junction solar cells, Journal of Applied Physics 32, 510 (1961).
  * Kaner  _et al._ [2020] N. T. Kaner, Y. Wei, Y. Jiang, W. Li, X. Xu, K. Pang, X. Li, J. Yang, Y. Jiang, G. Zhang, and W. Q. Tian, Enhanced shift currents in monolayer 2D GeS and SnS by strain-induced band gap engineering, ACS Omega 5, 17207 (2020), pMID: 32715206.
  * Cook  _et al._ [2017] A. M. Cook, B. M. Fregoso, F. de Juan, S. Coh, and J. E. Moore, Design principles for shift current photovoltaics, Nature Communications 8, 14176 (2017).
  * Wei  _et al._ [2021] Y. Wei, W. Li, Y. Jiang, and J. Cheng, Electric field induced injection and shift currents in zigzag graphene nanoribbons, Physical Review B 104, 115402 (2021).
  * Xu  _et al._ [2021] H. Xu, H. Wang, J. Zhou, Y. Guo, J. Kong, and J. Li, Colossal switchable photocurrents in topological Janus transition metal dichalcogenides, npj Computational Materials 7, 1 (2021).
  * Yuan  _et al._ [2014] H. Yuan, X. Wang, B. Lian, H. Zhang, X. Fang, B. Shen, G. Xu, Y. Xu, S.-C. Zhang, H. Y. Hwang, _et al._ , Generation and electric control of spin–valley-coupled circular photogalvanic current in WSe2, Nature Nanotechnology 9, 851 (2014).
  * Arora  _et al._ [2021] A. Arora, J. F. Kong, and J. C. W. Song, Strain-induced large injection current in twisted bilayer graphene, Physical Review B 104, L241404 (2021).
  * Kaplan  _et al._ [2022] D. Kaplan, T. Holder, and B. Yan, Twisted photovoltaics at terahertz frequencies from momentum shift current, Physical Review Research 4, 013209 (2022).
  * Ji  _et al._ [2019] Z. Ji, G. Liu, Z. Addison, W. Liu, P. Yu, H. Gao, Z. Liu, A. M. Rappe, C. L. Kane, E. J. Mele, and R. Agarwal, Spatially dispersive circular photogalvanic effect in a Weyl semimetal, Nature Materials 18, 955 (2019).
  * Chan  _et al._ [2017] C.-K. Chan, N. H. Lindner, G. Refael, and P. A. Lee, Photocurrents in Weyl semimetals, Physical Review B 95, 041104 (2017).
  * Young and Rappe [2012] S. M. Young and A. M. Rappe, First principles calculation of the shift current photovoltaic effect in ferroelectrics, Physical Review Letters 109, 116601 (2012).
  * Pal  _et al._ [2018] S. Pal, A. B. Swain, P. P. Biswas, D. Murali, A. Pal, B. R. K. Nanda, and P. Murugavel, Giant photovoltaic response in band engineered ferroelectric perovskite, Scientific Reports 8, 8005 (2018).
  * Peng  _et al._ [2020] Y. Peng, X. Liu, Z. Sun, C. Ji, L. Li, Z. Wu, S. Wang, Y. Yao, M. Hong, and J. Luo, Exploiting the bulk photovoltaic effect in a 2D trilayered hybrid ferroelectric for highly sensitive polarized light detection, Angewandte Chemie 132, 3961 (2020).
  * Ji  _et al._ [2010] W. Ji, K. Yao, and Y. C. Liang, Bulk photovoltaic effect at visible wavelength in epitaxial ferroelectric BiFeO3 thin films, Advanced Materials 22, 1763 (2010).
  * Mele  _et al._ [2000] E. J. Mele, P. Král, and D. Tománek, Coherent control of photocurrents in graphene and carbon nanotubes, Physical Review B 61, 7669 (2000).
  * Cheng  _et al._ [2015] J. L. Cheng, N. Vermeulen, and J. E. Sipe, Third-order nonlinearity of graphene: Effects of phenomenological relaxation and finite temperature, Physical Review B 91, 235320 (2015).
  * Zheng  _et al._ [2022] Z. Zheng, Y. Song, Y. W. Shan, W. Xin, and J. L. Cheng, Optical coherent injection of carrier and current in twisted bilayer graphene, Physical Review B 105, 085407 (2022).
  * Sun  _et al._ [2010] D. Sun, C. Divin, J. Rioux, J. E. Sipe, C. Berger, W. A. de Heer, P. N. First, and T. B. Norris, Coherent control of ballistic photocurrents in multilayer epitaxial graphene using quantum interference, Nano Letters 10, 1293 (2010).
  * Sun  _et al._ [2012] D. Sun, C. Divin, M. Mihnev, T. Winzer, E. Malic, A. Knorr, J. E. Sipe, C. Berger, W. A. De Heer, P. N. First, _et al._ , Current relaxation due to hot carrier scattering in graphene, New Journal of Physics 14, 105012 (2012).
  * Vandelli  _et al._ [2019] M. Vandelli, M. I. Katsnelson, and E. A. Stepanov, Resonant optical second harmonic generation in graphene-based heterostructures, Physical Review B 99, 165432 (2019).
  * Brun and Pedersen [2015] S. J. Brun and T. G. Pedersen, Intense and tunable second-harmonic generation in biased bilayer graphene, Physical Review B 91, 205405 (2015).
  * Lin  _et al._ [2014] K.-H. Lin, S.-W. Weng, P.-W. Lyu, T.-R. Tsai, and W.-B. Su, Observation of optical second harmonic generation from suspended single-layer and bi-layer graphene, Applied Physics Letters 105, 151605 (2014).
  * Cheng  _et al._ [2017] J. Cheng, N. Vermeulen, and J. Sipe, Second order optical nonlinearity of graphene due to electric quadrupole and magnetic dipole effects, Scientific reports 7, 43843 (2017).
  * Shan  _et al._ [2018] Y. Shan, Y. Li, D. Huang, Q. Tong, W. Yao, W.-T. Liu, and S. Wu, Stacking symmetry governed second harmonic generation in graphene trilayers, Science Advances 4, eaat0074 (2018).
  * Xiong  _et al._ [2021] Y. Xiong, L.-k. Shi, and J. C. Song, Atomic configuration controlled photocurrent in van der Waals homostructures, 2D Materials 8, 035008 (2021).
  * Chaudhary  _et al._ [2022] S. Chaudhary, C. Lewandowski, and G. Refael, Shift-current response as a probe of quantum geometry and electron-electron interactions in twisted bilayer graphene, Physical Review Research 4, 013164 (2022).
  * Chuang  _et al._ [2012] Y.-C. Chuang, J.-Y. Wu, and M.-F. Lin, Analytical calculations on low-frequency excitations in AA-stacked bilayer graphene, Journal of the Physical Society of Japan 81, 124713 (2012).
  * Sipe and Shkrebtii [2000] J. E. Sipe and A. I. Shkrebtii, Second-order optical response in semiconductors, Physical Review B 61, 5337 (2000).
  * Gao  _et al._ [2020] Y. Gao, Y. Zhang, and D. Xiao, Tunable layer circular photogalvanic effect in twisted bilayers, Physical Review Letters 124, 077401 (2020).
  * McCann and Koshino [2013] E. McCann and M. Koshino, The electronic properties of bilayer graphene, Reports on Progress in Physics 76, 056503 (2013).
  * Rangel  _et al._ [2017] T. Rangel, B. M. Fregoso, B. S. Mendoza, T. Morimoto, J. E. Moore, and J. B. Neaton, Large bulk photovoltaic effect and spontaneous polarization of single-layer monochalcogenides, Physical Review Letters 119, 067402 (2017).