# Enhancing Shift Current via Virtual Multiband Transitions

Sihan Chen<EMAIL_ADDRESS>Kadanoff Center for Theoretical Physics,
University of Chicago, Chicago, Illinois 60637, USA  Department of Physics,
California Institute of Technology, Pasadena CA 91125, USA  Swati Chaudhary
Department of Physics, The University of Texas at Austin, Austin, Texas 78712,
USA  Department of Physics, Northeastern University, Boston, Massachusetts
02115, USA  Department of Physics, Massachusetts Institute of Technology,
Cambridge, Massachusetts 02139, USA  Gil Refael  Department of Physics,
California Institute of Technology, Pasadena CA 91125, USA  Institute for
Quantum Information and Matter, California Institute of Technology, Pasadena
CA 91125, USA  Cyprian Lewandowski<EMAIL_ADDRESS>National High
Magnetic Field Laboratory, Tallahassee, Florida, 32310, USA  Department of
Physics, Florida State University, Tallahassee, Florida 32306, USA

###### Abstract

Finding materials exhibiting substantial shift current holds the potential for
designing shift current-based photovoltaics that outperform conventional solar
cells. However, the myriad of factors governing shift current response poses
significant challenges in designing devices that showcase large shift current.
Here, we propose a general design principle that exploits inter-orbital mixing
to excite virtual multiband transitions in materials with multiple flat bands
to achieve enhanced shift current response. We further explicitly relate this
design principle to maximizing Wannier function spread as expressed through
the formalism of quantum geometry. We demonstrate the viability of our design
using a 1D stacked Rice-Mele model. Then, we consider a concrete material
realization - alternating angle twisted multilayer graphene (TMG) - a natural
platform to experimentally realize such an effect. We identify a new set of
twist angles at which the shift current response is maximized via virtual
transitions for each multilayer graphene and highlight the importance of TMG
as a promising material to achieve an enhanced shift current response at
terahertz frequencies. Our proposed mechanism also applies to other 2D systems
and can serve as a guiding principle for designing multiband systems that
exhibit enhanced shift current response.

##  I Introduction

The bulk photovoltaic effect (BPVE) is a promising alternative source of
photocurrent to conventional p-n junction based photovoltaics [1, 2, 3, 4, 5,
6, 7]. One of the microscopic mechanisms behind the BPVE is the shift current,
which generates direct current upon electromagnetic radiation in
noncentrosymmetric materials [2, 3]. Intuitively, the physical origins of
shift current trace back to a real-space shift experienced by an electron
wavepacket upon an interband excitation driven by linearly polarized light [8,
9]. Shift current is a second-order optical response whose contribution arises
from direct optical transitions and transitions via a virtual state [2, 10,
11, 12]. Multiple factors like interband velocity matrix elements, density of
states, number of bands, band gap as well as the quantum geometry[2, 13, 12,
14, 15, 16, 17] determine the magnitude of the shift current response. The
goal of establishing guiding principles to maximize the shift current response
is an ongoing problem with direct technological ramifications [10, 18, 19, 20,
19].

In multi-band systems, virtual transitions through intermediary bands are an
additional transition process that can contribute to the magnitude of
nonlinear optical response. In addition to a direct transition from initial
band to final band, virtual transition scatters through an intermediary band.
This manifests in the structure of the shift current response, which can be
decomposed into direct and virtual transition components [2, 11]. Previously,
the seminal work [8] outlined the design principles to maximize the shift
current from direct transitions. This work was inspired by ferroelectrics and
orthorhombic monochalcogenides [21, 22, 23] (materials with large band edge
responsivities which exhibit BPVE in the visible and far-infrared range). For
the two-band model under consideration in the Ref. 8 work, the virtual
transitions are not present, and thus the focus of Ref. [8] lies primarily on
maximizing singularities in the joint density of states (JDOS). Conversely, in
the case of materials with multiple bands, for example alternating angle moiré
graphene[24, 25, 26] (see Fig.1a-d for a typical bandstructure), it becomes
imperative to explore an alternative approach to maximizing shift current
response (e.g. Ref. [1, 27]), in particular based on the significance of
contributions stemming from virtual transitions.

In this paper, we demonstrate a viable mechanism that leverages virtual
transitions in multiband systems to significantly enhance shift current
response, see Fig.1e. The enhancement is achieved through designing materials
with multiple bands close in energy, thereby increasing the number of possible
virtual transitions through the increased number of intermediary bands.
Increasing the number of closely spaced energy bands however can carry
implications also for the quantum texture of the electron states as
exemplified by the Berry connection sum rules [28]. We demonstrate that shift
current, whose magnitude is partly controlled by such quantum geometric matrix
elements (e.g. Refs. [29, 30, 15]), can take advantage of this mechanism.
Specifically we find that the transition rate between the initial and final
bands involved in the photo-absorption process is enhanced in the regime where
wavefunction of the involved states (including the virtual state) become more
delocalized.

These two design principles are the key results of our paper, which underlying
physical mechanisms we elucidate first using a 1D multi-chain Rice-Mele model
and then we focus on the multilayer moiré graphene materals[24, 25, 26]. Moiré
materials such as twisted multilayer graphene (TMG) have acquired much
attention due to their multiple flat bands exhibiting nontrivial topology and
exotic phenomena such as superconductivity and correlated-insulating behavior
[31, 32, 33]. Our focus on the moiré materials stems from the fact that they
are expected to host a large shift current response[29, 30], but we emphasize
that the general principles we discuss apply to other materials well.

The paper is organized as follows. In Sec. II, we provide the theory for the
shift current responses. Next, in Sec. III, we present a one-dimensional toy
model to illustrate the principles for enhancing the contribution of virtual
transitions in the shift current response. We apply these ideas to TMG in Sec.
IV and study the dependence of shift current on number of layers and the
displacement field which controls the hybridization between different bands.
In Sec. V, we discuss the implications of our results for TMG and other
multilayer systems.

Figure 1: Bandstructure of T2G, T5G, T7G, T8G, and shift current conductivity
of T5G with varying external displacement field strength. a-d) Band structure
of TMG system with external displacement field strength $U=20$ meV at physical
twist angles $\theta=0.8^{\circ},1.32^{\circ},1.63^{\circ},1.71^{\circ}$ for
$N=2,5,7,8$, respectively. e) Shift current conductivity for T5G at
$\theta=1.32^{\circ}$ with varying external displacement field strength. The
result of increased external displacement field leads to band mixing, which
leads to enhancement of shift current via virtual transition as well as a new
peak at low frequency. Non-vanishing signal as $\omega\to 0$ in this figure is
a consequence of finite Lorentzian broadening (see Appendix .9 for more
details)

##  II Theory of Shift Current

Shift current is a second-order DC response to linearly polarized AC fields
[2, 11]. It is characterized by a third-rank conductivity tensor [2] given by

$\displaystyle\textbf{J}^{\beta}=2\sigma_{\alpha\alpha}^{\beta}(0,-\omega,\omega)E^{\alpha}(\omega)E^{\alpha}(-\omega)$
(1)

where $\textbf{J}^{\beta}$ is the $\beta$-th component of current density,
$E^{\alpha}(\omega)$ is electric field with frequency $\omega$ in the $\alpha$
direction, with $\alpha,\beta$ denotes spatial indices $x,y$. The shift-
current conductivity [8, 2, 11] is given by the expression:

$\displaystyle\sigma_{\alpha\alpha}^{\beta}(0,-\omega,\omega)=\frac{\pi
e^{3}}{\hbar^{2}}\sum_{m,n}\int
d^{2}\textbf{k}f_{mn}\textbf{S}_{mn}^{\beta\alpha}|\textbf{A}_{mn}^{\alpha}|^{2}\delta(\omega-\varepsilon_{mn}),$
(2)

where it has the form of a shift vector
$\mathbf{S}_{mn}^{\beta\alpha}=\mathbf{A}_{mm}^{\beta}-\textbf{A}_{nn}^{\beta}+\partial_{k_{\beta}}\text{Arg}(\textbf{A}_{mn}^{\alpha})$
that characterizes the shift in the localization of Bloch wavefunction upon
transition from state $\ket{u_{n}}$ to state $\ket{u_{m}}$, weighted by the
transition amplitude given by the interband Berry connection
$\textbf{A}_{mn}^{\alpha}=i\bra{u_{m}}\ket{\partial_{k_{\alpha}}u_{n}}$. The
sum is over all energy bands, where
$\varepsilon_{mn}=\varepsilon_{m}-\varepsilon_{n}$ is the energy difference
between two states, and $f_{mn}=f_{m}-f_{n}$ is the difference in occupancy of
their energy levels.

To visualize the transition process and draw a distinction between direct and
virtual transitions, we use generalized sum rules to replace wavefunction
derivatives with sums over all states of matrix element derivatives [2, 10,
11]. The shift current integrand defined as
$R_{mn}^{\alpha\alpha\beta}=|\textbf{A}_{mn}^{\alpha}|^{2}\textbf{S}_{mn}^{\beta\alpha}$
is given by,

$\displaystyle\begin{split}R_{mn}^{\alpha\alpha\beta}&=\frac{1}{\varepsilon_{mn}^{2}}\text{Im}\left[\frac{h_{mn}^{\alpha}h_{nm}^{\beta}\Delta_{mn}^{\alpha}}{\varepsilon_{mn}}+w_{mn}^{\alpha\beta}h_{mn}^{\alpha}\right]\\\
&+\frac{1}{\varepsilon_{mn}^{2}}\text{Im}\left[\sum_{l\neq
mn}\left(\frac{h_{nm}^{\alpha}h_{ml}^{\beta}h_{ln}^{\alpha}}{\varepsilon_{ml}}-\frac{h_{nm}^{\alpha}h_{ml}^{\alpha}h_{ln}^{\beta}}{\varepsilon_{ln}}\right)\right].\end{split}$
(3)

where $h_{mn}^{\alpha}=\matrixelement{m}{\partial_{\alpha}H}{n}$,
$\Delta_{mn}^{\alpha}=h_{mm}^{\alpha}-h_{nn}^{\alpha}$, and
$w_{mn}^{\alpha\beta}=\matrixelement{m}{\partial_{k_{\alpha}}\partial_{k_{\beta}}H}{n}$
(See Appendix .1 for derivation). The first term represents a direct
transition from band $n$ to band $m$, and the second term represents virtual
transitions through an intermediary band [11]. In a two-band system in 1D,
only $w_{mn}^{\alpha\beta}$ term contributes as the first term in the direct
transition term becomes purely real and virtual transitions are not present in
two-band model. In TMG system when we expand momentum to linear order near the
$K/K^{\prime}$ points, $w_{mn}^{\alpha\beta}$ vanishes. Fig. 1b shows a
schematic depiction of direct and virtual transitions between different bands
in TMG. The virtual transition sums over all intermediary bands, and we expect
such a process to be the dominant contribution in systems containing multiple
bands of similar energies. The coupling between these bands can be controlled
by the external displacement field which can be used to enhance the shift
current conductivity as shown in Fig. 1e. In addition to discerning the direct
and virtual transition contributions, Eq. 3 is numerically more amenable as it
avoids dealing with the issue of gauge fixing when evaluating the derivative
of the wavefunction directly [10, 11].

##  III Multilayer Rice-Mele Model

In order to study the role of virtual transitions in controlling the magnitude
of the shift current response, we construct a toy model with multiple close
spaced bands. Specifically, we consider a 1D stacked Rice-Mele (RM) model to
demonstrate the interplay between virtual transitions and bandstructures in
multilayer systems. The Rice-Mele model is a prototypical model for one-
dimensional ferroelectrics [8] represented by the Hamiltonian

$\displaystyle
H=\sum_{i}\left[\left(\frac{t}{2}+(-1)^{i}\frac{\delta}{2}\right)(c_{i+1}^{\dagger}c_{i}+\text{h.c.})+(-1)^{i}\Delta
c_{i}^{\dagger}c_{i}\right],$ (4)

where $t$ is the hopping potential, $\delta$ parameterizes the difference in
hopping strength between the two neighboring sites with unit cell length $a$,
and $\Delta$ is the staggered onsite potential. The canonical Bloch
Hamiltonian in momentum space is given by [34]

$\displaystyle
H(k)=\sigma_{x}t\cos{\frac{ka}{2}}-\sigma_{y}\delta\sin{\frac{ka}{2}}+\sigma_{z}\Delta.$
(5)

To construct a multilayer Rice-Mele model, we denote the dependence of the
momentum-space Bloch Hamiltonian on its parameters as
$H_{i}=H_{i}(t_{i},\delta_{i},\Delta_{i})$ where $i$ is the layer index, and
with the full Hamiltonian for $N$ layers taking the following form:

$\displaystyle H_{\text{MRM},0}(k)=H_{1}(k)\oplus H_{2}(k)\oplus\cdots\oplus
H_{N}(k)$ (6)

As written above in Eq. 6, the multilayer Rice-Mele Hamiltonian is block-
diagonal thus transitions between the different blocks are impossible. To
enable virtual transitions between the RM sectors belonging to different
layers we introduce couplings between layers. Specifically, consider the case
of $N=2$ with band mixing, the 2-layer Rice-Mele (2RM) model has Hamiltonian
of the form

$\displaystyle H_{2\text{RM}}(k)=\begin{pmatrix}H_{1}(k)&\varepsilon V\\\
\varepsilon V^{\dagger}&H_{2}(k)\end{pmatrix},$ (7)

where $V$ is the operator that characterizes interactions and $\varepsilon$
parameterize the mixing strength. Fig. 2a shows the bandstructure of 2RM model
with $V=\sigma_{x}$, and
$(t,\delta_{1},\Delta_{1},\delta_{2},\Delta_{2})=(1,0.8,0.7,0.86,0.6)$ with
(without) mixing as depicted by solid (dashed) line. For this particular form
of coupling, the two low (high) energy bands hybridize which widens the gap
between black and red in low (high) energy sector. This effect manifest most
strongly at $ka=\pm\pi$ where the bands flatten out significantly in
comparison to the unmixed case as shown in Fig. 2a. Here we fixed $t$ to be
the same for both RM models and vary the hopping amplitude by tuning
$\delta_{i}$.

Figure 2: Bandstructure and shift current conductivity of 2RM model and its
dependence on model parameters $\delta_{2}$. a) Band structure of 2RM model
without band mixing (dashed line) and with band mixing of $\varepsilon=0.1$
(solid line). The red/black band refers the first/second RM model, with
intracell hopping strength $v_{i}=\frac{1}{2}(t_{i}+\delta_{i})$, intercell
hopping strength $w_{i}=\frac{1}{2}(t_{i}-\delta_{i})$, and sublattice offset
potential $\Delta_{i}$, where $i=1,2$ indexing the different RM layers. A
schematic description of the model is included, with $\varepsilon$ coupling
the different sublattice degree of freedom between the chain. The parameters
of 2RM are
$(t,\delta_{1},\Delta_{1},\delta_{2},\Delta_{2})=(1,0.8,0.7,0.86,0.6)$. b)
Shift current conductivity due to transitions between the red band (first RM
model) with and without band mixing. The presence of mixing enhances the peak
due to transition at band edge via virtual transitions through the second RM
bands. c) Bandstructure of the top two bands, with the orange dashed line
denoting chemical potential at filling $\nu=0.5$. d) Shift current
conductivity due to transitions from band $2$ to band $3$, showing a large
response at low frequency determined by the energy gap of band $2$ and $3$ at
$ka=\pm\pi$.

The coupled bands after mixing enable transitions to occur among all four
bands, allowing virtual transition to enhance the overall shift current
response. As shown in Fig. 2b, the transition from valence to conduction band
of the first RM model is enhanced after band mixing due to virtual transitions
through a different RM sector. The peaks in the response are determined by the
energy gaps at $ka=\pm\pi,0$, where the joint density of state is the largest.
Level repulsion as a result of band mixing widens the gap, consequently
shifting the location of the peak to a slightly higher frequency. In addition
to the selective enhancment at precise frequencies, the overall shift current
conductivity $\sigma(\omega)$ over whole frequency range is increased.
Specifically, using a “figure of merit”

$\displaystyle M=\left|\int\sigma d\omega\right|$ (8)

as a measure of the overall shift current response [27], we find that
$M_{\text{unmix}}\approx 9$ and $M_{\text{mix}}\approx 13.5$, demonstrating an
enhanced response.

The coupling-induced hybridization of bands not only flattens the dispersion
around $ka=\pm\pi$ but also engenders additional transition pairs throughout
the Brilliouin zone that further enhance the overall shift current response.
The total shift current conductivity is now obtained by summing over all
possible transition pairs, which include transitions from one RM sector to
another. The sign of conductivity depends on the specific transition pairs. On
the basis of quantity $M$, introduced in Eq. 8 above, some transition pairs
will contribute towards the cancellation of conductivity. Nevertheless, an
overall enhancement is still observed compared to models without band mixing
(See Appendix .5 for more details). Different transition pairs give rises to
peaks in conductivity at different frequencies range. To avoid cancellation
due to opposite signs, one could work at specific frequencies to excite the
corresponding transition pairs.

Furthermore, by studying the shift current response of 2RM model at filling
$\nu=0.5$, where the chemical potential lies in between the band edge of the
top two bands, see Fig. 2c, we find an additional shift current response
between the two bands 2 and 3 with $M_{32}\approx 42.4$ in the low frequency
range as shown in Fig. 2d. This response relies on band mixing and is much
larger than the response at charge neutrality considered in the previous
paragraph and in Fig. 2b. This response is dominated by direct transition as
virtual transitions are suppressed by the large energy gap. Only one peak
occurs at frequency determined by the energy gap difference at the band edge
because the states around $k=0$ do not contribute to transitions. The reasons
for such a large response are twofold: small energy gap and large joint
density of state, which can be inferred from Eq. 3 and Eq. 2, respectively.
The exclusion of states around $k=0$ would not drastically reduce the shift
current response, because the gap is much larger at $k=0$ in comparison to the
gap at the band edge.

We can understand the above results at $\nu=0.5$ with help of the Ref. [27].
As pointed out in Ref. [27], the energy gap is not the only criteria which
dictates the shift current response but band dispersion also plays an
important role. In a typical tight-binding two-band model, the structure of
the valence band fully determines the structure of the conduction band. This
unique relationship defines an upper limit on quantity $M$ in terms of energy
gap, bandwidth, and the range of hopping as discussed in Ref. [27]. These
three factors also dictate the extent of delocalization of the wavefunction in
a two-band model which was highlighted[27] to be the most important criteria
for enhanced shift current response. On the other hand, for a system with
energetically close multiple bands, there are additional physical principles
that control the limit of shift current response. In particular, many recent
works have demonstrated how quantum geometry which manifests in multi-band
systems can contribute to physical quantities, for example such as the
superfluid stiffness, beyond their nominal limits inherent in single-band or
two-band models [33, 35, 36, 37, 38].

Figure 3: Parameter choice that leads to enhanced shift current response. a)
Total shift current $M$ broken down into direct and virtual contributions. b)
Shift vector before and c) after band mixing, demonstrating contribution
mainly coming from band edge. d) FS metric, e) shift current integrand, and f)
shift vector at $ka=-\pi$, showing that enhanced response corresponds to a
larger spread in the Wannier wavefunction.

While the multi-band aspect of our system manifests explicitly through virtual
transition contributions to shift-current response, the specific choice of
model parameter is crucial in achieving enhanced response via virtual
transitions. Fig. 3a compares the magnitude of total shift current response
$M$ with and without band mixing of transition from band $0$ to band $3$
fixing $\delta_{1},\Delta_{1},\Delta_{2},\varepsilon$ and vary $\delta_{2}$,
breaking down the contribution to total response into direct and virtual
transition. We see that the effect of band mixing can lead to a reduction in
direct transition contribution to shift current. Simultaneously however there
can be a significant enhancement in the virtual transition contributions to
the shift current.

Specifically there exist parameter regime, where this increase in virtual
transition component is enough to offset the reduction in direct transition.
In this region, where we picked $\delta_{2}\approx 0.86$, we observe an
enhanced shift vector, $S_{30}(k)$, near the band edge, as shown in Fig. 3b,c.
We attribute this enhancement effect to the sudden increase in the spatial
delocalization of Wannier wavefunctions stemming from strong band
hybridization at those band parameter values.

This Wannier function delocalization can be studied more systematically with
help of the Fubini-Study (FS) metric [39, 35]. (See Appendix .2 for more
details) The FS metric defines a distance between momentum states, and its
integral over the Brillouin zone serves as a lower bound for Wannier function
localization [40, 38, 41]. Fig. 3d shows a sweep across the region of
$\delta_{2}$ that exhibits enhancement in total shift current via virtual
transitions, and we find a corresponding increase in the FS metric integral
$\int g(k)dk$, which implies a delocalization of Wannier function. This
observation is in agreement with previous works that demonstrate a
correspondence between larger delocalization and larger shift vector[1], and
consequently large shift current response [42]. The blue dashed line denotes
the maximum of the integrated FS metric and where $S_{30}$ from Fig. 3c
diverges at the band edge. The blue dashed line is where the band gap of the
top/bottom is minimized, and it also corresponds to optical zero ($A_{mn}(k)$
= 0). Hence, we identify a competition between the increasing of $S_{mn}(k)$
and decreasing of $|A_{mn}(k)|$, leading to maximum of shift vector integrand
$R_{mn}$ defined in Eq. 3 occurring at the parameter near optical zeros
(marked by orange dashed line), where the shift vector is still large, see
Fig. 3d-f.

We further note that the enhancement in FS metric $g(k)$ can be described by a
Lorentzian function. (See Appendix .3 for more details). This suggests that
the sudden increase in FS metric, and consequently the spread of Wannier
wavefunction, can be thought of as resonance in parameter space. It is this
resonantly enhanced delocalization of Wannier wavefunction that leads to the
enhanced shift current response. In momentum space, this enhancement is due to
virtual transition via intermediary bands. Specifically, the enhancement comes
from the small energy gap between the top two bands $E_{32}$ and the large
matrix element $\text{Im}\\{v_{20}\\}$ with a pole occurring at
$\delta_{2}=\sqrt{\delta_{1}^{2}+\Delta_{1}^{2}-\Delta_{2}^{2}}$ at the band
edge. (See Appendix .4 for more details).

The proposed mechanism above suggests to stack multiple 1D systems together
and induce band mixing. By careful tuning the system’s parameters, we find a
region of parameters where the Wannier function becomes maximally delocalized.
Although our analysis is constrained to the simple 1D model, we propose that
the identified design principles can extend to more realistic
models.Specifically, in the following section we demonstrate these principles
with the help of a realistic material model.

Figure 4: Bandstructure of T5G and its total shift current conductivity plot
at various filling, and breakdown of the enhancement effect via virtual
transition of two flat bands at $\nu=2$. a) A schematic depiction of
transition between flat to flat bands and from flat to dispersive bands, as
well as bandstructure of T5G with $U=20$ meV (solid line) and with $U=0meV$
(dashed line) at $\theta=1.32^{\circ}$. b) Shift current conductivity at
various system fillings, demonstrating the ability to optimize conductivity at
different frequency by tuning the system filling. c) Shift current
conductivity due to transitions between bands colored in red in a) at $\nu=2$.
Here we break down the total conductivity into virtual transition through flat
bands and virtual transition through dispersive bands. Virtual transitions
enhance both flat-flat and flat-dispersive contribution to shift current.

##  IV Twisted Multilayer Graphene

Twisted bilayer graphene (TBG) has been proposed to exhibit large shift
current response due to its nontrivial flatband topology [29, 30]. We thus
propose to use TBG as a viable building block to construct the above proposed
multilayer system. Specifically we will follow the stacking procedure
introduced in Ref. [24], that has been recently realized experimentally[25,
26]. The single-particle spectrum of TBG is described by an effective
continuum model (e.g. Ref. [43]). Consider Hamiltonian in the sublattice basis
$(A_{1},B_{1},A_{2},B_{2})$,

$\displaystyle H_{\theta}=\begin{pmatrix}H_{1}&U^{\dagger}\\\
U&H_{2}\end{pmatrix}.$ (9)

Let $l=1,2$ denote the bottom and top layer, respectively. The intralayer
Hamiltonian centered at Dirac cone $\mathbf{K}_{\xi}^{(l)}$ given by

$\displaystyle H_{l}=-\hbar
v[R(\pm\theta/2)(\mathbf{k}-\mathbf{K}_{\xi}^{(l)})]\cdot(\xi\sigma_{x},\sigma_{y})+\Delta_{l}\sigma_{z}$
(10)

where $\theta$ is twist angle relative to the origin, $+/-$ corresponds to
$l=1/2$, and $\xi=\pm 1$ labels the valley index of the original graphene
Dirac cone. We include the sublattice offset term $\Delta_{l}\sigma_{z}$ to
break inversion $C_{2z}$ symmetry, which is essential for non-zero shift
current response[7]. Experimentally, this is achieved by coupling TBG to
top/bottom layer substrate. The effective interlayer coupling $U$ has the form

$\displaystyle U$
$\displaystyle=\begin{pmatrix}U_{A_{2},A_{1}}&U_{A_{2},B_{1}}\\\
U_{B_{2},A_{1}}&U_{B_{2},B_{1}}\end{pmatrix}$ (11)
$\displaystyle=\begin{pmatrix}u&u^{\prime}\\\
u^{\prime}&u\end{pmatrix}+\begin{pmatrix}u&u{{}^{\prime}}w^{-\xi}\\\
u{{}^{\prime}}w^{\xi}&u\end{pmatrix}e^{i\xi\mathbf{G}_{1}^{M}\cdot\mathbf{r}}$
(12) $\displaystyle+\begin{pmatrix}u&u{{}^{\prime}}w^{\xi}\\\
u{{}^{\prime}}w^{-\xi}&u\end{pmatrix}e^{i\xi(\mathbf{G}_{1}^{M}+\mathbf{G}_{2}^{M})\cdot\mathbf{r}},$
(13)

where $w=e^{i2\pi/3}$, and
$\mathbf{G}=n_{1}\mathbf{G}_{1}^{M}+n_{2}\mathbf{G}_{2}^{M}$ is the moiré
reciprocal lattice vector, for $n_{1},n_{2}\in\mathbb{Z}$. We choose
$\mathbf{G}_{i}^{M}=R(-\theta/2)\mathbf{G}_{i}-R(\theta/2)\mathbf{G}_{i}$,
with the $\mathbf{G}_{1}=(2\pi/a)(1,-1/\sqrt{3})$ and
$\mathbf{G}_{2}=(2\pi/a)(0,2/\sqrt{3})$. Following the work done in [29], we
choose $u{{}^{\prime}}=90meV$ and $u=0.4u^{\prime}$, sublattice offset
$\Delta=5meV$, and set $\hbar v/a=2.1354eV$ and $a=0.246nm$ however the
qualitative behaviour of the shift current that will be discussed below does
not depend on these specific parameters of the model.

We alternatively twist each layer of graphene by an angle $\theta/2$ with
respect to each adjacent layer to extend the bilayer to a multilayer system
[24, 25, 26]. Coupling is only considered between adjacent layers. The
Hamiltonian in the layer basis is given by

$\displaystyle H_{n}=\begin{pmatrix}H_{1}&U^{\dagger}&0&0&\cdots\\\
U&H_{2}&U&0&\cdots\\\ 0&U^{\dagger}&H_{1}&U^{\dagger}&\cdots\\\
\vdots&\vdots&\vdots&\vdots&\ddots\end{pmatrix}$ (14)

where we treated tunneling terms between all layers to be the same, and the
angle of $H_{n}$ is determined by
$\theta_{n}=2\cos(\frac{\pi}{n+1})\theta_{\text{TBG}}$. In the multilayer
model, sublattice offset term is added at the top and bottom layers to
incorporate the effect of coupling to hBN substrate. This multilayer
Hamiltonian can be block-diagonalized into a direct sum of effective TBG pairs
under a unitary transformation [24]. To induce virtual transition among
different Hamiltonian layers, we consider applying an external displacement
field on the top and bottom gates. (See Appendix .6 for details) The
bandstructure of $N=2,5,7,8$ is shown in 1a-d, where we see $\lfloor
n/2\rfloor$ number of TBG-like bands with an additional Dirac cone if $n$ is
odd. Bandstructure of T5G is reproduced in Fig. 4a with an external
displacement field, including a schematic description of flat-to-flat band
transition and flat-to-dispersive band transition. We study the transition
between the lowest flat bands (colored in red) and observe an enhancement in
shift current via virtual transitions to nearby flat bands once the external
displacement field is included. To note, our TMG model has $C_{3z}$ symmetry,
and by group theoretical method we can show that only two components of the
conductivity tensor are independent (See Appendix .8 for details).

In our analysis, we neglect the role of electron-electron interactions, as our
goal is to demonstrate just the principle of the physical mechanism that
causes the shift current enhancement. The role of interactions in TBG and
moiré materials is extensive as studied both theoretically (e.g. Refs. [44,
45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58]) and experimentally
(e.g. Refs. [59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 69, 70, 71]).
However, the exchange-driven physics, particularly at high temperatures
($T>10$ K), should not affect the qualitative ideas discussed in this
manuscript. Moreover, as argued in Ref. [58], interactions can, in fact, give
rise to a self-generated “displacement” field, which will induce hybridization
between different energy sectors even in the absence of an external field. For
simplicity of the discussion in what follows, we focus on the non-interacting
model, leaving a detailed study of interactions on the shift-current
photoresponse of a multilayer moiré system to future works.

The physical twist angle in TMG determines the effective twist angles of TBG-
like bands, which in turn controls the energy difference between two flat band
pairs in each energy sector of the model. The magnitude and direction of
virtual transition are determined by the velocity matrix element and the
energy gap of bands involved in the optical transition process as implied by
Eq. 3. To receive maximum enhancement from virtual transitions, we choose
twist angles such that the energy gap between two flat bands is as small as
possible while the velocity matrix element is nonzero. The interplay between
velocity matrix element and energy gap determines the optimal twist angle for
each multilayer system. Crucially here, since our goal is not to achieve a
flat band, but rather to obtain closely spaced energy bands, the relevant
twist angles differ from the magic angles of these multilayer systems. For
example, from our numerics, we found $\theta_{5}=1.32^{\circ}$ for T5G helps
optimize the shift current reponse, while the magic angle for T5G occurs at
$\theta_{5,m}=1.73^{\circ}$, where perfect flat band emerges. (See Appendix .7
for more details) Correspondingly, there is an increase in shift current
response due to reduced band gap and large joint density of states, an
enhancement in direct transition. At the magic angle, the strongly interacting
nature of the system requires us to incorporate Coulomb interaction, which
enhances the resonance at non-interacting characteristic frequency due to
filling-dependent band renormalization [29]. In TBG away from the magic angle,
the shift current response is qualitatively the same, scales according to band
gap and moiré length. In T5G, however, an enhancement via virtual transition
to nearby flat bands occurs at angles where the two flat bands become close in
energy, in addition to a new resonance peak at low frequency.

In contrast to the simple 4-band model (the 2RM model) described in last
section, the TMG system has a contribution to shift current conductivity from
both direct transition from flat to flat bands and virtual transitions from
flat to dispersive bands coming from both the same and different TBG-like
sectors. In the absence of external displacement field, it is the virtual
transition to dispersive bands within the same TBG sector that dominates the
shift current response. With external displacement field, band mixing enables
transitions to access different TBG-like sectors in the form of virtual
transitions to either flat or dispersive bands of the other TBG-like sectors.
Shift current conductivity $\sigma_{xxy}$ at $\nu=2$ is shown in Fig. 4c.
After band mixing, virtual transition through intermediary bands significantly
enhances the response, in both flat-to-flat and flat-to-dispersive bands. We
notice a shift in the peak of shift current response to higher $\omega$ with
band mixing. Without mixing, the dominant contribution to shift current comes
from flat-to-dispersive band transition within the same TBG sector, with peak
determined by the energy difference between the flat and dispersive band at
the $\gamma$ point. After band mixing, virtual transition via nearby flat
bands becomes the dominant contribution, which enhances the response occurring
that peak near $\mu$, where van Hove singularity occurs. Our choice of T5G
angles gives rise to a larger energy gap at $\mu$ point between the two flat
bands than the energy gap between flat and dispersive band at $\gamma$ point,
resulting in a shift to higher $\omega$ in shift current peak.

Summing up all transition pairs reveals a giant peak at low frequency in shift
current conductivity due to transitions to nearby flat bands and Dirac cone
around the $\kappa$ points. As shown in Fig. 4b) for various fillings at
$U=20meV$, the peak of the response occurs at a frequency determined by the
energy gap of those bands at the $\kappa$ points, roughly determined by the
sublattice offset term $\Delta$. Previous work [72] also pointed out the
effect of low-frequency divergence in topological semi-metals, where they
attributed the divergence to singularity in quantum geometry. At a given
filling, in addition to the giant peak at low frequency, there is a secondary
peak occurring at frequency proportional to the energy gap of the two flat
bands at $\mu$, which corresponds to transitions among the flat bands pairs.
The appearance of two distinct peaks is a result of summing over different
transition pairs. More explicitly, transitions between top/bottom pairs of
flat band give rise to low-frequency divergence, whereas transition from
bottom to top flat bands give rise to the secondary peaks. This is another
unique feature to TMG system compared to TBG system studied in previous works
[29, 30].

##  V Discussion

In this work, we discussed new design principles for optimizing shift current
response in multi-band systems. In particular, we demonstrated a mechanism
based on enhancement of virtual transition contributions in systems with
multiple layers. The enhancement arises mainly from the mixing of bands
belonging to different layers. The mechanism is explained using a 1D toy model
where the stacking order provides an additional “knob” to tune the shift-
current response.

We have shown that in TMG systems, shift current response is enhanced by both
virtual transitions among nearby bands and the additional transition pairs
between additional bands. Virtual transitions enhance resonance at the
original peaks whereas transitioning among nearby bands around $\kappa$ point
produces another peak in conductivity at low frequency. Both phenomena of
enhancement occur over a wide range of fillings, with different fillings
corresponding to a slightly different enhanced peak in conductivity. The
effect of virtual transitions found in this paper relies on the existence of
multiple flat bands and applies when the number of layers $N\geq 5$. Compared
to TBG, T3G only has an additional Dirac cone at $\kappa$ point, and although
T4G has two flat band pairs, the second pair is too dispersive for the
enhancement effect to occur. This makes shift current response in TMG for
$N\geq 5$ qualitatively different from that of $N<5$. As we increase $N$, the
existence of additional flat band pairs can lead to further enhancement via
virtual transitions. The physical angle at which we would observe enhanced
shift current response is also expected to increase, making it potentially
easier to realize experimentally.

Much of the experimental attention in the studies of TMG has been focused on
the system’s behavior near magic angles (e.g. Refs. [59, 60, 61, 62, 63, 64,
65, 66, 67, 68, 69, 69, 70, 71]), where superconductivity and correlated
insulating behavior are revealed. In TMG systems, in addition to the set of
magic angles that give rise to superconductivity and correlated insulating
behavior, we find a different set of angles that will lead to a giant shift
current response. (See Appendix .7 for details) At these angles, we find
enhanced virtual transitions via nearby bands, as well as more delocalized
wavefunctions. These angles are relevant in order to achieve large nonlinear
optical response. These findings make TMG a promising platform for designing
devices that produce giant photocurrents at terahertz frequency.

Finally, we believe that this procedure of constructing multilayer systems can
be used as a general guide to enhance shift current response in other 2D
systems. Given a generic 1D or 2D system that already exhibits a large shift
current response, we can further enhance its response through a stacking
procedure that induces virtual transition between neighbouring energy bands.
It is particularly appealing to twisted materials as this procedure of twist
and stack is the commonly used technique of constructing this class of
materials. By varying the twist angle, one can find a setup where shift
current response is enhanced.

Motivated by Eq. 3, we see that to achieve a large shift current, we want to
minimize the band gap $\varepsilon_{mn}$. However, simply engineering bands
with a small band gap will also shift the frequency peak to small $\omega$,
making the current susceptible to thermal fluctuations at finite temperature
[72]. In addition, the derivation of shift current expression via perturbation
theory assumes small $|E/\omega|$, leading to worries concerning the breakdown
of perturbation theory in the small $\omega$ regime. By keeping the large band
gap between band $m$ and $n$ and engineering nearby bands $l$ close in energy
such that $\varepsilon_{ml}$ is small, we leverage virtual transitions to
enhance shift current response while maintaining transition at relatively
large $\omega$. Furthermore, the FS metric is also inversely proportional to
the energy separation of two bands. When a pair of bands with non-zero
interband Berry connection have a small energy gap, it leads to a large FS
metric, which also corresponds to a large localization length of Wannier
wavefunction and a large shift vector. This confirms the intuition that a
large shift current arises when there is a large shift in the center of charge
of the wavefunction [8, 9] in real space upon excitation to the conduction
band. Recently, nonlinear optical response has been formulated in the language
of quantum geometry [73], and it has been shown that shift current is related
to the quantum geometric connection. From this perspective, the meticulous
adjustment of model parameters positions the material in close proximity to
the geometric singularity in the band structure, leading to divergence of the
geometric quantities that govern the shift current response.

The intricate interplay between shift current and quantum geometry has enabled
its use as a probe of quantum geometry and interactions in topological
materials [29]. Here we show two examples of enhanced shift current response
via virtual transition near the band closing point. In topological materials,
gap closing is associated with topological phase transition. It would be
interesting to study shift current response in noncentrosymmetric multiband
topological materials near critical point to see if we observe a corresponding
enhanced response. This could potentially allow us to use shift current as a
probe of topological phase transitions.

One could also imagine our design principle to be utilized to enhance other
optical responses that demonstrate a geometric origin. Another common second-
order bulk photovoltaic response is the injection current, which is shown to
be related to quantum geometric tensor [73, 72]. Engineering bands near
geometrical singularity might also lead to large injection current response.

Furthermore, quantum geometric effects are not limited to optical responses
but also lead to many other exotic phenomena like flat-band superconductivity
[36, 37, 38, 74], undamped plasmons [75], non-reciprocity in plasmonics [76,
77] and in Landau-Zener tunneling [78]. The mechanism illustrated here for
enhancement of shift-current response using multi-band systems would
potentially be applicable in a wide variety of quantum geometric effects.

## Acknowledgement

We thank Roshan Krishna Kumar for useful discussions and collaboration on a
related project. S.C. acknowledges support from the Summer Undergraduate
Research Fellowship at Caltech. Swati Chaudhary acknowledges support from
National Science Foundation through the Center for Dynamics and Control of
Materials: an NSF MRSEC under Cooperative Agreement No. DMR-1720595. G.R.
expresses gratitude for the support by the Simons Foundation, and the ARO MURI
Grant No. W911NF-16-1- 0361 and the Institute of Quantum Information and
Matter. C.L. was supported by start-up funds from Florida State University and
the National High Magnetic Field Laboratory. The National High Magnetic Field
Laboratory is supported by the National Science Foundation through
NSF/DMR-2128556 and the State of Florida.

## References

  * Dai and Rappe [2023] Z. Dai and A. M. Rappe, Recent progress in the theory of bulk photovoltaic effect, Chemical Physics Reviews 4, 10.1063/5.0101513 (2023).
  * Sipe and Shkrebtii [2000] J. E. Sipe and A. I. Shkrebtii, Second-order optical response in semiconductors, Phys. Rev. B 61, 5337 (2000).
  * von Baltz and Kraut [1981] R. von Baltz and W. Kraut, Theory of the bulk photovoltaic effect in pure crystals, Phys. Rev. B 23, 5590 (1981).
  * Fridkin [2001] V. Fridkin, Bulk photovoltaic effect in noncentrosymmetric crystals, Crystallography Reports 46, 654 (2001).
  * Nastos and Sipe [2006] F. Nastos and J. E. Sipe, Optical rectification and shift currents in gaas and gap response: Below and above the band gap, Phys. Rev. B 74, 035201 (2006).
  * Pusch  _et al._ [2023] A. Pusch, U. Römer, D. Culcer, and N. J. Ekins-Daukes, Energy conversion efficiency of the bulk photovoltaic effect, PRX Energy 2, 013006 (2023).
  * Nagaosa [2022] N. Nagaosa, Nonlinear optical responses in noncentrosymmetric quantum materials, Annals of Physics 447, 169146 (2022).
  * Fregoso  _et al._ [2017] B. M. Fregoso, T. Morimoto, and J. E. Moore, Quantitative relationship between polarization differences and the zone-averaged shift photocurrent, Phys. Rev. B 96, 075421 (2017).
  * Wang and Qian [2019] H. Wang and X. Qian, Ferroicity-driven nonlinear photocurrent switching in time-reversal invariant ferroic materials, Science advances 5, eaav9743 (2019).
  * 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, 1 (2017).
  * Parker  _et al._ [2019] D. E. Parker, T. Morimoto, J. Orenstein, and J. E. Moore, Diagrammatic approach to nonlinear optical response with application to weyl semimetals, Phys. Rev. B 99, 045121 (2019).
  * Morimoto and Nagaosa [2016] T. Morimoto and N. Nagaosa, Topological nature of nonlinear optical effects in solids, Science advances 2, e1501524 (2016).
  * Orenstein  _et al._ [2021] J. Orenstein, J. Moore, T. Morimoto, D. Torchinsky, J. Harter, and D. Hsieh, Topology and symmetry of quantum materials via nonlinear optical responses, Annual Review of Condensed Matter Physics 12, 247 (2021).
  * Ma  _et al._ [2021] Q. Ma, A. G. Grushin, and K. S. Burch, Topology and geometry under the nonlinear electromagnetic spotlight, Nature materials 20, 1601 (2021).
  * Ahn  _et al._ [2022] J. Ahn, G.-Y. Guo, N. Nagaosa, and A. Vishwanath, Riemannian geometry of resonant optical responses, Nature Physics 18, 290 (2022).
  * Morimoto  _et al._ [2023] T. Morimoto, S. Kitamura, and N. Nagaosa, Geometric aspects of nonlinear and nonequilibrium phenomena, Journal of the Physical Society of Japan 92, 072001 (2023).
  * Wang  _et al._ [2022] H. Wang, X. Tang, H. Xu, J. Li, and X. Qian, Generalized wilson loop method for nonlinear light-matter interaction, npj Quantum Materials 7, 61 (2022).
  * Wang  _et al._ [2016] F. Wang, S. M. Young, F. Zheng, I. Grinberg, and A. M. Rappe, Substantial bulk photovoltaic effect enhancement via nanolayering, Nature communications 7, 10419 (2016).
  * Dong  _et al._ [2023] Y. Dong, M.-M. Yang, M. Yoshii, S. Matsuoka, S. Kitamura, T. Hasegawa, N. Ogawa, T. Morimoto, T. Ideue, and Y. Iwasa, Giant bulk piezophotovoltaic effect in 3r-mos2, Nature nanotechnology 18, 36 (2023).
  * Grinberg  _et al._ [2013] I. Grinberg, D. V. West, M. Torres, G. Gou, D. M. Stein, L. Wu, G. Chen, E. M. Gallo, A. R. Akbashev, P. K. Davies, _et al._ , Perovskite oxides for visible-light-absorbing ferroelectric and photovoltaic materials, Nature 503, 509 (2013).
  * Singh and Hennig [2014] A. K. Singh and R. G. Hennig, Computational prediction of two-dimensional group-iv mono-chalcogenides, Applied Physics Letters 105 (2014).
  * Gomes and Carvalho [2015] L. C. Gomes and A. Carvalho, Phosphorene analogues: Isoelectronic two-dimensional group-iv monochalcogenides with orthorhombic structure, Phys. Rev. B 92, 085406 (2015).
  * Young and Rappe [2012] S. M. Young and A. M. Rappe, First principles calculation of the shift current photovoltaic effect in ferroelectrics, Phys. Rev. Lett. 109, 116601 (2012).
  * Khalaf  _et al._ [2019] E. Khalaf, A. J. Kruchkov, G. Tarnopolsky, and A. Vishwanath, Magic angle hierarchy in twisted graphene multilayers, Physical Review B 100, 085109 (2019).
  * Zhang  _et al._ [2022] Y. Zhang, R. Polski, C. Lewandowski, A. Thomson, Y. Peng, Y. Choi, H. Kim, K. Watanabe, T. Taniguchi, J. Alicea, F. von Oppen, G. Refael, and S. Nadj-Perge, Promotion of superconductivity in magic-angle graphene multilayers, Science 377, 1538–1543 (2022).
  * Park  _et al._ [2022] J. M. Park, Y. Cao, L.-Q. Xia, S. Sun, K. Watanabe, T. Taniguchi, and P. Jarillo-Herrero, Robust superconductivity in magic-angle multilayer graphene family, Nature Materials 21, 877–883 (2022).
  * Tan and Rappe [2019a] L. Z. Tan and A. M. Rappe, Upper limit on shift current generation in extended systems, Phys. Rev. B 100, 085102 (2019a).
  * Xiao  _et al._ [2010] D. Xiao, M.-C. Chang, and Q. Niu, Berry phase effects on electronic properties, Rev. Mod. Phys. 82, 1959 (2010).
  * 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, Phys. Rev. Res. 4, 013164 (2022).
  * Kaplan  _et al._ [2022] D. Kaplan, T. Holder, and B. Yan, Twisted photovoltaics at terahertz frequencies from momentum shift current, Phys. Rev. Research 4, 013209 (2022).
  * Andrei  _et al._ [2021] E. Y. Andrei, D. K. Efetov, P. Jarillo-Herrero, A. H. MacDonald, K. F. Mak, T. Senthil, E. Tutuc, A. Yazdani, and A. F. Young, The marvels of moirématerials, Nature Reviews Materials 6, 201 (2021).
  * Balents  _et al._ [2020] L. Balents, C. R. Dean, D. K. Efetov, and A. F. Young, Superconductivity and strong correlations in moiré flat bands, Nature Physics 16, 725 (2020).
  * Törmä  _et al._ [2022] P. Törmä, S. Peotta, and B. A. Bernevig, Superconductivity, superfluidity and quantum geometry in twisted multilayer systems, Nature Reviews Physics 4, 528 (2022).
  * Rice and Mele [1982] M. J. Rice and E. J. Mele, Elementary excitations of a linearly conjugated diatomic polymer, Phys. Rev. Lett. 49, 1455 (1982).
  * Törmä [2023] P. Törmä, Essay: Where can quantum geometry lead us?, Phys. Rev. Lett. 131, 240001 (2023).
  * Verma  _et al._ [2021] N. Verma, T. Hazra, and M. Randeria, Optical spectral weight, phase stiffness, and t c bounds for trivial and topological flat band superconductors, Proceedings of the National Academy of Sciences 118, e2106744118 (2021).
  * Herzog-Arbeitman  _et al._ [2022] J. Herzog-Arbeitman, V. Peri, F. Schindler, S. D. Huber, and B. A. Bernevig, Superfluid weight bounds from symmetry and quantum geometry in flat bands, Phys. Rev. Lett. 128, 087002 (2022).
  * Xie  _et al._ [2020a] F. Xie, Z. Song, B. Lian, and B. A. Bernevig, Topology-bounded superfluid weight in twisted bilayer graphene, Phys. Rev. Lett. 124, 167002 (2020a).
  * Provost and Vallee [1980] J. P. Provost and G. Vallee, Riemannian structure on manifolds of quantum states, Communications in Mathematical Physics 76, 289–301 (1980).
  * Peotta and Törmä [2015] S. Peotta and P. Törmä, Superfluidity in topologically nontrivial flat bands, Nature Communications 6, 8944 (2015).
  * Julku  _et al._ [2020] A. Julku, T. J. Peltonen, L. Liang, T. T. Heikkilä, and P. Törmä, Superfluid weight and berezinskii-kosterlitz-thouless transition temperature of twisted bilayer graphene, Phys. Rev. B 101, 060505 (2020).
  * Tan and Rappe [2019b] L. Z. Tan and A. M. Rappe, Effect of wavefunction delocalization on shift current generation, Journal of Physics: Condensed Matter 31, 084002 (2019b).
  * Koshino  _et al._ [2018] M. Koshino, N. F. Yuan, T. Koretsune, M. Ochi, K. Kuroki, and L. Fu, Maximally localized wannier orbitals and the extended hubbard model for twisted bilayer graphene, Physical Review X 8, 031087 (2018).
  * Cea and Guinea [2021] T. Cea and F. Guinea, Coulomb interaction, phonons, and superconductivity in twisted bilayer graphene, arXiv:2103.01815 (2021).
  * Guinea and Walet [2018] F. Guinea and N. R. Walet, Electrostatic effects, band distortions, and superconductivity in twisted graphene bilayers, Proceedings of the National Academy of Sciences 115, 13174 (2018).
  * Pantaleón  _et al._ [2022] P. A. Pantaleón, V. o. T. Phong, G. G. Naumis, and F. Guinea, Interaction-enhanced topological hall effects in strained twisted bilayer graphene, Phys. Rev. B 106, L161101 (2022).
  * Goodwin  _et al._ [2020] Z. A. Goodwin, V. Vitale, X. Liang, A. A. Mostofi, and J. Lischner, Hartree theory calculations of quasiparticle properties in twisted bilayer graphene, Electronic Structure 2, 034001 (2020).
  * Rademaker  _et al._ [2019] L. Rademaker, D. A. Abanin, and P. Mellado, Charge smoothening and band flattening due to hartree corrections in twisted bilayer graphene, Phys. Rev. B 100, 205114 (2019).
  * Klebl  _et al._ [2021] L. Klebl, Z. A. H. Goodwin, A. A. Mostofi, D. M. Kennes, and J. Lischner, Importance of long-ranged electron-electron interactions for the magnetic phase diagram of twisted bilayer graphene, Phys. Rev. B 103, 195127 (2021).
  * Rademaker and Mellado [2018] L. Rademaker and P. Mellado, Charge-transfer insulation in twisted bilayer graphene, Phys. Rev. B 98, 235158 (2018).
  * Calderón and Bascones [2020] M. J. Calderón and E. Bascones, Interactions in the 8-orbital model for twisted bilayer graphene, Phys. Rev. B 102, 155149 (2020).
  * Xie and MacDonald [2020] M. Xie and A. H. MacDonald, Nature of the correlated insulator states in twisted bilayer graphene, Phys. Rev. Lett. 124, 097601 (2020).
  * Bernevig  _et al._ [2021] B. A. Bernevig, B. Lian, A. Cowsik, F. Xie, N. Regnault, and Z.-D. Song, Twisted bilayer graphene. v. exact analytic many-body excitations in coulomb hamiltonians: Charge gap, goldstone modes, and absence of cooper pairing, Phys. Rev. B 103, 205415 (2021).
  * Kang and Vafek [2019] J. Kang and O. Vafek, Strong coupling phases of partially filled twisted bilayer graphene narrow bands, Phys. Rev. Lett. 122, 246401 (2019).
  * Bultinck  _et al._ [2020] N. Bultinck, E. Khalaf, S. Liu, S. Chatterjee, A. Vishwanath, and M. P. Zaletel, Ground state and hidden symmetry of magic-angle graphene at even integer filling, Phys. Rev. X 10, 031034 (2020).
  * Lian  _et al._ [2021] B. Lian, Z.-D. Song, N. Regnault, D. K. Efetov, A. Yazdani, and B. A. Bernevig, Twisted bilayer graphene. iv. exact insulator ground states and phase diagram, Phys. Rev. B 103, 205414 (2021).
  * Kwan  _et al._ [2021] Y. H. Kwan, G. Wagner, T. Soejima, M. P. Zaletel, S. H. Simon, S. A. Parameswaran, and N. Bultinck, Kekulé spiral order at all nonzero integer fillings in twisted bilayer graphene, Phys. Rev. X 11, 041063 (2021).
  * Kolář  _et al._ [2023] K. c. v. Kolář, Y. Zhang, S. Nadj-Perge, F. von Oppen, and C. Lewandowski, Electrostatic fate of $n$-layer moiré graphene, Phys. Rev. B 108, 195148 (2023).
  * Cao  _et al._ [2018a] Y. Cao, V. Fatemi, A. Demir, S. Fang, S. L. Tomarken, J. Y. Luo, J. D. Sanchez-Yamagishi, K. Watanabe, T. Taniguchi, E. Kaxiras, R. C. Ashoori, and P. Jarillo-Herrero, Correlated insulator behaviour at half-filling in magic-angle graphene superlattices, Nature 556, 80 EP (2018a).
  * Cao  _et al._ [2018b] Y. Cao, V. Fatemi, S. Fang, K. Watanabe, T. Taniguchi, E. Kaxiras, and P. Jarillo-Herrero, Unconventional superconductivity in magic-angle graphene superlattices, Nature 556, 43 (2018b).
  * Yankowitz  _et al._ [2019] M. Yankowitz, S. Chen, H. Polshyn, Y. Zhang, K. Watanabe, T. Taniguchi, D. Graf, A. F. Young, and C. R. Dean, Tuning superconductivity in twisted bilayer graphene, Science 363, 1059 (2019).
  * Lu  _et al._ [2019] X. Lu, P. Stepanov, W. Yang, M. Xie, M. A. Aamir, I. Das, C. Urgell, K. Watanabe, T. Taniguchi, G. Zhang, A. Bachtold, A. H. MacDonald, and D. K. Efetov, Superconductors, orbital magnets and correlated states in magic-angle bilayer graphene, Nature 574, 653 (2019).
  * Oh  _et al._ [2021] M. Oh, K. P. Nuckolls, D. Wong, R. L. Lee, X. Liu, K. Watanabe, T. Taniguchi, and A. Yazdani, Evidence for unconventional superconductivity in twisted bilayer graphene, Nature 600, 240–245 (2021).
  * Sharpe  _et al._ [2019] A. L. Sharpe, E. J. Fox, A. W. Barnard, J. Finney, K. Watanabe, T. Taniguchi, M. A. Kastner, and D. Goldhaber-Gordon, Emergent ferromagnetism near three-quarters filling in twisted bilayer graphene, Science 365, 605 (2019).
  * Serlin  _et al._ [2020] M. Serlin, C. L. Tschirhart, H. Polshyn, Y. Zhang, J. Zhu, K. Watanabe, T. Taniguchi, L. Balents, and A. F. Young, Intrinsic quantized anomalous hall effect in a moiré heterostructure, Science 367, 900–903 (2020).
  * Nuckolls  _et al._ [2020] K. P. Nuckolls, M. Oh, D. Wong, B. Lian, K. Watanabe, T. Taniguchi, B. A. Bernevig, and A. Yazdani, Strongly correlated chern insulators in magic-angle twisted bilayer graphene, Nature 588, 610–615 (2020).
  * Arora  _et al._ [2020] H. S. Arora, R. Polski, Y. Zhang, A. Thomson, Y. Choi, H. Kim, Z. Lin, I. Z. Wilson, X. Xu, J.-H. Chu, K. Watanabe, T. Taniguchi, J. Alicea, and S. Nadj-Perge, Superconductivity in metallic twisted bilayer graphene stabilized by wse2, Nature 583, 379 (2020).
  * Choi  _et al._ [2019] Y. Choi, J. Kemmer, Y. Peng, A. Thomson, H. Arora, R. Polski, Y. Zhang, H. Ren, J. Alicea, G. Refael, F. von Oppen, K. Watanabe, T. Taniguchi, and S. Nadj-Perge, Electronic correlations in twisted bilayer graphene near the magic angle, Nature Physics 15, 1174 (2019).
  * Choi  _et al._ [2021] Y. Choi, H. Kim, C. Lewandowski, Y. Peng, A. Thomson, R. Polski, Y. Zhang, K. Watanabe, T. Taniguchi, J. Alicea, and S. Nadj-Perge, Interaction-driven band flattening and correlated phases in twisted bilayer graphene, Nature Physics 17, 1375–1381 (2021).
  * Kim  _et al._ [2023] H. Kim, Y. Choi, E. Lantagne-Hurtubise, C. Lewandowski, A. Thomson, L. Kong, H. Zhou, E. Baum, Y. Zhang, L. Holleis, K. Watanabe, T. Taniguchi, A. F. Young, J. Alicea, and S. Nadj-Perge, Imaging inter-valley coherent order in magic-angle twisted trilayer graphene, Nature 623, 942–948 (2023).
  * Nuckolls  _et al._ [2023] K. P. Nuckolls, R. L. Lee, M. Oh, D. Wong, T. Soejima, J. P. Hong, D. Călugăru, J. Herzog-Arbeitman, B. A. Bernevig, K. Watanabe, T. Taniguchi, N. Regnault, M. P. Zaletel, and A. Yazdani, Quantum textures of the many-body wavefunctions in magic-angle graphene, Nature 620, 525–532 (2023).
  * Ahn  _et al._ [2020] J. Ahn, G.-Y. Guo, and N. Nagaosa, Low-frequency divergence and quantum geometry of the bulk photovoltaic effect in topological semimetals, Phys. Rev. X 10, 041041 (2020).
  * Ahn  _et al._ [2021] J. Ahn, G.-Y. Guo, N. Nagaosa, and A. Vishwanath, Riemannian geometry of resonant optical responses, Nature Physics 18, 290–295 (2021).
  * Huhtinen  _et al._ [2022] K.-E. Huhtinen, J. Herzog-Arbeitman, A. Chew, B. A. Bernevig, and P. Törmä, Revisiting flat band superconductivity: Dependence on minimal quantum metric and band touchings, Phys. Rev. B 106, 014518 (2022).
  * Lewandowski and Levitov [2019] C. Lewandowski and L. Levitov, Intrinsically undamped plasmon modes in narrow electron bands, Proceedings of the National Academy of Sciences 116, 20869 (2019).
  * Arora  _et al._ [2022] A. Arora, M. S. Rudner, and J. C. Song, Quantum plasmonic nonreciprocity in parity-violating magnets, Nano Letters 22, 9351 (2022).
  * Dutta  _et al._ [2023] D. Dutta, A. Chakraborty, and A. Agarwal, Intrinsic nonreciprocal bulk plasmons in noncentrosymmetric magnetic systems, Phys. Rev. B 107, 165404 (2023).
  * Kitamura  _et al._ [2020] S. Kitamura, N. Nagaosa, and T. Morimoto, Nonreciprocal landau–zener tunneling, Communications Physics 3, 10.1038/s42005-020-0328-0 (2020).
  * Xie  _et al._ [2020b] F. Xie, Z. Song, B. Lian, and B. A. Bernevig, Topology-bounded superfluid weight in twisted bilayer graphene, Phys. Rev. Lett. 124, 167002 (2020b).
  * Su  _et al._ [1979] W. P. Su, J. R. Schrieffer, and A. J. Heeger, Solitons in polyacetylene, Phys. Rev. Lett. 42, 1698 (1979).
  * Zhang  _et al._ [2021] Y. Zhang, R. Polski, C. Lewandowski, A. Thomson, Y. Peng, Y. Choi, H. Kim, K. Watanabe, T. Taniguchi, J. Alicea, _et al._ , Ascendance of superconductivity in magic-angle graphene multilayers, arXiv:2112.09270 (2021).

## Supplementary Materials

###  .1 Deriving shift current integrand expression

For completeness, we include the derivation of the shift current integrand
used in the main text. First derived in [2], the expression of shift current
can be written in terms of dipole matrix element $r_{mn}^{\alpha}$ and its
generalized derivative $r_{mn;\beta}^{\alpha}$ is given by

$\displaystyle\sigma^{\alpha\alpha}_{\beta}(0;\omega,-\omega)=\frac{-i\pi
e^{3}}{\hbar^{2}}\int_{BZ}\sum_{mn}f_{mn}r_{mn;\beta}^{\alpha}r_{mn}^{\alpha}\delta(\omega_{mn}-\omega),$
(S.1)

where

$\displaystyle r_{mn}^{\alpha}$ $\displaystyle=A_{mn}^{\alpha}\qquad m\neq
n\text{ and $0$ otherwise}$ (S.2) $\displaystyle r_{mn;\beta}^{\alpha}$
$\displaystyle=\frac{\partial r_{mn}^{\alpha}}{\partial
k^{\beta}}-i(A_{mm}^{\beta}-A_{nn}^{\beta})r_{mn}^{\alpha}.$ (S.3)

The generalized sum rule for generalized derivative as provided in [10] is

$\displaystyle
r_{mn;\beta}^{\alpha}=-\frac{1}{i\varepsilon_{mn}}\left[\frac{h_{mn}^{\alpha}\Delta_{mn}^{\beta}+h_{mn}^{\beta}\Delta_{mn}^{\alpha}}{\varepsilon_{mn}}-w_{mn}^{\alpha\beta}+\sum_{l\neq
m,n}\left(\frac{h_{ml}^{\alpha}h_{ln}^{\beta}}{\varepsilon_{ln}}-\frac{h_{ml}^{\beta}h_{ln}^{\alpha}}{\varepsilon_{ml}}\right)\right],$
(S.4)

where $h_{mn}^{\alpha}=\matrixelement{m}{\partial_{k_{\alpha}}H}{n}$,
$w_{mn}^{\alpha\beta}=\matrixelement{m}{\partial_{k_{\alpha}}\partial_{k_{\beta}}H}{n}$.
Using $r_{nm}^{\alpha}=-ih_{nm}^{\alpha}/\varepsilon_{nm}$, the shift current
integrand as defined in the main text is then

$\displaystyle R_{mn}^{\alpha\alpha\beta}$
$\displaystyle=\text{Im}\left[{r_{mn;\beta}^{\alpha}r_{nm}^{\alpha}}\right]$
(S.5)
$\displaystyle=\frac{1}{\varepsilon_{mn}^{2}}\text{Im}\left[-\frac{h_{nm}^{\alpha}h_{mn}^{\alpha}\Delta_{mn}^{\beta}+h_{nm}^{\alpha}h_{mn}^{\beta}\Delta_{mn}^{\alpha}}{\varepsilon_{mn}}+w_{mn}^{\alpha\beta}h_{mn}^{\alpha}\right]$
(S.6) $\displaystyle+\frac{1}{\varepsilon_{mn}^{2}}\text{Im}\left[\sum_{l\neq
mn}\left(\frac{h_{nm}^{\alpha}h_{ml}^{\beta}h_{ln}^{\alpha}}{\varepsilon_{ml}}-\frac{h_{nm}^{\alpha}h_{ml}^{\alpha}h_{ln}^{\beta}}{\varepsilon_{ln}}\right)\right].$
(S.7)

Noting that $h_{nm}^{\alpha}h_{mn}^{\alpha}=|h_{mn}^{\alpha}|^{2}$ and
$\text{Im}[h_{mn}^{\alpha}h_{nm}^{\beta}]=-\text{Im}[h_{nm}^{\alpha}h_{mn}^{\beta}]$,
we recover the expression in the main text:

$\displaystyle R_{mn}^{\alpha\alpha\beta}$
$\displaystyle=\frac{1}{\varepsilon_{mn}^{2}}\text{Im}\left[\frac{h_{mn}^{\alpha}h_{nm}^{\beta}\Delta_{mn}^{\alpha}}{\varepsilon_{mn}}+w_{mn}^{\alpha\beta}h_{mn}^{\alpha}\right]+\frac{1}{\varepsilon_{mn}^{2}}\text{Im}\left[\sum_{l\neq
mn}\left(\frac{h_{nm}^{\alpha}h_{ml}^{\beta}h_{ln}^{\alpha}}{\varepsilon_{ml}}-\frac{h_{nm}^{\alpha}h_{ml}^{\alpha}h_{ln}^{\beta}}{\varepsilon_{ln}}\right)\right].$
(S.8)

###  .2 Fubini-Study metric

The quantum geometric tensor provide a tool for us to describe the geometry of
eigenstate. In the context of bandstructure, the quantum geometric tensor[39,
35] for a single band $m$ is

$\displaystyle G^{m}_{\mu\nu}(k)$
$\displaystyle=\bra{\partial_{\mu}m(k)}(1-\ket{m(k)}\bra{m(k)})\ket{\partial_{\nu}m(k)}$
(S.9) $\displaystyle=\sum_{m\neq
n}\bra{\partial_{\mu}m(k)}\ket{n(k)}\bra{n(k)}\ket{\partial_{\nu}m(k)}$ (S.10)
$\displaystyle=\sum_{m\neq
n}\frac{\matrixelement{m}{\partial_{\mu}H}{n}\matrixelement{n}{\partial_{\nu}H}{m}}{(E_{m}-E_{n})^{2}}$
(S.11)

Under this definition, we find

$\displaystyle
G^{m}_{\mu\nu}(k)=\text{Re}[G^{m}_{\mu\nu}(k)]+i\text{Im}[G^{m}_{\mu\nu}(k)]=g^{m}_{\mu\nu}(k)-\frac{i}{2}F^{m}_{\mu\nu},$
(S.12)

with the real part being the symmetric Fubini-Study (FS) metric and the
imaginary part being the antisymmetric Berry curvature[39, 35], given by

$\displaystyle g_{\mu\nu}^{m}$ $\displaystyle=\frac{1}{2}\sum_{m\neq
n}\frac{\matrixelement{m}{\partial_{\mu}H}{n}\matrixelement{n}{\partial_{\nu}H}{m}+\matrixelement{m}{\partial_{\nu}H}{n}\matrixelement{n}{\partial_{\mu}H}{m}}{(E_{m}-E_{n})^{2}},$
(S.13) $\displaystyle F_{\mu\nu}^{m}$ $\displaystyle=i\sum_{m\neq
n}\frac{\matrixelement{m}{\partial_{\mu}H}{n}\matrixelement{n}{\partial_{\nu}H}{m}-\matrixelement{m}{\partial_{\nu}H}{n}\matrixelement{n}{\partial_{\mu}H}{m}}{(E_{m}-E_{n})^{2}}$
(S.14)

As pointed out in [79], the FS metric characterizes the spread of Wannier
wavefunction. Following the definition of Ref.[79], we define the localization
of Wannier wavefunciton as

$\displaystyle
F=\sum_{n}\left[\matrixelement{0n}{\mathbf{{\bf{\widehat{r}}}}^{2}}{0n}-|\matrixelement{0n}{\mathbf{{\bf{\widehat{r}}}}}{0n}|^{2}\right],$
(S.15)

where $\mathbf{{\bf{\widehat{r}}}}$ is the position operator and
$\ket{\mathbf{R}n}$ is the Wannier state given by
$\ket{\mathbf{R}n}=\frac{1}{\sqrt{N}}\sum_{\mathbf{k}}e^{-i\mathbf{k}\cdot\mathbf{R}}\ket{\psi_{n\mathbf{k}}}$.
Then, FS metric provides a lower bound on $F$ given by

$\displaystyle F\geq\frac{\Omega_{c}}{(2\pi)^{2}}\int d^{2}k\tr g(k),$ (S.16)

with $\Omega_{c}$ being the area of the unit cell and
$g_{\mu\nu}(k)=\sum_{m}g^{m}_{\mu\nu}(k)$.

The FS metric defined above refers to a single band $m$ and sums over
contribution from all intermediate band $n$. Without summing over $n$, we
define

$\displaystyle Q^{mn}_{\mu\nu}(k)$
$\displaystyle=\bra{\partial_{\mu}m(k)}\ket{n(k)}\bra{n(k)}\ket{\partial_{\nu}m(k)}$
(S.17)
$\displaystyle=\frac{\matrixelement{m}{\partial_{\mu}H}{n}\matrixelement{n}{\partial_{\nu}H}{m}}{(E_{m}-E_{n})^{2}},$
(S.18)

we can consequently define interband FS metric and interband Berry curvature
as

$\displaystyle g^{mn}_{\mu\nu}(k)=\text{Re}[Q_{\mu\nu}^{mn}(k)],\qquad
F^{mn}_{\mu\nu}=-2\text{Im}[G^{mn}_{\mu\nu}(k)].$ (S.19)

In terms of projector notations, we can write

$\displaystyle
G_{\mu\nu}^{mn}=\Tr[P_{n}\partial_{\mu}P_{m}\partial_{\nu}P_{m}]=\Tr[P_{m}\partial_{\mu}P_{nm}\partial_{\nu}P_{mn}],$
(S.20)

where $P_{m}=\ket{m}\bra{n}$, and $P_{mn}=\ket{m}\bra{n}$.

In 1D, there is no interband Berry curvature since $G_{xx}^{mn}$ is real.
Furthermore, the transition probability amplitude can be identified with
interband FS metric in the direction of transition:

$\displaystyle|A_{mn}^{\alpha}|^{2}=-\bra{m(k)}\ket{\partial_{\alpha}n(k)}\bra{n(k)}\ket{\partial_{\alpha}m(k)}=g_{mn}^{\alpha\alpha}.$
(S.21)

This allows us to write the shift current expression of Eq. 2 as

$\displaystyle\sigma^{\beta}_{\alpha\alpha}(\omega)=\frac{\pi
e^{3}}{\hbar^{2}}\sum_{m,n}\int
d^{2}\mathbf{k}f_{mn}S^{\beta\alpha}_{mn}g_{mn}^{\alpha\alpha}\delta(\omega-\varepsilon_{mn}),$
(S.22)

where we establish a direct connection between shift current response and FS
metric. Since $\sum_{mn\alpha}\int d\mathbf{k}g_{mn}^{\alpha\alpha}$ gives a
lower bound for the spread of Wannier wavefunction, this confirms many of the
experimental evidence[1, 8] that materials with delocalized wavefunction tend
to exhibit large shift current response. In addition, viewing
$g_{mn}^{\alpha\alpha}$ as the $\alpha$ component of the metric, we can view
shift current conductivity tensor as the sum of shift vector at each k-state
over the BZ weighted by the interband metric along the direction of
transition.

###  .3 Resonance of the Quantum Geometric Tensor in Parameter Space

Near the parameter region where the band gap vanishes at the band edge, we
expect $g^{mn}_{xx}$ to diverge. We describe the rate of divergence by a
Lorentzian function

$\displaystyle L(x)=\frac{2A}{\pi}\frac{\Gamma}{4(x-x_{0})^{2}+\Gamma^{2}}+c,$
(S.23)

with constant $A,x_{0},\Gamma,c$. For the 2RM model considered in the main
text, we fix parameter $(t,\delta_{1},\Delta_{1},\Delta_{2},\varepsilon)$ and
vary $\delta_{2}$ to optimize for shift current enhancement via virtual
transitions, and we note that the divergence of FS metric at the band edge as
$\delta_{2}\rightarrow\delta_{2,0}$ is well-characterized by a Lorentzian
function, where $\delta_{2,0}$ is where the gap is minimized.

Figure S1: The FS metric $g^{mn}_{xx}(k)$ at the band edge $ka=-\pi$ for
$m=3$, $n=2,1,0$, and the fit using a Lorentzian function, with fitting curve
in red and exact calculation in blue dots. Here we used the same parameters as
in the main text, namely
$(t,\delta_{1},\Delta_{1},\Delta_{2},\varepsilon)=(1,0.8,0.7,0.6,0.1)$.

We plot the component of $g^{mn}$ at the band edge for the upper conduction
band $m=3$ and $n=0,1,2$ in Fig. S1a,b,c. For $g^{32}(k)$, the Lorentzian
function well-characterizes the divergence at the band edge as $\delta_{2,0}$
is approached. For $g^{31}(k)$ and $g^{30}(k)$, the tails show deviations from
Lorentzian behavior, but near the resonance pole, they can be described by the
Lorentzian function. It is this divergence structure that gives rise to
enhanced $g_{3}(k)=\sum_{n\neq 3}g^{3n}(k)$ near $\delta_{2,0}$.

This evokes the resonance phenomena but in parameter space. As we vary the
hopping parameters, there exists a critical hopping strength such that the FS
metric diverges at the band edge. This enhances the sum of FS over all
k-states, leading to Wannier wavefunction being “resonantly” more delocalized.
Consequently, we observe an enhancement in shift current response via virtual
transitions.

Now, we demonstrate this connection more precisely by explicitly putting
$g^{mn}$ into the form of a Lorentzian. The exact energy eigenvalues of 2RM
model is given by

$\displaystyle a$
$\displaystyle=\Delta_{1}^{2}+\Delta_{2}^{2}+\delta_{1}^{2}\sin^{2}{\left(\frac{k}{2}\right)}+\delta_{2}^{2}\sin^{2}{\left(\frac{k}{2}\right)}+2t^{2}\cos^{2}{\left(\frac{k}{2}\right)},$
(S.24) $\displaystyle b$
$\displaystyle=4t^{2}\cos^{2}{\left(\frac{k}{2}\right)}+\left(\Delta_{1}-\Delta_{2}\right)^{2}+\left(\delta_{1}-\delta_{2}\right)^{2}\sin^{2}{\left(\frac{k}{2}\right)},$
(S.25) $\displaystyle c$
$\displaystyle=\left(\Delta_{1}-\Delta_{2}\right)\left(\Delta_{1}+\Delta_{2}\right)+\left(\delta_{1}-\delta_{2}\right)\left(\delta_{1}+\delta_{2}\right)\sin^{2}{\left(\frac{k}{2}\right)},$
(S.26) $\displaystyle E$
$\displaystyle=\pm\frac{\sqrt{2}}{2}\sqrt{a+2\varepsilon^{2}\pm\sqrt{4\varepsilon^{2}b+c^{2}}}.$
(S.27)

Recall that the FS metric tensor defined in Appendix .2 in 1D is given by

$\displaystyle
g^{mn}(k)=\frac{\matrixelement{m}{\partial_{k}H}{n}\matrixelement{n}{\partial_{k}H}{m}}{(E_{m}-E_{n})^{2}}.$
(S.28)

Using $m=3$, $n=2$ as an example and focus on the band edge at $k=\pm\pi$, we
put the above FS in the form of Lorentzian explicitly. Noting that there are
no poles in the numerator, so we focus on the denominator term, which gives

$\displaystyle(E_{3}-E_{2})^{2}=a+2\varepsilon^{2}-\sqrt{(a+2\varepsilon^{2})^{2}-4\varepsilon^{2}b+c^{2}}\approx\frac{1}{2}\frac{4\varepsilon^{2}b+c^{2}}{(a+2\varepsilon^{2})},$
(S.29)

where we note that $4\varepsilon^{2}b+c^{2}\ll(a+2\varepsilon^{2})^{2}$ near
the pole. The pole occurs in $4\varepsilon^{2}b+c^{2}$. We expand $\delta_{2}$
around the maximum of $g^{32}(k)$, denoted by $\delta_{2,0}$, to second order,
and we find

$\displaystyle
4\varepsilon^{2}b+c^{2}=\tilde{a}\left(\delta_{2}-\delta_{2,0}+\frac{\tilde{b}}{2\tilde{a}}\right)^{2}-\frac{\tilde{b}^{2}}{4\tilde{a}}+\tilde{c}+\mathcal{O}((\delta_{2}-\delta_{2,0})^{2}),$
(S.30)

with

$\displaystyle\tilde{a}$
$\displaystyle=-2(\delta_{1}^{2}+\Delta_{1}^{2}-\Delta_{2}^{2}-3\delta_{2,0}^{2}-2\varepsilon^{2})$
(S.31) $\displaystyle\tilde{b}$
$\displaystyle=-4(\delta_{2,0}(\delta_{1}^{2}+\Delta_{1}^{2}-\Delta_{2}^{2}-\delta_{2,0}^{2})+2(\delta_{1}-\delta_{2,0})\varepsilon^{2})$
(S.32) $\displaystyle\tilde{c}$
$\displaystyle=(\delta_{1}^{2}+\Delta_{1}^{2}-\Delta_{2}^{2}-\delta_{2,0}^{2})^{2}+4((\Delta_{1}-\Delta_{2})^{2}+(\delta_{1}-\delta_{2,0})^{2})\varepsilon^{2}.$
(S.33)

By the assumption that $\delta_{2,0}$ is the pole, we require
$\frac{\tilde{b}}{2\tilde{a}}=0$. Alternatively, the value of $\delta_{2,0}$
is obtained by solving for $\frac{\tilde{b}}{2\tilde{a}}=0$. The FS metric
$g^{32}$ at $k=\pm\pi$ as a function of hopping strength $\delta_{2}$ can be
written in the form of a Lorenzian

$\displaystyle
g^{32}(\delta_{2})=\frac{2A}{\pi}\frac{\Gamma}{4(\delta_{2}-\delta_{2,0})^{2}+\Gamma^{2}},$
(S.34)

with the width and amplitude given by

$\displaystyle\Gamma=2\sqrt{\frac{\tilde{c}}{\tilde{a}}},\qquad
A=\frac{4\pi(a+2\varepsilon^{2})|\matrixelement{3}{\partial_{k}H}{2}|^{2}}{\tilde{a}\Gamma}.$
(S.35)

A similar manipulation follows for the case of $n=1$ and $n=0$. The $\Gamma$
terms acts like a damping term that restricts the divergences of FS at the
band edge. This $\Gamma$ damping term usually occurs in systems with many
parameters, and may be absent in models with fewer parameters. Consider the
Su–Schrieffer–Heeger (SSH) model[80], with Hamiltonian given by

$\displaystyle H(k)=h_{x}\sigma_{x}+h_{y}\sigma_{y},\qquad
h_{x}=v+w\cos(k),\qquad h_{y}=w\sin(k),$ (S.36)

where $v$ is the intralayer hopping amplitude and $w$ is the interlayer
hopping amplitude. Fixing interlayer hopping $w$ and vary intralayer hopping
$v$, the FS metric at $k=\pm\pi$ is given by

$\displaystyle g(v)$ $\displaystyle=\frac{w^{2}}{4(v-w)^{2}},$ (S.37)

where it diverges without damping at $v=w$, the point when energy gap closes.
Adding a mass term $\Delta$ to break the chiral symmetry, we obtain a single-
layer Rice-Mele Model, and the damping in FS metric becomes $4\Delta$.

###  .4 Analytical Results of 2RM model

In the main text, we established that it is the enhancement in virtual
transitions near the critical parameter region that leads to enhanced shift
current response. Inspecting Eq. 3 shows us the two factors that control of
the magnitude of the system are energy gap $\varepsilon_{mn}$ and matrix
element $h_{mn}=\matrixelement{m}{\partial_{k}H}{n}$. In particular, we expect
$\varepsilon_{32}$ to dominate in the denominator when the gap is small, and
we will show that it is divergence in $\text{Im}\\{h_{20}\\}$ that dominates
in the numerator.

The Hamiltonian of 2RM model takes the form $H=H_{0}+V$, with

$\displaystyle
H_{0}=\begin{pmatrix}d^{(1)}_{z}&d^{(1)}_{x}-id^{(1)}_{y}&0&0\\\
d^{(1)}_{x}+id^{(1)}_{y}&-d^{(1)}_{z}&0&0\\\
0&0&d^{(2)}_{z}&d^{(2)}_{x}-id^{(2)}_{y}\\\
0&0&d^{(2)}_{x}+id^{(2)}_{y}&-d^{(2)}_{z}\end{pmatrix},\quad
V=\begin{pmatrix}0&0&0&\varepsilon\\\ 0&0&\varepsilon&0\\\
0&\varepsilon&0&0\\\ \varepsilon&0&0&0\end{pmatrix},$ (S.38)

where $H$ is in the layer basis $\\{\ket{1A},\ket{1B},\ket{2A},\ket{2B}\\}$,
and superscript $(1/2)$ denotes the first/second RM model. Exactly solving for
the eigenstate and computing matrix element will result in complicated
expressions that shed little intuition. Instead, we will use perturbation
theory and expand near the band edge. A direct application of non-degenerate
perturbation theory leads to inaccurate results near the band edge; hence, we
first perform a basis transformation to the unperturbed eigenstate
$\\{\ket{\psi_{1,A}},\ket{\psi_{2,A}},\ket{\psi_{1,B}},\ket{\psi_{2,B}}\\}$.
The Hamiltonian $\tilde{H}=\tilde{H_{0}}+\tilde{V}$ becomes

$\displaystyle\tilde{H_{0}}=\begin{pmatrix}-\varepsilon_{1}(k)&f_{A}(k)&0&0\\\
f_{A}^{*}(k)&-\varepsilon_{2}(k)&0&0\\\ 0&0&\varepsilon_{1}(k)&f_{B}(k)\\\
0&0&f^{*}_{B}(k)&\varepsilon_{2}(k)\end{pmatrix},\quad\tilde{V}=\begin{pmatrix}0&0&0&g_{1}(k)\\\
0&0&g_{2}(k)&0\\\ 0&g_{2}^{*}(k)&0&0\\\ g_{1}^{*}(k)&0&0&0\end{pmatrix},$
(S.39)

where $\varepsilon_{1/2}$ is the eigenvalue of the first/second unperturbed RM
model, and

$\displaystyle f_{A}(k)$
$\displaystyle=\matrixelement{\psi_{1,A}}{V}{\psi_{2A}}=\frac{\varepsilon}{2\sqrt{\varepsilon_{1}\varepsilon_{2}}}\left(-e^{i\varphi_{k}^{(2)}}\sqrt{\varepsilon_{1}-d_{z}^{(1)}}\sqrt{\varepsilon_{2}+d_{z}^{(2)}}-e^{-i\varphi_{k}^{(1)}}\sqrt{\varepsilon_{1}+d_{z}^{(1)}}\sqrt{\varepsilon_{2}-d_{z}^{(2)}}\right),$
(S.40) $\displaystyle f_{B}(k)$
$\displaystyle=\matrixelement{\psi_{1,B}}{V}{\psi_{2B}}=\frac{\varepsilon}{2\sqrt{\varepsilon_{1}\varepsilon_{2}}}\left(e^{-i\varphi_{k}^{(1)}}\sqrt{\varepsilon_{1}-d_{z}^{(1)}}\sqrt{\varepsilon_{2}+d_{z}^{(2)}}+e^{i\varphi_{k}^{(2)}}\sqrt{\varepsilon_{1}+d_{z}^{(1)}}\sqrt{\varepsilon_{2}-d_{z}^{(2)}}\right),$
(S.41) $\displaystyle g_{1}(k)$
$\displaystyle=\matrixelement{\psi_{1,A}}{V}{\psi_{2B}}=\frac{\varepsilon}{2\sqrt{\varepsilon_{1}\varepsilon_{2}}}\left(-e^{i\varphi_{k}^{(2)}}\sqrt{\varepsilon_{1}-d_{z}^{(1)}}\sqrt{\varepsilon_{2}-d_{z}^{(2)}}+e^{-i\varphi_{k}^{(1)}}\sqrt{\varepsilon_{1}+d_{z}^{(1)}}\sqrt{\varepsilon_{2}+d_{z}^{(2)}}\right),$
(S.42) $\displaystyle g_{2}(k)$
$\displaystyle=\matrixelement{\psi_{2,A}}{V}{\psi_{1B}}=\frac{\varepsilon}{2\sqrt{\varepsilon_{1}\varepsilon_{2}}}\left(-e^{i\varphi_{k}^{(1)}}\sqrt{\varepsilon_{1}-d_{z}^{(1)}}\sqrt{\varepsilon_{2}-d_{z}^{(2)}}+e^{-i\varphi_{k}^{(2)}}\sqrt{\varepsilon_{1}+d_{z}^{(1)}}\sqrt{\varepsilon_{2}+d_{z}^{(2)}}\right),$
(S.43)

where $\varphi^{(1/2)}_{k}=\arctan(d^{(1/2)}_{y}/d^{(1/2)}_{x})$. Now, we can
apply the non-degenerate perturbation theory safely. In particular, the first
order correction in state is given by

$\displaystyle\ket{n}=\ket{n_{0}}+\sum_{k\neq
n}\frac{\matrixelement{k_{0}}{V}{n_{0}}}{E_{n}-E_{k}}\ket{k_{0}},$ (S.44)

with nonzero matrix element only in

$\displaystyle\matrixelement{1_{0}}{\tilde{V}}{3_{0}}=\frac{1}{2\sqrt{E_{d}^{(A)}E_{d}^{(B)}}}\Big{(}$
$\displaystyle-
g_{1}(k)e^{i\varphi_{k}^{(B)}}\sqrt{E_{d}^{(A)}+d_{3}^{(A)}}\sqrt{E_{d}^{(B)}+d_{3}^{(B)}}-g_{2}(k)e^{-i\varphi_{k}^{(A)}}\sqrt{E_{d}^{(A)}-d_{3}^{(A)}}\sqrt{E_{d}^{(B)}-d_{3}^{(B)}}\Big{)},$
(S.45)
$\displaystyle\matrixelement{2_{0}}{\tilde{V}}{3_{0}}=\frac{1}{2\sqrt{E_{d}^{(A)}E_{d}^{(B)}}}\Big{(}$
$\displaystyle
g_{1}(k)e^{i\varphi_{k}^{(B)}}\sqrt{E_{d}^{(A)}-d_{3}^{(A)}}\sqrt{E_{d}^{(B)}+d_{3}^{(B)}}-g_{2}(k)e^{-i\varphi_{k}^{(A)}}\sqrt{E_{d}^{(A)}+d_{3}^{(A)}}\sqrt{E_{d}^{(B)}-d_{3}^{(B)}}\Big{)},$
(S.46)
$\displaystyle\matrixelement{1_{0}}{\tilde{V}}{4_{0}}=\frac{1}{2\sqrt{E_{d}^{(A)}E_{d}^{(B)}}}\Big{(}$
$\displaystyle-
g_{1}(k)e^{i\varphi_{k}^{(B)}}\sqrt{E_{d}^{(A)}+d_{3}^{(A)}}\sqrt{E_{d}^{(B)}-d_{3}^{(B)}}+g_{2}(k)e^{-i\varphi_{k}^{(A)}}\sqrt{E_{d}^{(A)}-d_{3}^{(A)}}\sqrt{E_{d}^{(B)}+d_{3}^{(B)}}\Big{)},$
(S.47)
$\displaystyle\matrixelement{2_{0}}{\tilde{V}}{4_{0}}=\frac{1}{2\sqrt{E_{d}^{(A)}E_{d}^{(B)}}}\Big{(}$
$\displaystyle
g_{1}(k)e^{i\varphi_{k}^{(B)}}\sqrt{E_{d}^{(A)}-d_{3}^{(A)}}\sqrt{E_{d}^{(B)}-d_{3}^{(B)}}+g_{2}(k)e^{-i\varphi_{k}^{(A)}}\sqrt{E_{d}^{(A)}+d_{3}^{(A)}}\sqrt{E_{d}^{(B)}+d_{3}^{(B)}}\Big{)},$
(S.48)

where

$\displaystyle d_{0}=\frac{\varepsilon_{1}(k)+\varepsilon_{2}(k)}{2}$
$\displaystyle,\quad d_{1}^{(A/B)}=\text{Re}\\{f_{A/B}(k)\\},\quad
d_{2}^{(A/B)}=\text{Im}\\{f_{A/B}(k)\\},\quad
d_{3}=\frac{\varepsilon_{1}(k)-\varepsilon_{2}(k)}{2},$ (S.49) $\displaystyle
E_{d}^{(A/B)}$
$\displaystyle=\sqrt{\left(d_{1}^{(A/B)}\right)^{2}+\left(d_{2}^{(A/B)}\right)^{2}+d_{3}^{2}},\quad\varphi_{k}^{(A/B)}=\arctan(d_{2}^{(A/B)}/d_{1}^{(A/B)})$
(S.50)

In particular, there is no matrix element connecting the same block of
$\tilde{H}_{0}$ where the difference in energy is small, thus avoiding the
divergence in a small energy gap near gap closing. Comparing the perturbation
theory result with exact diagonalization indeed shows the validity of this
approach. The $\tilde{H}_{0}$ now becomes an effective 2-band model, which
shift current response is shown in Fig. 2d. Comparing with tight-binding
2-band model, this effective 2-band model can achieve a large joint density of
states and dispersive band structure simultaneously.

Now, we have the ingredients necessary to compute matrix element $h_{mn}$. For
our consideration, take $m=3,n=0$. At the band edge, comparing the magnitude
of direct and virtual transition contribution to shift current amounts to
comparing $\text{Re}\\{w_{30}\\}$ and
$\text{Re}\\{h_{32}h_{20}\\}=\text{Im}\\{h_{32}\\}\text{Im}\\{h_{20}\\}$ (in
the integrand we are taking the Im of the expression, and
$\text{Re}\\{h_{03}\\}=0$ at band edge). We find

$\displaystyle\text{Re}\\{w_{30}\\}$
$\displaystyle=\matrixelement{3_{0}}{\partial_{k}^{2}\tilde{H}}{0_{0}}+\frac{\matrixelement{3_{0}}{\tilde{V}}{0_{0}}}{E_{30}}\left[\matrixelement{0_{0}}{\partial_{k}^{2}\tilde{H}}{0_{0}}-\matrixelement{3_{0}}{\partial_{k}^{2}\tilde{H}}{3_{0}}\right],$
(S.51) $\displaystyle\text{Im}\\{h_{32}\\}$
$\displaystyle=\frac{\matrixelement{1_{0}}{\partial_{k}\tilde{H}}{2_{0}}\matrixelement{3_{0}}{\tilde{V}}{1_{0}}}{E_{31}}+\frac{\matrixelement{0_{0}}{\partial_{k}\tilde{H}}{2_{0}}\matrixelement{3_{0}}{\tilde{V}}{0_{0}}}{E_{30}}+\frac{\matrixelement{3_{0}}{\partial_{k}\tilde{H}}{1_{0}}\matrixelement{1_{0}}{\tilde{V}}{2_{0}}}{E_{21}}+\frac{\matrixelement{3_{0}}{\partial_{k}\tilde{H}}{0_{0}}\matrixelement{0_{0}}{\tilde{V}}{2_{0}}}{E_{20}},$
(S.52) $\displaystyle\text{Im}\\{h_{20}\\}$
$\displaystyle=\matrixelement{2_{0}}{\partial_{k}\tilde{H}}{0_{0}}+\mathcal{O}(1/E^{2}),$
(S.53)

where we see $\text{Im}\\{h_{20}\\}$ dominate over $\text{Im}\\{h_{32}\\}$
since it is not suppressed by energy. In addition, find
$\text{Im}\\{h_{32}\\}$ is maximized near minimum gap point

$\displaystyle\delta_{2,0}=\pm\sqrt{\delta_{1}^{2}+\Delta_{1}^{2}-\Delta_{2}^{2}}.$
(S.54)

The negative pole is suppressed by other factors, such as direct transition,
and it is only near the positive pole where we see enhancement in the overall
shift current response.

Near the point of minimum energy gap at the band edge, the small
$\varepsilon_{32}$ in the denominator and diverging $\text{Im}\\{h_{20}\\}$ in
the numerator of Eq. 3 combine to generate an enhanced response in virtual
transition. Nevertheless, we remark that the peak shift current response is
not located precisely at gap-closing; rather, it is in close proximity to this
point, due to the vanishing of transition probability amplitude at the gap-
closing point.

###  .5 Total Shift Current Response in 2RM

The presence of multiple bands allows for additional transition pairs to
happen. The total number of transition pairs is the combination of 2 bands out
of all the possible bands, and their contributions need to be summed over when
computing the total shift current response. The total shift current response
for $\nu=0.5$ is given in Fig. S2:

Figure S2: Total shift current response for $\nu=2$ summing up all transition
pairs.

The band index is numbered based on Fig. 2a in the main text. Most of the
contribution from $\nu=0.5$ arises from the transition between band $2$ and
$3$, enhanced due to the small band gap. In addition, there are transitions
from the first RM model (band $0$ to band $3$), giving rise to two peaks,
corresponding to transition mostly from $ka=\pm\pi$ and $ka=0$, respectively.
The transition from the second RM model (band $1$ to band $2$) contributes
only one peak with contribution mostly from $ka=0$, which is due to
suppression in joint density of state. Summing up all transition pairs also
leads to an enhancement in shift current response, compared to the case of
simply stacking two pairs of RM chains together (without band mixing).

###  .6 Band Mixing in TMG

The multilayer Hamiltonian constructed using ATMG model decouples into pairs
of TBG-like bands, but we expect band mixing to occur in actual materials [81,
58]. There are two ways to induce band mixing, either through the application
of an external displacement field or through a self-generated displacement
field.

External displacement field The effect of the external displacement field is
applied through tuning external gate on the device. It modifies the
hamiltonian through the addition of the following term in the layer basis:

$\displaystyle\Delta H=\frac{U}{N-1}\begin{pmatrix}\frac{N-1}{2}&&&&\\\
&\frac{N-1}{2}-1&&&\\\ &&\ddots&&\\\ &&&-\frac{N-1}{2}+1&\\\
&&&&-\frac{N-1}{2}\end{pmatrix}$ (S.55)

where $U$ is the strength of the displacement field, in unit of $meV$. This
leads to band mixing among TBG-like sectors, allowing transition to take
place.

###  .7 T5G Shift Current and Dependence on Twist Angles

The physical twist angle of the TMG system controls the energy gap of nearby
flat bands. At the magic angle where the perfect flat band pair emerges, we
expect a large shift current response due to the vanishing of the energy gap.
In addition, away from the magic angle in multilayer systems, we anticipate a
substantial shift current response stemming from enhanced virtual transition
via nearby flat bands. Here we show the dependence of total shift current
response $\int dw\sigma_{xxy}(\omega)$ on the physical angle $\theta_{5}$ of
T5G. In addition, we compute the integral of FS metric over k-state as a
function of twist angle, as shown in Fig. S3. We indeed observe an enhanced
shift current response near angles where FS metric diverges, corresponding to
a more delocalized Wannier wavefunction.

Figure S3:  T5G FS metric of top flat bands and overall shift current response
and their dependence on physical twist angle. The orange dashed line shows the
angle that leads to max FS, and correspondingly largest overall shift current
measured by metric $M$. The blue dashed line shows the angle that is used in
the main text to illustrate enhancement via virtual transitions. The red line
denotes the magic angle, where we also observe an enhancement in the overall
shift current response, which we attribute to the increased joint density of
states.

In the main text, we illustrate virtual transition enhancement in T5G with
$\theta_{5}=1.32^{\circ}$ since the flat bands belong to distinct TBG-like
sectors, allowing us to better compare the effect of band mixing on
transitions within the same band. As shown in Fig. S3b, this angle is not the
physical angle that leads to maximized response. The ideal angle occurs when
$\theta_{5}=1.28^{\circ}$, where the two flat bands are tightly intertwined.

Figure S4: T5G band structure occurring at $\theta=1.28^{\circ}$ with
$U=20meV$ and the corresponding shift current response at various fillings.

At $\theta=1.28^{\circ}$, the two pairs of TBG have a similar effective twist
angle, but with the top and bottom flat bands flipped with respect to each
other. The slight difference in dispersion relations of the top and bottom
flat bands results in the intertwining of the bands, as shown in Fig. S4.
Instead of studying transitions between single band pairs, observing an
enhancement in shift current necessitates summing over all transition pairs.
Fig. S4 b) shows that virtual transitions among the flat bands indeed lead to
an enhancement response, and we also observe that $\theta=1.28^{\circ}$
generates a larger shift current response than $\theta=1.32^{\circ}$.

To search for the optimal twist angles for a higher number of layers, we would
want the TBG band pairs to have similar effective angles. A more systematic
search can be done by computing the FS metric for the flat bands and
pinpointing the angle that gives rise to the largest $\int d^{2}k\tr g$. We
adopt this procedure and obtained the optimal angles for $N=6,7,8,9$ layers of
twisted graphene, and the result is shown in Fig. S5.

Figure S5: Optimal angle for shift current in TMG a) Integral of trace of FS
metric for the top flat band of T6G, T7G, T8G, T9G, with the orange dashed
line indicating the optimal angle. b) Their corresponding band structure at
optimal angles and c) their conductivity $\sigma_{xxy}(\omega)$ for a wide
range of fillings.

The conductivity includes transitions among flat bands and include transition
to dispersive bands via virtual transitions. We see that peaks in $\int
d^{2}k\tr(g)$ indeed corresponds to the angles where two flat bands are the
closest.

###  .8 T5G Symmetry

The TMG model has $C_{3z}$ symmetry, which restricted
$\sigma^{\beta}_{\alpha\alpha}$ to two independent components. A detailed
argument is discussed in [29] for TBG, so here we simply restate the result
for TMG. We have

$\displaystyle\sigma^{y}_{xx}=-\sigma^{y}_{yy}=\sigma^{y}_{xy}=\sigma^{x}_{yx}\neq
0,$ (S.56)
$\displaystyle\sigma^{x}_{yy}=-\sigma^{x}_{xx}=\sigma^{y}_{yx}=\sigma^{y}_{xy}\neq
0.$ (S.57)

In the main text, we focused on $\sigma^{y}_{xx}$ component.

###  .9 Finite Broadening at Low-Frequency

We point out that the non-zero shift current conductivity at low-frequency is
due to finite Lorentzian broadening. In computing shift current conductivity,
we need to evaluate $\delta(\omega-\varepsilon_{mn})$ in Eq. 2. We choose to
approximate the delta function using

$\displaystyle\delta(x)\approx\frac{\gamma}{\pi(x^{2}+\gamma^{2})},$ (S.58)

where $\gamma$ characterizes the width of the function. We choose $\gamma$ to
be the average energy spacing of the participating bands between nearby
k-points in the Brillouin zone. This large choice of $\gamma$ ensures lack of
spurious (due to a finite mesh size) oscillations in the conductivity from
Fig. 1e at finite frequencies, however it also renders our conductivity finate
as $\omega\to 0$ where the contributions stem from the two closely spaced flat
bands. With higher mesh resolution, conductivity as $\omega\to 0$ vanishes as
required.