Twisted photovoltaics at terahertz frequencies from momentum shift current

The bulk photovoltaic effect (BPVE) converts light into a coherent dc current at zero bias, through what is commonly known as the shift current. This current has previously been attributed to the displacement of the electronic wave function center in real space, when the sample is excited by light. We reveal that materials like twisted bilayer graphene (TBG) with a flatband dispersion are uniquely suited to maximize the BPVE because they lead to an enhanced shift in the momentum space, unlike any previously known shift current mechanism. We identify properties of quantum geometry, which go beyond the quantum geometric tensor, and are unrelated to Berry charges, as the physical origin of the large BPVE we observe in TBG. Our calculations show that TBG with a band gap of several meV exhibits a giant BPVE in a range of 0.2-1 THz, which represents the strongest BPVE reported so far at this frequency in a two-dimensional material and partially persists even a room temperature. Our paper provides a design principle for shift current generation, which applies to a broad range of twisted heterostructures with the potential to overcome the so-called"terahertz gap"in THz sensing.


I. INTRODUCTION
The bulk photovaltaic effect (BPVE) [1,2] refers to the generation of a dc current from a homogeneous solid upon irradiation with light. The so-called shift current [3][4][5] is one of the most important mechanisms for the BPVE which has been appreciated for a long time for its potential in photovoltaic and photonic applications. So far, the BVPE has been studied predominantly for near-infrared and optical frequencies, but it may prove useful for overcoming the so-called 'terahertz gap' [6] in the field of terahertz photonics [7]; that is the perceived lack of robust, tunable and broad-frequency materials for terahertz detection. Since typical excitation energies in twisted bilayer graphene (TBG) are situated in the THz regime [8][9][10][11][12][13][14][15][16][17][18][19][20], we are motivated to study the BPVE effect in this particular material system, as it will involve resonant transitions which occur within this energy window.
Our central insight is depicted schematically in Fig. 1: The shift current appears as the result of a wavefunction displacement in either real space or momentum space. The former interpretation as a real-space shift of wave function centers upon excitation has been proposed since the earliest works [3], and it has recently gained renewed interest in studies of Weyl semimetals [21,22], because the shift is believed to get enhanced by the Berry curvature near Weyl points [23][24][25][26][27][28][29][30]. Very recently, we proposed to instead view the shift current as a result of anomalous quasiparticle acceleration, which is determined by quantum geometric properties [31]. In the latter interpretation, the small effective mass in a Weyl semimetal leads to strong acceleration in the field of light and hence generates large photocurrent, consistent with the established shift current interpretation [26]. However, materials with a flat dispersion represent the other extreme (large quasiparticle mass) compared to Weyl bands (nearly massless quasiparticles). Using the language of a transient po- * binghai.yan@weizmann.ac.il Figure 1. A schematic comparison of the wavefunction in real-space and momentum-space of electronic states at the Fermi surface for three types of gapless materials. For a Weyl semimetal, states are well localized at few points in the momentum space but quite extended in the real space. In contrast, for a flat band all momenta are equally occupied, leading to sharply localized peaks of the wavefunction in real space. The ordinary metal with a large Fermi surface represents the intermediate region between the Weyl semimetal and the flat-band system. If the shift current is viewed as a result of the anomalous acceleration, the Weyl semimetal exhibits a large real-space shift while the flat band displays a major momentum-space shift. larization, it is therefore not obvious how a large shift current can emerge for a flat-band system. In contrast, in the new semiclassical acceleration picture, a state can be displaced (i.e. accelerated) in either real-space and momentum-space, (cf. Fig. 1). Although the shift is indeed small in real-space due to the large quasiparticle mass, for flat bands the acceleration (i.e. the displacement in momentum space) can still be very large.
TBG itself has attracted an immense amount of attention following the recent discovery of correlated insulating arXiv:2101.07539v3 [cond-mat.mes-hall] 23 Mar 2022 and superconducting phases at small twist angle (magic angle) θ ∼ 1 • [32,33]. While the root cause for this exciting behavior is widely accepted to be tied to a highly quenched band structure [34][35][36], opinions differ about the mechanisms associated with the various correlated phases which have since been documented in the system [8][9][10][11][12][13][14][15][16][17][18][19][20]. An important stepping stone towards the understanding of these phases is the characterization of the quantum geometric and dispersive features of the flat bands. Due to inversion symmetry breaking, TBG and similar twisted bilayers are known to exhibit a number of intriguing nonlinear phenomena [37][38][39][40][41][42][43][44] such as the BPVE and nonlinear anomalous Hall effect [45,46]. These nonlinear probes are an important tool for the investigation of quasiparticle properties in flat bands as they are not suppressed by a vanishing Fermi velocity [31].
In this work, we focus on the shift current generation below 10meV in magic angle TBG across the single particle gaps. In a numerical analysis, we find a large photocurrent response for light in the low terahertz regime. This is consistent with the acceleration picture outlined above. To further support such a connection, we express the response in terms of the real-space shift current and momentumspace shift current. Due to the flat-band dispersion, the dominant contribution is constituted by the momentumspace shift, which we explicitly show to be dissimilar to an ordinary dispersive material (e.g., MoS 2 ). We point out that while flat-band materials are expected to be optically active with a large joint density of states (jDOS) [47], this alone could not account for the large response, as it exceeds that of other 2D materials with a large jDOS at the band edge. Distinct from the nonlinear anomalous Hall effect, the shift current of TBG originates from properties of the quantum geometry of the band structure that are unrelated to the Berry curvature dipole (BCD). The magnitude and tunability of the photocurrent, the broadness of the resonances, and the terahertz frequency range in which they are all observed hold promise of the utilization of TBG devices in terahertz technologies. We further explore the effect of heterostrain, as it is typically observed in TBG [48].

A. Shift current theory
In the clean limit, the Bloch states produce an intrinsic dc-current response to linearly polarized light [3]. This nonlinear conductivity σ aa;c (s) is usually formulated as [3,4] where m, n are the band indices, ε mn = ε m − ε n the band energy difference, f mn = f (ε m ) − f (ε n ) the Fermi distribution function difference, and r a mn ≡ m|r a |n the dipole transition matrix element, integration over the Brillouin zone is taken as k = 1 (2π) 2 dk x dk y , while a, c represent the light field and current directions, respectively. S mn is the so-called shift vector, where (r c mm −r c nn ) is the shift of the wave function centers, which is gauge dependent. The second term in Eq. 2 is the phase derivative of r c mn , which ensures gauge-invariance for the shift vector. S mn has then been interpreted as the real-space shift of the mass center of a quasiparticle upon excitation from band m to n [3]. As we will show now, the phase term in (2) gives rise to a shift in momentum space. Therefore, the shift current is actually the result of a shift of the quasiparticle excitation in both real space and momentum space, which is naturally captured using the language of the anomalous quasiparticle acceleration [31]. To this end, we express the nonlinear conductivity in the form of geometric properties of the band structure, where we use the symmetrized derivative λ ab nm = 1 2 ∂ ka r b nm + ∂ k b r a nm . Additionally, we defined the inter-band Berry curvature, Ω ac nm = i l =(n,m) (r a nl r c lm −r c nl r a lm ). For the shift current, this Berry curvature is weighted by the transition element δ(ω ± ε nm ), arising from optical selection rules. The derivative of this object gives the interband BCD that appears in Eq. (6) (the derivative does not act on the delta function). The two pieces R aac shift and K aac shift correspond to the contributions from the first term and the second term, respectively, of the shift vector in Eq. (2). The second term in K aac shift is the Berry curvature dipole, and the first term is related to the other geometric quantity, the quantum metric, which characterizes the quantum distance between two states. Thus, K aac shift originates directly from the quantum geometry of the wave function in k-space, and it represents the anomalous acceleration in momentum space [31]. λ ab nm encodes the skewness of the acceleration: The Berry curvature dipole involves derivatives of the type ∂ c (r a mn r b nm −r b mn r a nm ), which does not cover the entire motion in momentum space. The remaining terms involving λ ab in Eq. (5) can be connected to skew-symmetric derivatives of the structure (∂ c r a mn )r b nm − r a mn (∂ c r b nm ) + (∂ c r b mn )r a nm − r b mn (∂ c r a nm ). Note that in the special case of a two-band Dirac dispersion, it can be shown that the general expressions for the shift current present here can be connected by a pull-back mapping to the Christoffel symbols of the geometric connection on the generalized Bloch sphere [49]. While this reinforces the association of the shift current with a semiclassical acceleration, for the multi-band case discussed here, the more convenient expression is the one presented in Eq. (3) in terms of λ ab (for another reparametrization cf. [50]).

B. Shift current of TBG
We model TBG using a modified form of the Bistrizer-Macdonald continuum model [36,38]. We attach two monolayer graphene sheets, with first rotating them with respect to one another with an angle θ = 1.05 o , and then introducing inter-layer coupling. By construction, this model is endowed with C 3z symmetry, and since both monolayers are inversion symmetric, the model as a whole has inversion symmetry. Inversion symmetry breaking is introduced by a staggered potential ∆, as it typically arises from encapsulation of the bilayer in hBN. On the other hand, C 3z is not necessarily broken because it requires some uniaxial strain of magnitude ε in the bottom layer. The parameters used here all correspond to the original values used in the Bistrizter-Macdonald construction. A detailed description of the model, as well as the results when corrugation effects are included, are given in the supplementary material. We emphasize that such a non-interacting band structure will of course not capture the important effects of many-body correlations in TBG [51]. However, since the BVPE across the gap relies on photon energies which greatly exceed the onset temperature of the ordered phases that form in TBG, the qualitative features of the shift current reported here can therefore be expected to also apply for a wide range of temperatures as long as the interacting system remains in the normal state. The resulting band structure in the mini-Brillouin zone of valley K is shown in Fig. 2a, both for the unstrained and strained TBG. The flat bands near half filling (µ = 0 meV) are gapped at Γ and K due to inversion symmetry breaking. Fig. 2b shows the Berry curvature in the mini-Brillouin zone in unstrained TBG, which is C 3 -symmetric. We note that the Berry curvature near the K and K -points in the original dispersion have opposite signs and cancel.
While the shift vector S in Eq. (2) itself is not expected to be further decomposable into gauge invariant pieces, the expression that enters the shift current can indeed be decomposed. Writing σ (s) = σ (s1) + σ (s2) , with Here, Eq. (7) contains the contributions from the Berry curvature dipole ∂ a Ω ac , which vanishes in the presence of C 3z symmetry. This symmetry renders only two components of σ ab;c independent. For simplicity we focus in the main text on σ yy;x , which encodes the transverse nonlinear conductivity response. The other independent component is σ xx;y (Supplementary Figs. 2 and 3); the symmetry analysis can be found in the SI. In Fig. 3a we present the total shift current in the unstrained case, for three values of the chemical potential which reside within the three gaps that are opened by the staggered potential. Within our numerical accuracy, the current is indeed found to have no contributions from the Berry curvature dipole, i.e. σ yy;x (s2) = σ yy;x (s) . Irrespective of this, for frequencies which cross the single-particle gap the shift current reaches giant values nearly 200 µAV −2 nm, far exceeding predicted values for other non-magnetic materials [52,53].
As we pointed out, from the quasiparticle shift it is not immediately obvious where such a giant nonlinear response could originate from. For this reason, we examine the different contributions in Eq. (8). Fig. 2c shows the momentum-space structure of the integrand in Eq. (8) at µ = 0, ω = 1.4meV, i.e. just above the band gap of the flat bands. The largest contributions are from the residual dispersion around Γ and from the flat parts of the band structure around K s , both of which are fairly broad in momentum space. Fig. 3b illustrates the relative contribution of the pieces K aa;c shift and R aa;c shift in Eq. (8), as a function of the distance from the gap at 0. Over a large range of frequencies, it holds that K aa;c shift R aa;c shift , meaning that the shift current is produced mostly by the phase of the Berry connection and not the shift of wavefunction centers. For comparison, we performed the same breakdown for the much more dispersive, gapped material MoS 2 , which has very similar symmetry properties [ Fig.  3b]. In this latter case, K aa;c shift R aa;c shift near the band edge, which supports the notion that the shift current is resulting from a displacement in both real space and momentum space, with the latter being greatly enhanced for a flatband dispersion. As we mentioned in the beginning, these statements might seem questionable upon regauging, a possible shortcoming on which we comment in the discussion.

C. Effect of strain
In the presence of uniaxial strain, the symmetry group of TBG is reduced to C 1 . All components of the conductivity are now independent, but for clarity we continue to examine only the transverse contribution σ yy;x . Since the Berry curvature dipole contribution as given by Eq (7) is no longer zero, we present both its contribution (Fig. 4a) and the remaining terms in Eq. 8 for σ aa;c (s2) (Fig. 4b). The total shift current is depicted in Fig. 4c, with the relative sizes of R aa;c shift and K aa;c shift shown in the inset. At moderate uniaxial strain of size = 0.1%, the conductivity σ (s1) due to the Berry curvature dipole becomes comparable in size to σ (s2) , but it has consistently the opposite sign. Thus, while the total conductivity unsuprisingly increases due to the reduced symmetry in the system, it has less accentuated resonances for the transitions at the band edges, with values up to 300 µAV −1 nm across the flatband gap. More unexpectedly, the shift current does not seem to profit from a net Berry curvature dipole, as this contribution either subtracts from the remaining current, or is almost negligible for transitions across the flatbands. Our results establish that a large anomalous acceleration in the quasiparticle motion can arise even if ∂ ka Ω ac = 0, i.e. its existence does not rely on the presence of a finite Berry curvature dipole in the system. This is an important distinction between the bulk photovoltaic effect and the anomalous Hall effect in TBG. We further observe that in all cases the shift current at frequencies ω 2 meV, corresponding to transitions between dispersive bands, is negligible small compared to the resonances around the band edges. In other words, the giant shift current shown in Fig. 3 is not tied to the topological properties of the flat bands but rather to their non-dispersive nature. We note that the maximal shift current obtained for TBG is larger than previously reported values for comparable, non-magnetic two-dimensional materials by a factor of 5 or more [53][54][55]. Figure 3c also shows the shift current evaluated at room temperature, T = 300K (dashed lines). While the transitions across the gaps of the dispersive bands are strongly suppressed, the large density of states in the flat bands supports a strong signal for ω = 4meV, with an amplitude of still 5% of its value at T = 0. We note that these conclusions contain the effect of disorder broadening through the inclusion of a finite quasiparticle relaxation rate [56].

D. Real and momentum space displacements
In introducing R shift and K shift , we are able to distinguish the sources of displacement that the quasiparticle suffers. These quantities are, however, not measurable observables of the system. In the following we explain  (3) into the R yy;x shift and K yy;x shift pieces, respectively. Solid lines depict TBG, while the dashed lines show MoS2, for comparison. The shift current is essentially zero for x < 1, has a large resonance around x ≈ 1 and then decays quickly for values above that. The normalized fraction does not show qualitative changes at x ≈ 1 but consistently conforms with the semiclassical intuition that the shift receives contributions from both the dispersive acceleration and the anomalous acceleration, where the latter is dominant for a flat-band dispersion. The maximal conductivity for µ = −5.0meV is attained at ω = 1.4meV, with value σ yy;x = +198µAV −2 nm, and for µ = 6.0meV it is −168µAV −2 nm. why such a decomposition is nonetheless insightful. For this, first recall that all types of topological bands (including flat bands) have no uniquely defined center-of-mass coordinate within the unit cell, because the topological nature of the bands prevents such an assignment. However, from this it does not follow that the momentum space integral of R shift can take arbitrary values, because it contains much more specific information about the relative positional difference between two bands, summed for all momenta. Indeed, it was already pointed out a long time ago [4,57] that for an arbitrary band structure one cannot generally expect to find a gauge such that r a mm − r a nn consistently vanishes for all momenta and all bands.
Drawing from these observations, we therefore suggest that a useful indicator for band flatness is that the integrated positional difference r a mm − r a nn between Bloch wavefunctions can be made substantially smaller than the integrated phase contribution ∂ ka arg r mn . We further conjecture that for highly dispersive bands a similar statement should hold about the smallness of the integrated phase contribution. A paradigmatic example in the latter case is a two-band semimetal with one band crossing. There, a mostly smooth gauge is at the same time periodic (i.e. without phase jump at the Brillouin zone boundaries), thus completely eliminating the phase contribution from the integrated shift vector. This is the expected result for a quasiparticle with vanishing effective mass which is changing position in an applied electric field. We remark that the difficulties in separating realspace and momentum-space effects of the acceleration into gauge invariant pieces are intrinsic to the more complicated semiclassical motion arising at second order in the applied field. In particular, the analogous splitting of the quasiparticle velocity into the regular (dispersive) and anomalous velocity has the important distinction that these two components of the velocity are orthogonal to each other, making them linearly independent. Such a decomposition is not straightforward for the acceleration, because it describes the changes to both regular and anomalous velocity components in both normal and perpendicular direction, thus mixing them. This being said, we believe that the conclusions outlined above can be made more rigorous by deriving a lower bound for R shift and K shift , which can then serve as useful indicator for mechanism of shift current generation in a given systemeither by a displacement in real space or one in momentum space. This will be the subject of a future work. Recent advances in ab-initio modelling of TBG [58,59], coupled with observations of the importance in optical properties of the off-diagonal components of the position operator in the Wannier basis [60], suggest that a fully ab-initio approach might alter the relative magnitude of R shift . We stress that this is not the case for TBG near the magic angle. As Ref. [60] found, for atomically localized Wannier orbitals, the contributions to the off-diagonal parts of the position operator decay quickly beyond nearest-neighbor (NN) coupling. This occurs at the scale a 0 , where a 0 is the monolayer lattice constant, while the continuum model's position operator scales with L m , L m being the Moiré unit cell length. Consequently, all such corrections, near the magic angle -where L m a 0 -are expected to be negligible and will not affect our results.

III. DISCUSSION
The shift vector S has previously been connected to the real-space shift of the center-of-mass coordinate between two eigenstates upon excitation from the conduction band m into the valence band n. To give some intuition, we expand the real-space representation of the periodic eigenfunctions |u nk in terms of local Wannier orbitals |w nR with center coordinate R [61] Then, in momentum space the Berry connection is given by Evaluating the shift vector S c mn = r c mm −r c nn +∂ kc arg r c mn based on this representation yields r a mm − r a nn = R a mm − R a nn for the direct difference. R a mm refers to the component of the center of the m-th Wannier function, in the a-direction. This is supplemented by the phase derivative ∂ ka arg r mn , whose integral over the Brillouin zone is a multiple of 2π. If only two Wannier orbitals have a significant overlap, the modulus |S| is clearly bounded by |R a mm − R a nn | < a, with lattice constant a. If several orbitals overlap, the phase factors in the sum Eq. (10) become important, with slope of growth in momentum space being at most a. Then, the phase derivative is expected to contribute similarly at O(a) to the shift vector. This already indicates that interpretation of the shift current as a result of the wavefunction shift is narrow to some extent. For ease of illustration, imagine a set of Landau levels in symmetric gauge. Their center-of-mass coordinate R can be moved around freely in exchange for acquiring an additional phase factor. Indeed, for generic flat bands it is to be expected that there is a gauge choice which makes the shift S only depend on the phase, because the center of the Wannier functions can be repositioned with an appropriate gauge transformation. This is inconsistent with the interpretation of the shift current as the real space shift of the wavefunction center upon absorption of a photon, as the wavefunction only suffers a phase shift.
If the shift current is instead viewed as the anomalous acceleration that a quasiparticle undergoes due to the interaction with the electric field at second order, the photogalvanic response follows as a straightforward generalization of the linear response formalism involving the anomalous velocity [31,62], thus removing the direct inference of a current from a real-space displacement. Instead, both R shift and K shift appear as the result of the same acceleration that changes both the position and the wavevector of the quasiparticle. As shown in the last section, this is consistent with our numerical findings for R shift and K shift using Bloch wavefunctions for a mostly smooth gauge choice within each band [cf. Fig. 3b].
While one might object against inspecting gaugedependent quantities, we emphasize that R shift and K shift can still contain valuable information about the quasiparticle dynamics in the sense that while they are not unique, this does not at all imply that they are arbitrary. Most importantly, based on our results we conjecture there exist nonzero lower bounds for both R shift and K shift which allow to uniquely identify the shift in terms of a real space or momentum space displacement. These bounds generalize previous constraints on the shift vector (viewed as the real space displacement of the quasiparticle wavepacket) [54] and extend to include the momentum-space shift introduced in this work. We also remark that the general principles outlined here, inferred from the interplay of K shift and R shift , as for the magnitude and resonant features of the shift current, are valid even when extrinsic effects are included, such as corrugation and strain (see Appendix). This is because the very nature of the effect relies on flat bands, which survive beyond the original Bistrizer-Macdonald parametrization.
Finally, we comment on the stability of the shift cur-rent signal under experimental conditions. While disorder broadening and thermal broadening present the dominant limitations for the quasiparticle lifetime, more elaborate disorder effects like skew scattering [63,64] might in principle affect the shift current. In TBG we do not expect skew scattering to play an important role because it is enhanced only in fairly clean systems which additionally feature an asymmetric dispersion near the band edges.
Both are conditions which are not met in TBG. Another important source of disorder are variations in the twist angle, which are known to be present in TBG devices [65]. Experimental results for the linear conductivity do not indicate qualitative changes as a function of the twist angle [66], indicating that angle disorder has only a limited influence on transport. Additionally, if time-reversal is broken, for example by some magnetically ordered state or by the presence of additional relaxation channels, ballistic currents may appear [67]. For all these reasons, and based on the line shape, an exciting application of the strong nonlinear signal could be to employ the shift current spectrum to determine the direct band gap in TBG in the normal state. We note that for circular polarized light, TBG exhibits a chiral photogalvanic effect [68]. Also, in a calculation including magnetic order at filling 3/4, Ref. [37] has reported a large photogalvanic current at for frequencies of 30 meV and above. This result, however, should be read in the context of time-reversal breaking. As previously shown [31,69,70], a system with broken time-reversal symmetry will generate a large injection current, that will scale inversely as γ −1 , where γ is the carrier relaxation rate. Conversely, the shift current discussed here is actually lifetime independent (γ 0 ), so both effects are experimentally distinct. The decomposition of the nonlinear current at 3/4 filling is discussed separately in App. F, where the shift current is much smaller than the injection current.
Regarding the possible application of our results for THz sensing, we note that in contrast to existing proposals for Terahertz devices, our proposal relies on the BPVE, and can thus circumvent several issues exhibited for example by Schottky diodes [71] and Bolometric devices [72]. While the former has been known to generate a broadband response in the THz-range, the conversion efficiency is low and of the order of a few percent, requiring amplification and complicated electrical circuitry. The internal p-n junction used in Schottky diodes requires sensitive doping, and is not as tunable as van der Waals systems (such as TBG) and operates under a bias field, which is highly sensitive to temperature. Bolometric devices add the additional complication of requiring thermal or mechanical junctions which typically have sub-optimal noise characteristics, thus reducing the applicability of such devices. We stress that setups which rely on the BPVE do not require external biases, amplification, thermal or mechanical junctions, but rather produce a current through resonant absorption of light.
In summary, we report a giant photogalvanic current for TBG which is irradiated by linear polarized light in the terahertz range. The magnitude of this shift current exceeds any previous reported numbers in comparable two-dimensional materials [53][54][55]. Our mechanism for the giant response, which is due to a momentum-space shift and acceleration of quasiparticles is a new facet of nonlinear current generation. This goes beyond ordinary mechanisms for photoconductivity enhancement, such as a large jDOS. The resonance profile we observe in both the strained and unstrained cases suggest that TBG is a promising candidate for THz detection and circuits, even at room temperature. Since the shift current is robust against thermally excited carriers, because it is a coherent bulk effect driven by the quantum geometry [26,49,54], we believe it may substantially improve detection capabilities for terahertz radiation. It might also increase the photovoltaic efficiency for energy harvesting, and have further applications in medical imaging, single-photon circuits and novel electronic devices [73,74]. As the root cause behind the large response we identified the anomalous acceleration due to the skew symmetric properties of the quantum geometry of the band structure as encoded in λ ab mn , which always appear in the shift current, but are greatly amplified in TBG due to the flat-band dispersion. The latter also turned out to be particularly important for retaining a large shift conductivity at higher temperatures. Our findings present a new design principle for shift current generation which is particularly suitable for twisted heterostructures. We expect the transverse dc-current reported here to be accessible using current samples and measurement techniques [75]. The line of reasoning developed in this work can potentially shed light on the quantum geometry of the band structure in similarly twisted van-der-Waals materials with nearly flat bands, for example MoTe 2 or WSe 2 [76].
The continuum model for TBG [36,38] is well known and only repeated here for convenience of the reader. It is constructed by joining together two monolayer graphene layers at zero effective separation between them. We choose the following real space unit vectors, for each graphene layer, In this description, the A and B sublattices are located, respectively, at v A = (0, 0), v B = d(0, 1). The reciprocal lattice vectors are, and the Brillouin zone corners hosting the low energy states are at K u = ± 4π 3d √ 3 2 , 1 2 . Within the BM model, the bilayer system is symmetric under C 3 by construction, and inversion and time-reversal (when both K u valleys of the original graphene monolayers are included). In order to break inversion symmetry, we introduce a coupling to a substrate (for example, hBN), which lifts inversion symmetry but leaves C 3 symmetry intact. This allows for a finite shift current which is entirely independent of the Berry curvature dipole, because the latter is set to zero by C 3 symmetry. The Hamiltonian of the bilayer system is therefore given by, where t,b,tb denote the top, bottom layer, and interlayer hopping respectively. The top layer has the following continuum Hamiltonian, for a given momentum q, where s, u designate the spin and valley degrees of freedom; i.e., s =↑, ↓, u = ±1, for the K, K valleys of the original graphene monolayers. a t/b,s,u is the annihilation operator for an electron with spin s, and valley u, on the A/B sub-lattices of the top/bottom layers. R ± is the rotation matrix for the top/bottom layers, given by: R ± = R ± θ 2 = cos( θ 2 ) ∓ iσ y sin( θ 2 ), acting on the sub-lattice space, with θ ≈ 1.05 • denoting the twist angle. The Fermi velocity is v f = 0.596eV nm [77], and σ = (σ x , σ y , σ z ) are Pauli matrices. The momentum q is measured relative to the valley centered at K u . Inversion symmetry breaking and uniaxial strain are introduced in the bottom layer. Throughout this work, the strain is applied along the zigzag direction of the bottom graphene sheet. The Hamiltonian of the bottom layer then takes the form, Here, is the uniaxial strain matrix, which has the form = −1 0 0 ν . A is the pseudo-gauge field resultant from the application of strain [78,79]. ∆ = 17 meV, is the staggered potential generated by alignment with an hBN layer, which is applied to the bottom layer only, mimicking realistic symmetry breaking in experiments. The staggered potential term is chosen in such a way as to break inversion symmetry (P), and C 2x [80]. The evolution of the band structure with the staggered potential for θ = 1.05 o is presented in Fig. 5. Strain is introduced through the parameter = 0, 0.1%, representing the strain-less and strained cases, respectively. With strain H tb is changed accordingly, as discussed in Sec. B

Appendix B: Continuum model under strain
The Bistritzer-MacDonald continuum model contains an inter-layer coupling term, which couples two mo- 2 , 1 2 , are the Moiré lattice vectors with q 0 = 8π sin( θ 2 ) 3 √ 3d , and d = 1.42Å is the carbon-carbon bond length in graphene [36]. In this framework, the interlayer coupling is included via the H tb term in Eq. A3 as, T 2,u (q, q ) + T 3,u (q, q ))a b,s,u (q ).
In the presence of the strain defined in the main text, the coupling matrices (acting on the valley index) become, Throughout, we take t = 0.33 eV. Recent work on relaxation of twisted graphene bilayers [81] suggests that interlayer hopping, t is modified due to resultant corrugation of the graphene layers. The effect of this on the response is minor, as shown in the SM Sec. E. Accordingly, the lattice vectors of the Moiré superlattice are deformed in the presence of strain. These have the form, The Dirac points transform under strain as, We introduce a pseudo-gauge field which stems from the underlying two-center approximation for the tunneling matrix [82]. It is given by, with ν = 0.165, β = 1.57 obtained from monolayer graphene. Finally, the diagonalization of the Hamiltonian H in Eq. A3 is accomplished by recasting it in the form, where now A s,u (q) = [a b,s,u (q), a t,s,u (q + q 1,u ), a t,s,u (q + q 2,u ), a t,s,u (q + q 3,u )] T is the infinite-component operator vector satisfying the constraints on momentum transfer. h u has the following truncated structure, after applying the momentum transfer relations ensuring non-zero T tunnelling, with H (i) t,u = H t,u (q + q i,u ). In this work, we used 81 sites in reciprocal space for the construction of the hamiltonian, which results in a Hamiltonian which is 324 × 324. The integrals appearing in Eqs. (7), (8) are computed using a discretized grid of 600 × 600 in the (k x , k y ) plane of the mini Brillouin zone. 10 3 frequency point samplings in ω are carried out uniformly. Convergence is checked against the case with C 3 symmetry, where verification is done by comparing σ xx;x with −σ yy;x ; and σ yy;y with −σ xx;y . All equalities were verified to within 5%. We stress that both original graphene valleys K ± are included in every calculation (see SI for model construction details), and, since the expressions appearing in Eqs. 3-6 are timereversal symmetric, the effect of evaluating both valleys is to double the overall result. This was also verified numerically to the stated accuracy. The delta functions of Eqs. (7), (8) are broadened with a width Γ = 0.02meV for T = 0K and Γ = 0.1meV at T = 300K. This corresponds to transport lifetimes observed in bilayer graphene in the clean limit [56].
Appendix C: Symmetry properties of the response In the main text, we observe that the response tensor σ ab;c has only two independent components, and that the Berry curvature dipole, ∂ a Ω bc vanishes. We proceed to prove this. The generator of C 3z is given by [83], The Berry curvature dipole (BCD), D abc is a gauge invariant material property of the system, and has the form D abc = ∂ ka Ω bc , making it a rank-3 pseudotensor. Firstly, we observe that this quantity is antisymmetric in (b, c). Applying Neumann's principle, we enforce D abc = αβγ M aα M bβ M cγ D αβγ . Using the antisymmetry of D abc = −D acb , we focus only on nontrivial components, Since a = x, y, we obtain the following set of equations, This set admits only the solution D xxy = D yxy = 0, as required, demonstrating that the BCD is zero, under C 3z . For the general rank-3 symmetric conductivity tensor σ ab;c , we derive analogous symmetry constraints under C 3z . This results in two independent components overall, σ xx;x = −σ yy;x = −σ xy;y = −σ yx;y (C2) σ yy;y = −σ xx;y = −σ yx;x = −σ xy;x . (C3) We further note that under linear-polarized light, the conductivity tensor exhibits a special permutation symmetry, σ ab;c = σ ba;c .

Appendix D: Additional data for the transverse components
For completeness, we provide the remaining transverse component σ xx;y of the shift current, which agrees qualitatively and in part quantitatively with the component σ yy;x discussed in the main text. We recall that in general, with C 3z symmetry the conductivity tensor has only 2 independent components, therefore the longitudinal components can be deduced straightforwardly from the data presented here.
With finite strain, all terms in Eq. (3) of the main text contribute to the current. For consistency, we again present the transverse component, σ xx;y , for 3 chemical potential values, in the same way as in Fig. 3 of the main text. With the introduction of strain, the Berry curvature dipole contribution is no longer zero (Fig. 8a). Although the net conductivity is not substantially enhanced by the introduction of strain (cf. Fig. 8c), a broad resonance in σ xx;y appears at chemical potential µ = 6.7meV. Note that while the introduction of strain induces a large Berry curvature dipole, it is still smaller than the remaining contributions according to Eq. (8), meaning that the sign of the total conductivity is determined by the latter part of the response.
To examine whether our results depend of the sign of the applied strain (i.e., whether the strain is compressive or tensile), we show in Fig. 9 the conductivity σ yy;x for strain with negative (i.e., compressive) magnitude, = −0.001. While the detailed frequency dependence is indeed sensitively dependent on the strain, both the magnitude of the response and ts resonance structure are very comparable. One may also be tempted to compare the magnitude of the shift current of ∼ 200µAnmV −2 with the nonlinear anomalous Hall signal, which is expected to be around 10x larger [43] for the same value of strain. However, we emphasize that the latter appears in response to a static electric field, rendering such a comparison moot; the mechanisms for the nonlinear Hall effect and for shift current generation are unrelated. We note that for µ = 6.7meV, the response profile is exceptionally broad, and is the conductivity is almost constant for the range ω = 2 − 6meV. Where u is the ratio of AA to AB domain tunneling. The case of u = 1 reproduces the original B-M parameters, while u = 0 would correspond to the so-called "chiral limit" hypothesized to occur under certain conditions in TBG. We stress that the latter has not been observed experimentally. In what follows, we consider the parameters suggested by Koshino et al. [81], and adopt u = 0.81, as the value extracted from ab-initio calculations. The substitution of u is repeated for all T 1,u , T 2,u , T 3,u matrices. The coupling to the hBN substrate is unchanged and remains at ∆ = 17meV. In Fig. 10 (b) we plot the band structure with the modified band structure with strain (ε = 0.1%), and without strain (ε = 0%). For completeness, we enclose alongside it, in Fig. 10(a) the band structure with the original B-M parameterization. While certain differences are discernible (such as a shift in the energy distance to the dispersive bands), flat bands appear as they did for the original B-M model. Following the differences observed in the band structure, we calculate the shift current with the different value of u. Clearly, the flat band contribution (seen at µ = 0) produces a large shift-current response in the region ω < 10meV, as seen indeed Fig. 11 (b); compare this with the original B-M parameters 11 (a), which similarly show this, albeit at a slightly different frequency ω. The dispersive bands enter at higher frequencies as expected -ω > 15meV -but their shape resembles the one found in Fig. 3 (and shown in Fig.  11(a)). This precisely underscores the point made in the main text that the parameters of the continuum model affect the results qualitatively, but the salient features of the shift current response and our suggested design principle remain unaffected. For completeness, we include the results with strain, in Fig. 11. Here once more we observe that the main effect of the altered parameters is manifest in the magnitude of the response, but not in the main conclusions we provided in the main text, namely, the robustness of the shift current response stemming from flat bands, which are preserved even in the presence of strain. We find an enhancement of the conclusions we have derived the main text: while in the presence of strain the dispersive bands contribute to the shift current response mainly via the presence of a non-vanishing Berry curvature dipole, the flat band response is once again dominated solely by the momentum-space shift current K shift we have put forward as the mechanism behind shift current generation. In summary, the tuning of interlayer tunneling does not substantially affect the quantitative results we have observed in the main text. Reducing u to   Figure 13. Photoconductivity for TBG at 3/4 filling. The red curve indicates the total photoconductivity (shift and injection). The blue curve denotes the shift current contribution. The yellow curve shows in the injection current contribution. The response is overwhelmingly determined by the injection current. The value taken for the intraband relaxation is γ = 0.1meV. The delta function is broadened by Γ = 0.02meV. u = 0.81 (as suggested by ab-initio studies on the properties of graphene bilayers) does not modify our conclusions regarding the robustness of the shift current stemming from flat bands, and the importance of the momentum space picture for understanding the source of the giant response in TBG. We have shown that this is independent of the precise parameterization of the continuum model for the TBG, and therefore, could serve as a vital design principle in a new generation of photovoltaic devices.