Transport signatures of Van Hove singularities in mesoscopic twisted bilayer graphene

Magic-angle twisted bilayer graphene exhibits quasi-flat low-energy bands with Van Hove singularities close to the Fermi level. These singularities play an important role in the exotic phenomena observed in this material, such as superconductivity and magnetism, by amplifying electronic correlation effects. In this work, we study the correspondence of four-terminal conductance and the Fermi surface topology as a function of the twist angle, pressure, and energy in mesoscopic, ballistic samples of small-angle twisted bilayer graphene. We establish a correspondence between features in the wide-junction conductance and the presence of van Hove singularities in the density of states. Moreover, we identify additional transport features, such as a large, pressure-tunable minimal conductance, conductance peaks coinciding with non-singular band crossings, and unusually large conductance oscillations as a function of the system size. Our results suggest that twisted bilayer graphene close the magic angle is a unique system featuring simultaneously large conductance due to the quasi-flat bands, strong quantum nonlinearity due to the Van Hove singularities and high sensitivity to external parameters, which could be utilized in high-frequency device applications and sensitive detectors.

Most of the experimental and theoretical research on TBLG has so far focused on observables in the thermodynamic limit and on transport in macroscopic samples in the semiclassical regime [25,26]. The quantum transport studies have so far addressed specific questions such as the angle-dependent minimal conductivity, disorder effects, and emergent magnetic textures in driven TBLG [27][28][29], as well as transport in crossed graphene nanoribbons, where the scattering region is smaller than the magic-angle moiré unit cell [30][31][32]. However, it remains an outstanding challenge to understand the effects of quasi-flat bands, VHSs and the rich variety of possible Fermi surface topologies on the quantum transport in mesoscopic TBLG samples, where the transport characteristics can be probed in an energy-resolved fash-ion. Recently, first steps towards this direction have been taken by identifying quantum transport signatures of VHSs in TBLG in the regime of intermediate twist angles 3 • θ 10 • assuming semimetallic leads [33]. In the semimetallic leads the DOS goes to zero at the quasi-flat band energies thereby hindering the identification of quantum transport signatures related to the VHSs within the quasi-flat bands close to the magic angle θ ∼ 1 • . Here, we go further into this direction by studying four-terminal conductance in mesoscopic, ballistic TBLG samples, containing approximately one million sites in the scattering region, around the first magic angle. Importantly, we obtain a higher energy resolution by using metallic leads and, therefore, we are able to study the effects of the VHSs and of the Fermi surface topology of the quasi-flat bands on the quantum transport.
We find that the low-energy quantum transport in TBLG close to the magic angle is affected by several factors. We demonstrate that by tuning the twist angle or pressure to flatten the energy bands the system can support considerably larger minimal conductance than monolayer and Bernal-stacked bilayer graphene devices [34][35][36][37][38]. We further link energy-dependent conductance signatures to different VHSs in the bulk DOS and to non-singular band crossings, and we observe unusually large conductance oscillations as a function of the system size. Our findings put forward that TBLG close the magic angle is an exceptional system combining large conductance originating from quasi-flat bands with strong quantum nonlinearity from VHSs and high sensitivity to external parameters. We propose that these properties could be utilized in compact high-frequency devices and sensitive detectors. Such applications could utilize the quantum twisting microscope technology offer-ing the possibility to locally and continuously tune pressure and twist angle in TBLG devices [39].
The structure of the paper is as follows. In Sec. II we introduce the four-terminal geometry for the transport calculations and the TBLG model used in this work. We continue by presenting in Sec. III the minimal conductance as a function of the twist angle, pressure, and the system size. In Sec. IV we identify the effects of VHSs and the Fermi surface topology on the energy dependence of the conductance. Finally, we discuss possible highfrequency device applications in Sec. V and we summarize our results in Sec. VI.

II. MODEL AND SETUP FOR TRANSPORT CALCULATIONS
A. Twisted bilayer graphene There are two common configurations for untwisted bilayer graphene (BLG) [40,41]: In AA-stacked BLG, corresponding atoms from different layers are on top of each other. In Bernal AB-stacked BLG, on the other hand, one of the layers is shifted relative to the other, such that some atoms from one layer lie at the centers of the hexagons formed by the other layer and vice versa. Both configurations have in common that their primitive lattice vectors are identical to those of single-layer graphene. In TBLG, starting from one of these configurations, the two layers are rotated relative to each other by an angle θ around a fixed point in space [42]. As a consequence, a moiré pattern emerges, which breaks the translational symmetry of the individual graphene layers. For certain angles, however, the two layers form a periodic moiré honeycomb superlattice, whose primitive lattice vectors and lattice constant are angle-dependent. These commensurate twist angles have the following form [43] cos θ = 3m 2 + 3mr + r 2 /2 where m and r are coprime positive integers. Based on this notation, the moiré lattice constant is [44] a = a 0 2 sin θ 2 r gcd(r, 3) , with the lattice constant a 0 = 2.46Å of single-layer graphene and gcd(p, q) denotes the greatest common divisors of the integers p and q. Generally, the smaller the twist angle, the larger the moiré lattice constant. At the first magic-angle θ = 1.05 • (m = 31 and r = 1), for instance, we have a ≈ 15 nm.

B. Four-terminal transport setup
In this study, we consider ballistic quantum transport through a TBLG region formed by two crossed graphene ribbons, as illustrated in Fig. 1. Such a setup could be realized in the lab with state-of-the-art experimental techniques [45] and the recent advances allow also to tune the twist angle and pressure locally and continuously [39]. We choose the ribbons to both have armchair terminations to avoid contributions from edge states to the electronic transport, as we aim to study bulk signatures in this work. This results in a parallelogram-shaped overlap region between the two ribbons. In our setup, the top ribbon can be twisted by an angle θ around the center of the bilayer region. The ribbons are placed such that for θ = 0 the two layers are stacked in an AA fashion in the overlap region. We note that we obtain qualitatively similar results as shown in this paper if we start from an ABstacked overlap region instead (not shown). Moreover, the continuation of the ribbons outside the bilayer region defines four semi-infinite, monolayer graphene leads, which we use for our transport calculations (see Fig. 1).
In contrast to previous work [33], we use metallic leads by tuning the chemical potential outside the scattering region far away from the Dirac-point energy (see below).

C. Tight-binding model
In the literature, various modelling approaches have been used to study the low-energy properties of TBLG, such as continuum models [2,46,47], ab-initio calculations [30][31][32], or tight-binding models [1,[48][49][50][51]. Here, we aim to investigate multi-terminal electronic transport in mesoscopic samples with a focus on the small-angle regime. Continuum models are long-wavelength lowenergy theories and, therefore, cannot fully capture the details at length and energy scales relevant to transport in mesoscopic systems. In the case of ab-initio calculations, it is computationally expensive to model samples sufficiently large to overcome finite-size effects at small twist angles, and thereby the effects attributable to the electronic properties of the bulk are necessarily obscured. Tight-binding models, on the other hand, are able to accurately capture the low-energy electronic properties of TBLG over a wide range of twist angles [49,50] while having the advantage of being computationally cheaper to scale up. For this reason, we adopt a tight-binding model approach in this work.
The generic form of a Hamiltonian for a weaklycoupled bilayer system, such as TBLG, is H = H 1 + H 2 + H 12 , where H 1,2 are the Hamiltonians of the individual layers and H 12 contains interlayer-coupling terms [33,49,50]. For the individual layers and the leads, we use the common nearest-neighbor tight-binding Hamiltonian of graphene [40] where c † imσ (c imσ ) creates (annihilates) a p z electron with spin σ =↑, ↓ at lattice site r i of the m-th layer, t = 3.09 eV is the nearest-neighbor hopping amplitude [49], and µ is the chemical potential. We set µ = −2 eV in the leads, such that the leads are metallic, and µ = 0 throughout the scattering region.
For the interlayer part of the Hamiltonian, we follow Ref. 49 by using where r ij = |r i − r j | is the in-plane distance between two lattice sites in different layers at positions r i and r j , respectively, and t (r) is the isotropic interlayer hopping integral given by with the nearest-neighbor interlayer coupling V 0 ppσ = 0.39 eV, the distance between the graphene layers d 0 = 3.35Å, and the decay parameter λ [49]. The form of the interlayer hopping term is based on a Slater-Koster approximation for the overlap integrals between the p z orbitals in different layers [52]. It is generally found that longer-range interlayer hopping terms have to be taken into account to capture the electronic bands of TBLG in a wide range of angles [49]. In accordance with Ref. 49, we use λ = 0.27Å, which reproduces the well-known band structures of AA-and AB-stacked BLG. Making use of the rapidly decaying nature of the hopping integral t(r) in Eq. (5), we further neglect interlayer terms with r > 5Å, which is sufficient to accurately capture the bands of TBLG at the first-magic angle.
We note that we neglect interaction effects in our model. These are generally important for the groundstate properties of TBLG when the Fermi level is within the quasi-flat bands and, thus, correlated phases emerge. However, this would lead to a reconstruction of the energy bands thereby obscuring the origin of the enhanced interactions, such as VHSs. We therefore restrict our study to a non-interacting description of the system, which is a good approximation as long as the Fermi level is tuned away from the quasi-flat bands. The transport can still be studied in an energy-resolved manner by tuning the voltage bias. In this case, only nonequilibrium quasiparticles are occupying the flat-band states so that interaction effects are not expected to be as important as in the case of equilibrium flat-band systems. Alternatively, it is also possible to screen the interactions so that a non-interacting description becomes more accurate [53,54]. Furthermore, the correlated phases appear only at low temperatures.
For the four-terminal transport calculations of the device shown in Fig. 1 we use the quantum transport Python package kwant [55]. To be able to capture bulk effects in the transport calculations, the size of the bilayer region needs to be at least on the order of several moiré unit cells. Using an interlayer-hopping cut-off as explained above enables us to efficiently study systems with close to 10 6 lattice sites. In particular, for Fig. 2 we perform calculations for samples of 3.2 × 16 magic-angle moiré unit cells, which is equivalent to 40 nm × 200 nm. For Fig. 4 we increased the size of the samples to 4 × 20 magic-angle moiré unit cells (50 nm × 250 nm). For such samples, the bottom-layer nanoribbon represents a long and narrow junction, whereas the top-layer nanoribbon realizes a short and wide junction (see Fig. 1). Our setup also allows us to study the effect of pressure which changes the value of the interlayer coupling.
We further aim to draw a connection between transport signatures and spectral features of the bulk, in particular VHSs. Therefore, for commensurate twist angles, we impose periodic boundary conditions on the moiré unit cell and calculate the band structure. This enables us to extract the bulk DOS and the Fermi surface of the system at fixed energies, which we use to pin down the VHSs. We note, however, that the observed transport signatures do not depend on whether the system is commensurate or incommensurate. Throughout this work, we align E = 0 with the energy of the Dirac points at the K and K points of the moiré Brillouin zone (BZ). For incommensurate twist angles, we interpolate linearly between the Dirac-point energies corresponding to the closest commensurate twist angles.

III. CONDUCTANCE AT THE DIRAC-POINT ENERGY
For a multi-terminal setup, the differential conductance G is a tensor defined through G ij = dI i /dV j , where I i is the current at the i-th lead and V j is the voltage at the j-th lead of the system. Here, we consider a system with four leads, so our conductance tensor is represented by a 4 × 4 matrix. Moreover, our setup has an approximate C 2 symmetry, which is only weakly broken along the boundaries of the system due to its terminations. Hence, the conductance tensor G ij has only four independent components, which we choose to be the intralayer conductances G 12 (short and wide junction) and G 34 (long and narrow junction), as well as the interlayer conductances G 13 and G 14 .
First, we compute the conductance G ij as a function of the twist angle θ at the Dirac-point energy E = 0, which we dub minimal conductance of the quasi-flat band system in analogy with monolayer graphene. In fact, for twist angles θ > 1.05 • the wide junction conductance typically has a local minimum at this energy. We note, however, that at smaller twist angles the conductance is no longer necessarily minimal at E = 0 for reasons explained in Sec. IV. We show our results in Fig. 2(a). At large angles, θ 4 • , the conductance shows a universal behavior independent of the twist angle. The shortjunction conductance G 12 approaches a value corresponding to the minimal Dirac-point conductivity of singlelayer graphene, namely G 12 W bottom /W top ≈ 4e 2 /πh, which is in agreement with previous results in the literature [27,34,35]. The long-junction conductance G 34 , on the other hand, approaches the quantized value of 2e 2 /h. We find that it corresponds to a single spindegenerate, propagating bulk mode confined to the bottom layer, whose presence is attributed to the particular width of the nanoribbon [56]. At the same time, the interlayer conductances G 13 and G 14 vanish. Hence, the two nanoribbons are effectively decoupled and the conductance tensor decomposes into two independent twoterminal conductances. This is in agreement with previous results [27].
For small angles, θ 4 • , the interlayer conductances become nonzero clearly indicating coupling between the two nanoribbons. This is also reflected in the behavior of the long-junction conductance G 34 , which is suppressed, because the propagating bulk mode in the bottom layer can now interfere with the corresponding mode in the top layer thereby lowering the conductance. Close to the magic angle, G 34 becomes nearly zero. On the contrary, the short-junction conductance G 12 is strongly enhanced as we get closer to the magic angle and deviates considerably from the single-layer value. The behavior of the different conductance channels can be attributed to the formation of quasi-flat energy bands in small-angle TBLG: the presence of the flat band enhances the DOS around E = 0, while the Fermi velocities of the corresponding bulk modes are decreased. A measurable enhancement of conductance requires a sufficiently large number of lead modes, which is why the wide-junction conductance G 12 shows the largest effect. Therefore, G 12 is a suitable quantity to probe how the van Hove singularities affect the transport properties through their enhancement of the DOS. We note that there are also other factors influencing the conductance beyond the bulk DOS, as we will discuss below.
The wide-junction conductance G 12 is generally enhanced towards smaller angles, but shows a nonmonotonic behavior. It features two pronounced maxima, one at θ = 1.14 • and one at the magic angle θ = 1.05 • , and increases again for θ 1.0 • . We discuss this behavior in more detail below in the context of the correspondence between conductance and DOS.
In Fig. 2(b), we show the minimal conductance also as a function of the nearest-neighbor interlayer coupling V 0 ppσ at the magic angle θ = 1.05 • . This parameter can be controlled by applying pressure [16,57]. Overall, all components of the conductance tensor G ij exhibit the same qualitative behavior as in the case of a variation of the twist angle. In particular, the conductances for vanishing interlayer coupling V 0 ppσ ≈ 0 are the same as in the large twist-angle regime confirming the picture of effectively decoupled nanoribbons. In the large-coupling regime, we further observe the same enhancement of the interlayer conductances and of the wide-junction conductance, while the long-junction conductance is suppressed. Close to the magic angle, the similarities between the variations of these two parameters, θ and V 0 ppσ , even show a good quantitative agreement, as can be seen by comparing the insets of Figs. 2(a) and 2(b). More generally, we find that also the corresponding bulk energy bands evolve in a similar way. Our results suggest that a TBLG sample at an incommensurate angle can be approximated by a sample at the closest commensurate angle in combination with applied pressure. Below, we will use this insight to continuously trace the evolution of VHSs as a function of the twist angle.
Finally, in Fig. 3 we compare the minimal conductances of our setup for small and large twist angles as a function of the system size with fixed aspect ratio.  large angles, the minimal conductance converges quickly to the universal single-layer graphene value, which is in agreement with an effective decoupling of the layers. On the contrary, as we get closer to the magic angle the minimal conductance is enhanced and develops pronounced oscillations as a function of the system size. Notably, the minimal conductance in this regime is much larger than the universal value of Dirac-point conductance for a single layer indicating that both layers now contribute to the electronic transport.

IV. CONDUCTANCE SIGNATURES OF VAN HOVE SINGULARITIES
Next, we compare the conductance G ij (E) of the TBLG device to the bulk spectrum of TBLG focusing on the energy regime around the quasi-flat bands. Fig. 4 shows our results for a few selected cases. To better link the features in the conductance to the DOS and the dispersion of the energy bands, we show the evolution of the Fermi surface at energies around these features in Fig. 5. More results can be found in the Supplemental Material (SM) [58]. We note that with decreasing twist angle, additional energy band crossings lead to a plethora of features in tiny regions of the moiré BZ. Even though we see some of these features as peaks in our finite-momentum grid calculations of the DOS, we do not observe conductance signatures that can be attributed to them for the system sizes considered. In the following, we will therefore restrict the discussion regarding van-Hove singularities to the most pronounced peaks in the computed DOS. We have checked that all of the discussed peaks correspond to ordinary van Hove singularities with a logarithmic divergence.
A. Twist angles θ ≥ 1.12 • We begin by discussing the conductance in the smallangle regime leading up to the first pronounced maximum in the minimal conductance G 12 (E = 0) at θ = 1.14 • [see Fig. 2(a)]. Fig. 4(a) shows the conductance, the DOS, and the energy bands at the closest commensurate angle θ = 1.12 • . As we decrease the twist angle to small values, the energy bands of the system are reconstructed due to the enlargement of the moiré unit cell. This reconstruction flattens the Dirac cones at the K and K points of the moiré BZ and eventually decouples the four corresponding spin-degenerate bands from the rest of the bulk bands thereby forming isolated quasi-flat bands. Close to the Dirac point energy E = 0, we observe two pronounced peaks in the wide-junction conductance G 12 that align with the two main peaks of the bulk DOS. These peaks originate from VHSs in the bulk energy bands. As illustrated in Fig. 5 [see panels (α)] for one of the singularities, each of them corresponds to a double saddle point involving two energy bands crossing along the ΓK lines of the BZ. At these points, the topology of the Fermi surface changes constituting a Lifshitz transition leading to a measurable transport signature. Similar signatures arising from Lifshitz transitions associated with van Hove singularities have been experimentally probed in the context of untwisted bilayer graphene [59][60][61].
In agreement with the literature, the DOS goes to zero at the Dirac point energy, which is similar to monolayer graphene. In monolayer graphene the conductance still takes a finite value, as discussed above. Here, however, the conductance G 12 at E = 0 is large and exceeds the monolayer and the Bernal-stacked bilayer value [38] considerably, as we already pointed out above. This is mainly due to the strong coupling of the layers in this regime, so the Dirac points of both layers now contribute to the conductance. Another reason for the enhancement is the broadening of the DOS peaks in the conductance of the device: as the bands become flatter, the two broadened VHS features move closer together and their overlap at E = 0 is enhanced. These are factors contributing to the increase of G 12 (E = 0) towards smaller angles seen in Fig. 2(a).
We further observe two smaller peaks around the two main peaks in the wide-junction conductance G 12 . In contrast to the main peaks, these do not correspond to distinct features in the DOS and, therefore, do not correspond to VHSs. Nevertheless, we find that they instead correspond to a band crossing at the M points in the BZ (see Appendix A). These band crossings constitute Lifshitz transitions that are not associated with VHSs.
B. Twist angles 1.12 • > θ ≥ 1.05 • As we lower the twist angle further from θ = 1.12 • down to the first magic angle at θ = 1.05 • , the bandwidth of the quasi-flat bands is reduced by nearly one order of magnitude [see Fig. 4(b)]. This enhances the DOS overall and, consequently, leads to a generally larger wide-junction conductance G 12 within the energy range of the isolated quasi-flat bands. We still observe two pronounced main peaks close to E = 0, which correspond to the same type of double-saddle point VHSs as before [see panels (β) in Fig. 5]. The band crossings responsible for the additional conductance features at θ = 1.12 • have moved to the band edges of the isolated quasi-flat bands and no longer stand out as much as before. We merely see small kinks in the G 12 conductance close in energy to these features. At the Dirac point energy, the DOS is no longer zero because of other parts of the moiré bands crossing this energy. As before, the wide-junction conductance G 12 is largely enhanced at E = 0 and almost one order of magnitude larger than the monolayer and Bernal-stacked bilayer value. This is due to several factors: the strong coupling between the layers, the broadening of the VHS features, and also additional states crossing the Dirac-point energy. Besides the two main peaks in the DOS, we note that there are also two smaller features in the DOS slightly below the Dirac point energy. We find that these correspond to ordinary VHS originating from saddle points close to the center of the BZ (see SM [58]). They are not clearly visible in the wide-junction conductance G 12 , because their contributions merge with the broadened features associated with the other VHSs. Note that these are the same singularities that are visible in the DOS of Fig. 4(a) in between and close to the two main peaks discussed in Sec. IV A. Fermi surface topology. However, we find that only the VHS indicated by (δ) in Fig. 5 close to the Dirac-point energy leads to a clear feature in the wide-junction conductance G 12 . The two VHSs indicated by ( ) and (ζ) above the Dirac points are too close in energy to be resolved separately and, together, they only cause a small kink in G 12 . We note that the band crossing at M , leading to a small conductance feature at larger angles [see Fig. 4(a)], is close in energy to these VHSs and, therefore, might also contribute to the kink in G 12 . Similarly, the related band crossing at M below the Dirac point is close in energy to the VHS at (γ) and, therefore, we attribute them together to a small bump in the G 12 conductance within this energy range.
In addition, we find another pronounced peak in the G 12 conductance roughly at the Dirac-point energy. This feature is relatively far in energy from all VHSs and, therefore, does not correspond to a singularity. Nevertheless, we note that it is close in energy to a tangential band touching point along the ΓM high-symmetry lines of the moiré BZ (see Appendix A).
Generally, we observe discrepancies between the energetic positions of features in the conductance and the associated spectral features of up to 0.1 meV for the TBLG devices considered. For angles θ ≥ 1.05 • , such discrepancies are small with respect to the width of the quasi-flat bands. For smaller angles, the features of the quasi-flat bands move even closer together leading to a larger relative discrepancy in this regime. In Appendix B, we ana-lyze how the system size affects the broadening and the energetic positions of the conductance features.

D. Evolution of van Hove singularities between commensurate angles
Instead of lowering the twist angle, we can alternatively increase the nearest-neighbor interlayer coupling V 0 ppσ to approximate the bulk spectrum at the next commensurate twist angle. We demonstrate this in Fig. 4(d) showing that the bulk energy bands and the bulk DOS are almost identical to those in Fig. 4(c). Also the conductances are in good qualitative agreement despite minor quantitative deviations. This makes it possible to unravel the evolution of VHSs between two commensurate twist angles.
We have used the interlayer hopping to continuously tune between the energy bands and DOS spectra corresponding to the commensurate angles θ = 1.05 • and 1.02 • (see SM [58]). Starting from the magic angle θ = 1.05 • [see Fig. 4(b)], the two VHSs below the Dirac point energy close to Γ merge, which is accompanied by changes in the Fermi surface topology close to the Γ point. We observe that a single VHS emerges from the two original VHSs, which moves further up in energy. It passes through the Dirac-point energy and eventually merges with the large double-saddle point VHS feature. The energy bands associated with the two singu-larities reconnect and form new singularities, which subsequently split in energy. These are the two pronounced VHS features we observe at θ = 1.02 • for positive energies. On the other hand, the double-saddle point VHS below the Dirac-point energy we observe at the magic angle θ = 1.05 • splits into two separate VHSs in the DOS.
A similar analysis can be performed to study the evolution of energy bands and VHSs between any two commensurate twist angles provided that the angle difference is sufficiently small. In particular, this applies to the small angle regime below the magic angle where the distance between adjacent commensurate angles rapidly goes to zero.

V. NONLINEAR TRANSPORT AND DYNAMICS
The predicted nonlinear current response to an applied voltage, following from the energy dependence of the conductance, suggests that TBLG systems can be utilized in various kinds of generators, frequency multipliers, frequency mixers, parametric amplifiers and detectors of electromagnetic radiation [62][63][64][65][66][67][68]. Here, we briefly comment on some aspects regarding potential applications of TBLG devices in the light of our transport results.
The simplest approach to describe the photon-assisted current I(t) as a response to a time-dependent voltage is based on the formula where J n (x) are Bessel functions, the summations are from −∞ to ∞ and β k = eV ω k / ω k . Here, I S (eV dc ) is the static current-voltage characteristic and K S is related to it by the Kramers-Kronig relation where P denotes the Cauchy principal value. The formula in Eq. (7) is a generalization of the Tien-Gordon-Tucker relations [62][63][64] to arbitrary polychromatic fields, and it has been utilized in the description of photonassisted transport in various quantum systems operating in sequential tunneling and miniband transport regimes as well as in Josephson junctions and exciton condensates [63][64][65][69][70][71]. From Eq. (7) it is easy to notice that for the type of applications discussed above it would be preferable to have large conductance and strongly nonlinear I S (eV dc ) characteristics. Namely, the magnitudes of the Fourier components of I(t) depend, in addition to the amplitudes of the driving fields, on the overall conductance. Additionally, it follows from the properties of the Bessel functions that for linear I S (eV dc ) characteristics the current response would be I(t) = G 0 V (t), where G 0 is the conductance. Thus, linear I S (eV dc ) characteristics are unsuitable for the use in generators, frequency multipliers, frequency mixers, parametric amplifiers or the type of detectors discussed in Refs. [63][64][65][66][67][68]. We also point out that for some applications it is advantageous if the energy scale of nonlinearities in the I S (eV dc ) characteristics roughly matches the photon energies ω k . According to our calculations, TBLG devices close to the magic angle typically have almost an order of magnitude larger conductance than monolayer graphene systems of similar size. Moreover, our calculations demonstrate that there exists a strong quantum nonlinearity due to the van Hove singularities. Furthermore, we point out that the energy scales of the nonlinearities vary within the range 0.1 − 1 meV, so that TBLG devices are promising also for high-frequency applications in the range of frequencies from 100 GHz to 1 THz. This frequency range is particularly important from the technological perspective because of the lack of compact solid-state technologies for generating and detecting THz radiation.

VI. CONCLUSION
We have studied low-energy, four-terminal conductance of twisted bilayer graphene in mesoscopic samples formed by crossing a narrow and a wide graphene nanoribbon. At large twist angles, the interlayer conductance vanishes because the two junctions effectively decouple. At small twist angles, however, the interlayer conductance becomes nonzero indicating coupling between the graphene sheets, whereas the narrow-junction conductance goes to zero. Notably, the wide-junction conductance is largely enhanced due to the formation of quasi-flat energy bands. We have shown that the conductance can also be controlled with pressure.
We find that the low-energy quantum transport close to the magic angle is affected by several factors. We have identified VHSs within the quasi-flat bands contributing to conductance features through their effect on the bulk DOS. We observe that such features are generally pronounced around the center of the quasi-flat bands, whereas they are suppressed closer to their edges. Moreover, we find additional peaks in the conductance not corresponding to features in the bulk DOS. These peaks coincide in energy with non-singular band crossings located at high-symmetry points and lines of the bulk BZ. Finally, we also observe a broadening of the various con-ductance features. As a consequence of the broadening and of the strong coupling between the layers, the quasiflat bands exhibit a large minimal conductance at the Dirac-point energy exceeding the value of single-layer and Bernal-stacked bilayer graphene considerably. We have further studied the continuous evolution of VHSs corresponding to different commensurate twist angles by tuning the interlayer coupling, which could experimentally be realized by applying pressure.
VHSs play a prominent role in the exotic phenomena of magic-angle twisted bilayer graphene by amplifying electronic correlation effects. Our results indicate that the VHSs strongly influence also the transport properties of this system, but we find that there are additional factors affecting them. Moreover, our calculations provide an estimate on the sample sizes required to resolve signatures of the flat-band VHSs in the conductance. Specifically, we have considered bilayer samples of size 50 nm×250 nm. Even though our study neglects interaction effects, our results are also directly relevant experimentally: in a transport-measurement setup the appearance of correlated phases could be avoided, for instance, by tuning the Fermi level away from the quasi-flat energy bands and by using voltages to probe energies within the flat bands, or by operating at temperatures above the critical temperatures of the correlated phases. In the case of extremely accurate low-temperature transport experiments the interactions could also be screened [53,54]. We have neglected the effects of disorder as they are expected to be unimportant for the system sizes considered. We have further disregarded lattice relaxation effects, which may play a role for twist angles 1 • due to the possible formation of alternating AB and AA stacking domains instead of the moiré pattern [72].
As an interesting future research direction, we propose the investigation of TBLG devices for the use as compact solid-state frequency multipliers, frequency mixers, parametric amplifiers and detectors operating at THz frequencies. The high sensitivity of the conductance to external parameters suggests that TBLG devices could be utilized also in various types of sensitive detectors. Figure S1(a) shows a band-crossing transition at the M point in the BZ. The corresponding energy of the crossing point is highlighted by the dashed black line below E = 0 in Fig. 4(a) of the main text. With increasing energy, the Fermi surface extends and intersects with itself at the M point. A similar band-crossing transition takes place above E = 0 where the Fermi surface intersects with itself at the M point (not shown). We find the same band-crossing transitions at the M point also in Figs. 4(b-d).
In Figs. 4(c) and (d), we further identified a feature in the wide junction conductance close to E = 0, which is far in energy from any VHS. Instead, it is close in energy to a band crossing point, whose energetic position is highlighted by the dashed black line in Fig. 4(c) close to E = 0. The corresponding evolution of the Fermi surface is illustrated in Fig. S1(b). With increasing energy, the circular feature along ΓM close to the center of the BZ grows, touches another Fermi line tangentially, and then crosses through it.
Appendix B: Effects of system size on the conductance In this section, we analyze how a change of the device dimensions affects the conductance. In particular, we independently vary the width W bottom of the long bottom junction and the width W top of the wide top junction of our TBLG device (see Fig. 1). We focus on the wide junction conductance G 12 and discuss, as an example, a device with a twist angle of θ = 1.02 • and interlayer hopping V 0 ppσ = 390 meV [compare to Fig. 4(c)].
We show our results in Fig. S2. We find that decreasing the wide junction width W top suppresses the conductance overall, which is due to the reduced DOS in the top leads. All the conductance features are still visible. Nevertheless, there is no effect on the energetic position and width of the conductance features. On the contrary, a change of the long-junction width W bottom affects the broadening and also the energetic position of the peaks. As W bottom gets smaller, the features become flatter and more smeared out. In particular, we observe that the kink-like features close to the edges of the flat band disappear for a slightly narrower system.