Transparent mirror effect in twist-angle-disordered bilayer graphene

When light is incident on a medium with spatially disordered index of refraction, interference effects lead to near-perfect reflection when the number of dielectric interfaces is large, so that the medium becomes a"transparent mirror."We investigate the analog of this effect for electrons in twisted bilayer graphene (TBG), for which local fluctuations of the twist angle give rise to a spatially random Fermi velocity. In a description that includes only spatial variation of Fermi velocity, we derive the incident-angle-dependent localization length for the case of quasi-one-dimensional disorder by mapping this problem onto one dimensional Anderson localization. The localization length diverges at normal incidence as a consequence of Klein tunneling, leading to a power-law decay of the transmission when averaged over incidence angle. In a minimal model of TBG, the modulation of twist angle also shifts the location of the Dirac cones in momentum space in a way that can be described by a random gauge field, and thus Klein tunneling is inexact. However, when the Dirac electron's incident momentum is large compared to these shifts, the primary effect of twist disorder is only to shift the incident angle associated with perfect transmission away from zero. These results suggest a mechanism for disorder-induced collimation, valley filtration, and energy filtration of Dirac electron beams, so that TBG offers a promising new platform for Dirac fermion optics.

Transmission of light through a medium with random refractive index is a well studied problem in optics [1][2][3]. Of particular interest is the problem of "transparent mirrors", in which the comprising elements are all transparent, but when arranged in a many-layered stack they form a medium with near-perfect reflection. Even though the first theoretical treatment of a similar problem was proposed two centuries ago by Fresnel [4], it was only relatively recently that Berry and Klein showed that the problem could be understood by analogy with Anderson localization [5]. The corresponding disorder-averaged transmission intensity decays exponentially with the number of dielectric interfaces.
The linearity of the energy dispersion for Dirac electrons offers an analogy with optics. This analogy has inspired previous authors to explore numerous analog optical phenomena in graphene, including but not limited to Veselago lensing [6], sub-wavelength diffraction [7], and the Goos-Hänchen effect [8]. For the same reason, the Dirac equation with velocity disorder strikingly resembles the problem of transmission of light through a medium with random refractive index [9][10][11]. The analogy with optics seemingly suggests that for a series of parallel one-dimensional domains of random Fermi velocity, one should expect exponential decay of the transmission coefficient. However, unlike in optics, Dirac electrons experience Klein tunneling, which implies that there is no reflection at normal incidence [12]. This key difference suggests the idea of an angle-dependent localization length, which diverges as the incidence angle goes to zero. For an incident electron beam with a wide range of incident angles, the angle-averaged transmission may have a decay that is slower than exponential.
Recent experiments mapping the spatial variation of twist angle in twisted bilayer graphene (TBG) provide a conspicuous example of a Dirac system with disordered velocity [13][14][15][16][17]. Since the Fermi velocity in TBG depends sensitively on the twist angle [18][19][20], the formation of domains with different twist angle indicates a Fermi velocity disorder. While the effect of scalar and vector disorder potential on the trans-port of Dirac electrons has been explored theoretically [21][22][23][24][25][26], the consequences of a spatially random Fermi velocity has not been addressed thoroughly. In this paper we study a model of electron transmission through a Dirac material with spatially modulated Dirac velocity. In particular, we work with the setup shown in Fig. 1, where the velocity is modulated across parallel, quasi-one-dimensional domains with uniform thickness D. We assume step-like changes in velocity, which corresponds to the limit where D is long compared to the electron wavelength and the thickness of the domain wall is small compared to the wavelength. Within this setup we consider two scenarios: It is worth mentioning that there are other methods to achieve Fermi velocity modulation in a Dirac material, in arXiv:2008.05481v1 [cond-mat.mes-hall] 12 Aug 2020 addition to variations in twist angle. Namely, one can use (i) strain engineering [27,28] (ii) dielectric screening [10,29,30] or (iii) a synthetic Dirac material with a super-lattice potential [31][32][33][34]. Manipulating Fermi velocity via mechanical strain is a particularly well-studied method. Since strain is a tensor, in general it can change the velocity in anisotropic ways. Strain also shifts the location of the Dirac cone in momentum space so that Klein tunneling is not necessarily preserved, although it is preserved in the case of uniaxial strain in the zigzag direction [35]. Dielectric screening by a substrate permits the modulation of Fermi velocity near the Dirac point via renormalization of electron-electron interactions. One can also make a synthetic Dirac dispersion by arranging atoms or creating an imposed potential energy with hexagonal symmetry.
Our main results can be summarized as follows. In Sec. I, we show that the transmission of two-dimensional massless Dirac fermions through a medium with quasi onedimensional random Fermi velocity can be exactly mapped onto the problem of conventional Anderson localization in one dimension. For small incident angle φ 0 , the localization length ξ varies as ξ ∼ 1/φ 2 0 [see Fig. 2(a)]. The diverging localization length at normal incidence is the manifestation of Klein tunneling in this system. The angle-averaged transmission in this system decays as N −1/2 , where N is the number of interfaces. In Sec. II, twist angle disorder in TBG is mapped onto the problem of a Dirac Hamiltonian subjected to both velocity disorder and step changes of a magnetic vector potential. When the momentum of the incident electron is large compared to the shift in reciprocal space produced by the vector potential, Klein tunneling is recovered almost completely, although the incidence angle associated with perfect transmission is shifted away from zero [ Fig. 2(b)]. The resulting transmission similarly decays as N −1/2 up to some very large N. In Sec. III, we relate our results for the electron transmission to the electrical conductivity and Fano factor. While we discuss our results everywhere in the context of graphene, it is worth mentioning that the results of Sec. I also hold true for any electronic system that exhibit properties which can be well described by the Dirac equation [36].

I. MASSLESS DIRAC ELECTRONS WITH RANDOM FERMI VELOCITY
In this section, we derive the transmission coefficient for Dirac electrons incident on the quasi-one-dimensional pattern of random velocity depicted in Fig. 1. We assume, for the moment, that only the velocity is modulated as a function of the coordinate, so that electrons are described by a massless Dirac equation with position-dependent velocity. In Sec. II, we consider a more realistic model for twist-angledisordered bilayer graphene, which includes the effects of the random gauge field that describes the shift of the Dirac point in momentum space.
The massless Dirac equation with position dependent ve- locity v F (r) is given by [10] where is the reduced Planck constant, σ is the vector of Pauli matrices, and E is the electron energy. The wave function Ψ (r) is a two-component spinor If we define the auxillary spinor Ψ (r) ≡ √ v F (r)Ψ (r), the Dirac equation becomes We consider the case where the Fermi velocity varies only in the x-direction, i.e. we assume v F (r) = v F (x). The translational invariance along the y-axis allows us to write the following ansatz Substituting this ansatz into Eq. (3), we get the following coupled differential equations We are interested in a situation where the Fermi velocity is piece-wise constant; let it be v F locally. In this situation, we can manipulate the above coupled first order differential equation to get the following uncoupled second order differential equation The solution to this equation is a linear superposition of plane waves moving in the left and right directions: where q ≡ E v F 2 − q 2 y is the wavevector along the xdirection. The angle φ between the x and y momenta of the electron satisfies Plugging ψ A into Eq. (6) and using the definition of φ we can write ψ B as Now we can write the whole Dirac auxiliary spinor Ψ as We have added an overall normalization factor of N. In this paper the normalization is chosen such that j x = v F Ψ † σ X Ψ, which is the probability current of Dirac electrons along xdirection, is normalized to |a| 2 − |b| 2 . This choice is N = 1 / √ 2 cos φ.
Before considering the transmission through an arbitrary number of interfaces, we first examine the case of a single interface [see Fig. 3] at an arbitrary position x = l. The Fermi velocity across the interface is given by where Θ (x) is the Heaviside step function. Using the conservation of y-momentum across the interface we obtain the analog of Snell's law: If v 1 > v 0 , then there exists a critical angle φ c = arcsin(v 0 /v 1 ) such that φ > φ c corresponds to total internal reflection. At such large angles the electron wave function crosses the barrier only through an evanescent mode, which decays exponentially in amplitude with x − l. In this paper we work in the limit where the domain width D is sufficiently large that these evanescent modes can be neglected (In the numerical results presented below and in Fig. 2 The boundary condition for the Dirac Hamiltonian at the interface is given by the continuity of both components of the auxiliary spinor Ψ. We therefore arrive at the transfer matrix M which relates the amplitudes a 0 , b 0 on the left region to the amplitudes a 1 , b 1 of the right region: where The complex numbers α and β are given by We can calculate the transmission probability across a single interface by defining the reflection and transmission amplitudes, r and t, respectively [see Fig. 3], such that The corresponding transmission probability T 1 = 1 − |r| 2 is With the help of Snell's law the transmission probability can be rewritten as where which is the consequence of Klein tunneling. This expression also gives T 1 = 1 when u 1 = 1, which is the condition of perfect transmission when there is no barrier at all. When there are multiple interfaces, the transfer matrices for each interface follow the composition rule, such that for N interfaces We can use this idea to calculate the transmission coefficient when there are two interfaces present. In this case we have three consecutive regions with velocities v 0 , v 1 , and v 2 respectively. The transmission amplitude t 12 across the two interfaces combined can be found to be which is equivalent to the following geometric series Here, r 1 is the reflection amplitude for the waves incident from the right on the first interface if it was the only one present. Equation (22) provides a simple geometric interpretation of Eq. (21): when we combine two interfaces the net transmission amplitude is the coherent sum of all the multiply reflected waves in the region between the two interfaces [5]. The net transmission coefficient, T 2 , across both interfaces is where T 2 and R 2 represent the transmission and reflection probability, respectively if only the second interface was present. Here, Φ 12 = 2q 1 D is the total phase accumulated during one complete internal reflection in the sandwiched region between the two interfaces. Let us calculate the average of the logarithm of the transmission coefficient across two interfaces. (Throughout this paper, we use "log" to denote the natural logarithm.) As we explain below, the transmission T N is not self-averaging in the limit of large N, and instead we deal with log (T N ), which is self-averaging. The average is taken over the values of the random velocities, u 1 = v 1 /v 0 and u 2 = v 2 /v 0 . We assume the following box distribution for the variables u i : Taking the logarithm of both sides of Eq. (23) and averaging over the disorder yields log (T 2 ) = log (T 1 ) + log (T 2 ) where . denotes the disorder average. The phase, Φ 12 , depends very sensitively on the value of u 1 when D is large compared to the electron wavelength, which is the case we are considering. On the other hand, R 1 and R 2 vary slowly with the velocity. Consequently, the term e˙ι Φ 12 oscillates rapidly with u 1 and averages to zero as we integrate over different values of the velocity, giving log (T 2 ) log T 1 + log T 2 . Equation (25) directly implies that the disorderaveraged logarithm of the transmission is additive when multiple interfaces are present. Generalization of this equation to many interfaces is now straightforward, since in the logaverage the interfaces are effectively decoupled: where In the second line of Eq. (26), there is an approximate sign because for the very first term log T 1 , averaging is done only over the value of u 1 . However, for every other interface i > 1, the averaging of log T i is performed over both u i−1 and u i . Equation (26) demonstrates that the average logarithm of the transmission is a self-averaging quantity. The exponential of the log-averaged transmission gives the typical transmission, which corresponds to the median outcome in an ensemble average: Here we have defined the localization length ξ such that D/ξ = − log (T i ) . The exact expression for the inverse localization length is given by Defining the reflection coefficient for the ith interface as R i = 1 − T i and Taylor expanding for small φ 0 gives This result implies that When φ 0 → 0, the localization length diverges, as mentioned above. We confirm this expression for ξ using a numerical simulation of the transfer matrix. The result is plotted in Fig.  4(a) [see also Fig. 2(a)], and is compared with the exact result derived in Eq. (26). When electrons are incident with a range of different ymomenta, the more physically-relevant quantity is the initialangle-averaged typical transmission, This quantity corresponds to the average intensity of Dirac electrons transmitted across a sample, given an a initial wave-packet of varying initial angle of incidence. As we discuss in Sec. III, this angle-averaged transmission is relevant for determining the electrical conductance. For simplicity we have assumed a uniform distribution of the initial angle. For large N, the leading order behavior of T typical can be found by extending the limits to infinity, which gives This equation constitutes one of the central results of our paper. While, for a given incident angle, the mapping to Anderson localization implies an exponential decrease of the transmission with distance, the divergence of the localization length for small angles gives an angle-averaged typical transmission that decays only as a slow power law, N −1/2 . This dependence is verified by a numerical simulation in Fig.  4(b). One implication of this result is that, for an incoming electron beam with a wide dispersion of incident angles, only a small fraction 24 πw 2 N is transmitted to a depth N. These transmitted electrons are narrowly collimated to have nearnormal transmission, and this collimation becomes increasingly tight as N is increased.
The statistical properties of the transmission probability can be understood by noting that, for a fixed incident angle φ 0 , there is an exact analogy between our problem and the problem of Anderson localization in a disordered 1D wire. The resulting distribution of T N is log-normal for sufficiently large N [37] It is also possible to derive a more general differential equation that describes P (T N ) at not-too-large N. The corresponding calculations are discussed in Appendix A and the references therein. Twisted bilayer graphene consists of stacked sheets of graphene with a rotational mismatch between the individual layers. Twisting creates a moiré superlattice of periodicity, which for small twist angles can be more than two orders of magnitude larger than the lattice constant of monolayer graphene. The corresponding moiré brillouin zone (MBZ) is relatively small. The MBZ is depicted schematically in Fig. 5, along with the brillouin zones (BZs) of individual graphene layers. Coupling between the two twisted layers leads to a hybridized dispersion relation with a renormalized Fermi velocity, which for certain "magic angles" can be as dramatic as reducing v F to a small percentage of its bare value [18,19,38]. Near these magic angles the Fermi velocity depends sensitively on the twist angle, so that even a mild twist angle disorder corresponds to a large velocity disorder. The twist angle also changes the size of the MBZ and shifts the Dirac points in momentum space, which is evident from Fig. 5. As explained by Padhi et al. [17], when the Dirac points shift in momentum space across a domain wall it eliminates the possibility of Klein tunneling for lowenergy electrons. Reference 17 focuses on this effect and neglects the change in velocity across the domain wall.
In fact, the question of transmission across a twist angle domain wall is crucially determined by the details of how the twist angle is altered. For example, if one of the two layers is fixed while the other is free to twist, then one of the two Dirac points of the MBZ also remains fixed in momentum space. (As illustrated in Fig. 5, only one of the two larger BZs is rotated in this case, so that only one of the two Dirac points in the MBZ is shifted.) Correspondingly, this situation produces Klein tunneling at one of the two nodes of the MBZ. If this situation persists for every twist angle domain (i.e., the same is layer always twisting), then the results from the previous section remain valid for one of the two valleys, including the diverging localization length at φ 0 = 0. On the other hand, the Dirac point corresponding to the other valley wanders from one region of the sample to another, and, as we show below, the maximum of transmission for this valley occurs at a different, nonzero angle. This situation implies a disorder-induced valley filtration in which electrons from only one of the two valleys are transmitted at normal incidence.
A more general scenario would involve both layers twisting from one domain to another, with a well-defined mean global twist angle. In the remainder of this section we assume that all twist angle changes happen symmetrically (as in Ref. 17), and we describe the Hamiltonian for TBG reported in Ref. 20.
Let v F be the Fermi velocity of pristine graphene, and K be the difference in momentum between its two Dirac points. Let us define the mean global twist angle of the TBG system to be θ, and we will denote the mean global position of the K point of the MBZ as K θ (and similarly the K point of the MBZ as K θ ). The length of the MBZ is given by |∆K θ | = 2K sin θ 2 , which corresponds to the separation between K θ , and K θ . Using the continuum description [20], the renormalized Fermi velocity is v 0 = 2 v 2 F |K θ | /t ⊥ , wherẽ t ⊥ is the interlayer coupling strength between the individual layers.
Now consider an arbitrary domain with twist angle θ + , where is the fluctuating part of the twist angle which can be either positive or negative. To leading order in , the K θ+ point of the MBZ in the new region is located at K θ + δ, where The new Fermi velocity is found to be The Hamiltonian in the vicinity of this Dirac point is given by where P = ( p x , p y ), and p x and p y are the momenta operators defined with respect to the Γ point of the monolayer graphene (which corresponds to the lab frame of reference, and is fixed regardless of twist angle, see Fig. 5). Let us define new momenta operators q x and q y , which measure the momenta with respect to the mean global K point of the TBG With this redefinition we have Comparing the above Hamiltonian with the Hamiltonian of massless Dirac electron in a magnetic field, we can identify that the shift in the Dirac point is equivalent to the application of vector potential A = −eδ/c produced by an external magnetic field [39] or associated with lattice defects/strain [40]. Translational symmetry along the y-direction suggests the following auxillary spinor as an ansatz Following a calculation similar to that of Sec. I, we arrive at the following solution for ψ A and ψ B : where and φ represents the angle between the y-momenta and the x-momenta of the electron with respect to the local K point. It is straightforward to show that q satisfies the following energy-momentum relationship: Now let us turn our attention into the case of a single domain wall located at In the analogy with a magnetic vector potential, the shift δ is equivalent to a sheet of magnetic field located at the domain wall and pointing perpendicular to the TBG, with magnitude It is interesting to see that transitional symmetry along the y-direction implies that the shift of the Dirac point in the xdirection is not a physically relevant perturbation. Just as the vector potential can affect physical observables only through its curl B = ∇ × A, so only the shift in the Dirac point in the direction parallel to the domain wall can enter the transmission coefficient.
Using the conservation of q y , it is possible to relate sin φ 0 and sin φ 1 : Defining This equation represents a modified version of Snell's law for an interface separating two twist angle domains. Here, ∆ 1 is a dimensionless quantity, which is a proxy for how much the K point of the MBZ shifts along the y-axis when passing across the domain wall. Due to the presence of this additional term −u 1 ∆ 1 in Snell's law, the symmetry between opposite incidences angles, φ 0 and −φ 0 , is broken, and electrons with opposite incidences angles are refracted by different amounts. When ∆ 1 1, so that the shift in the Dirac point is small compared to the electron's momentum, one can expect near-Klein tunneling. As we show below, the leading effect of small ∆ 1 is to shift the angle of perfect transmission away from zero. On the other hand, if ∆ 1 1, the transmitted electron states are gapped and exist only as an evanescent modes. In this latter case, transmission through a single barrier is exponentially small in D.
By constructing the transfer matrix similar to the previous section (see Appendix B), one can find the reflection probability of a single domain wall to be which is same as Eq. (18) except that the relation between φ 0 and φ 1 is modified [see Eq. (51)]. In the remainder of this section we will consider the case where ∆ 1 1, so that the electron momentum is relatively large. As we estimate below, in TBG this limit typically corresponds to electron energies E 0.1 meV. In this case of small ∆ 1 , one can make a Taylor expansion of Eq. (52) around φ 0 = 0 and ∆ 1 = 0, which gives This equation suggests that perfect transmission happens at φ 0 = (u 1 ∆ 1 )/(1 − u 1 ) rather than at normal incidence. Since this angle depends on the local Fermi velocity and the random vector potential, it varies across the sample. This situation is similar to that of the Brewster angle in optics, for which perfect transmission of polarized light across an interface occurs at an angle that depends on the local refractive indices. (Including higher-order terms in ∆ 1 in fact gives a finite reflection coefficient R ∼ u 4 1 ∆ 4 1 even at the optimal transmission angle. ) We can now calculate the statistical properties of the transmission through many domains of twist angle, as in the setup of Fig. 1. In addition to the random velocities u i associated with each domain, here we also account for the random shift of the Dirac point represented by the variables ∆ i . Since in TBG the randomness in u i and ∆ i arise from the same variation in twist angle and are therefore not independent, we consider the situation where λ and β are dimensionless experimental parameters.
(For a general situation where the Fermi velocity disorder and random energy shift are independent, refer to Appendix C.) We assume a uniform distribution for the twist angle variation i : Before proceeding to discuss the transmission, it is worth pausing to make an estimate of the typical numerical values of the parameters λ, β and µ. In TBG the first magic angle is ∼ 1.1 • and current experiments show a typical variation in twist angle of ∼ 10% [13], which gives µ ≈ 0.004. Within the continuum description used here [20], β = cot(θ/2)/2 [see Eq. (36)], so that using θ ∼ 1.1 • gives β ≈ 50. Thus, in typical TBG samples near the magic angle, the dimensionless random energy shift ∆ i 1 for all electron energies above 0.1 meV, as mentioned above. For example, if ∆ i ∼ 0.01 − 0.1, then λ is in the range 1 -25 for µ = 0.004.
In order to calculate the inverse localization length, D/ξ ≈ R i , we first Taylor expand R i in terms of i , i+1 and φ 0 and then proceed to take the expectation value over the random twist angle. This calculation gives where φ min = λβ 2λ 2 + β 2 and R min = (57) The finite reflection probability R min > 0 arises because the Brewster-like angle associated with transmission varies from one interface to another, and therefore on average there is no Klein tunneling. Instead, the disorder average gives a minimum of reflection at a nonzero angle φ min . Using the experimental parameters presented earlier in this section, we get an estimate of R min that is extremely small, in the range 10 −6 − 10 −9 . This small value of R min implies a very long localization length for electrons incident at the angle of minimum reflection, φ min . Correspondingly, the collimation of a dispersed electron beam is centered around a value φ min that depends on the electron energy E ∝ 1/λ [see Fig. 2(b)]. This result implies that twist angle disorder may be used not only to collimate electrons but also to filter them in energy, since the direction of the transmitted electron beam depends on its energy.
As in the previous section, the value of the transmission T N is log-normal distributed across an ensemble of samples [see Eq. (34)], with the mean value of log (T N ) given by N times Eq. (56). The analytical prediction for the mean of log (T N ) is plotted in Fig. 6(a), along with numerical results obtained using the transfer matrix approach. Using Eq. (56), we can calculate the angle-averaged typical transmission [see also Eq. (32)], which gives This expression is plotted in the Fig. 6(b) along with numerical results.

III. CONDUCTANCE AND FANO FACTOR
In the preceding sections we calculated the typical transmission coefficient T typical (φ 0 ) for electrons as a function of the initial incident angle φ 0 , considering both Fermi velocity disorder and a more general model of twist angle disorder in TBG. In this section we briefly discuss the implications of these results for the electrical conductance and the Fano factor.
According to the Landauer formula, the electrical conductance of a sample with width W is given by [35] where G 0 is given by the ballistic conductance 4e 2 h × W. For the small angles φ 0 that dominate T typical at large N, cos φ 0 ≈ 1. Hence, the conductance is effectively proportional to the angle-averaged transmission coefficient T typical , and G inherits the same power-law decay [see Eq. (33)], G ∝ 1/ √ N. One can also calculate the conductivity of the sample, given by Thus, σ ∼ √ N. A similar growth of the conductivity with the square root of system size was predicted by Ref. 24 for the case of massless Dirac electrons with uncorrelated, onedimensional scalar disorder potential.
Another experimentally relevant quantity is the Fano factor F associated with the shot noise, defined by The average of log (T N ) is plotted as a function of N for three different initial incidence angles: φ 0 = 0, φ min , and 0.20. The localization length increases when φ 0 increases from 0, reaches a maximum at φ 0 = φ min and then decreases again. (b) The average of the typical transmission over incidence angle is plotted as a function of log (N). The analytical curve is plotted using Eq. (58), which corresponds to the N → ∞ limit, and the numerical data is generated using the transfer matrix in Eq. (B1). The parameters used for the simulation are: µ = 0.004, λ = 4 and β = 50.
For large N, where the transmission is similarly dominated by small incident angles φ 0 , we arrive at This same value of the Fano factor has been derived for a smooth ballistic p − n junction in monolayer graphene by Cheianov et al. [41], where it arises from a similar dependence of the transmission on the incident angle.

IV. SUMMARY AND CONCLUSION
In this paper, inspired by the recent experimental observation of random domains of twist angle in TBG, we have considered the effect of Fermi velocity disorder on the transport of massless Dirac fermions. We have focused in particular on the model of quasi-one-dimensional disorder depicted in Fig. 1, which has a close analogy with the physics of a "transparent mirror" in optics [5], which itself can be mapped onto Anderson localization. Unlike the optical problem, however, the localization length in our problem diverges as the angle of incidence approaches zero [Eq. 31], and for such angles the material is transparent due to Klein tunnelling. Since the inverse localization length has a quadratic dependence on the incidence angle up to leading order, the angle-averaged typical transmission of the system decays with the system size in a power-law manner: (33)]. These results have a direct relationship to the conductance of the sample, which follows the same power law behavior [Eq. (59)]. The corresponding Fano factor is 1 − 1/ √ 2. Whether these results apply directly to twist angle disorder in TBG depends, in principle, on the details of how the twist angle varies from one domain to another. If one of the two layers is perfectly fixed while the other is allowed to twist in a spatially random manner, then Klein tunneling persists for only of the two valleys, leading to a disorder-induced valley filtration in which only one of the two valleys is transmitted at normal incidence. On the other hand, in the more realistic case where both layers undergo random twisting, one must consider the spatial modulation of both the Fermi velocity and a random gauge field that describes the shifting of Dirac points in momentum space. This gauge field leads to a modification of Snell's law [Eq. (51)] and a near-perfect transmission of electrons at a non-zero incidence angle, similar to the Brewster angle in optics. Since the Brewster-like angle varies from one domain wall to another, the overall reflection probability has a nonzero minimum [Eq. (56)], which occurs at non-zero incidence. Making numerical estimates based on recent TBG experiments, we found this minimum reflection coefficient to be very small, which suggests essentially perfect transmission at a nonzero angle. Thus, TBG may also perform disorder-induced collimation of incoming electrons, but to a nontrivial angle that depends on the energy of electrons.
It is also worth pointing out that the modified Snell's law in TBG [Eq. (51)] implies that for suitable values of u 1 and ∆ 1 , it is possible to have negative refraction. This situation gives rise to Veselago lensing, as in p-n junctions in monolayer graphene [6].
Of course, the model considered in this work is only a highly simplified case, for which the disorder is onedimensional. The corresponding translational symmetry in the y-direction has allowed us to use a simple transfer matrix approach and to map the problem onto one-dimensional Anderson localization. A more realistic model with random two-dimensional domains does not admit a similarly simple transfer matrix description. Still, the results presented here, in particular the phenomena of disorder-induced valley filtration and collimation, may provide inspiration for new directions in Dirac fermion optics based on tailored twist patterns. Here we provide more details about the statistical properties of the transmission probability. In particular, we can obtain the variance of log (T N ). If we define ∆ log (T N ) as log T N+1 T N , then by Eq. (23) where Φ NN+1 represents the phase accumulated during one complete internal reflection path between interface N and interface N + 1. Thus, for small incidence angle, R 1 1, we have to leading order in R 1 . Averaging this expression over random velocities, and employing the same random phase approximation used in Eq. (25), we arrive at This is same as the result obtained in Eq. (26). The variance of ∆ log (T N ) can be found to be var ∆ log (T N ) ≈ 2 R 1 .
This expression implies that the variance in transmission is dominated by the randomness in the phase Φ NN+1 across a domain, rather than in the values of T N or T 1 . For N interfaces, the corresponding variance in the transmission is var log (T N ) ≈ 2N R 1 .
A more elaborate calculation can be used to show that there is a correction of −π 2 /3 to the variance in Eq. (A5) [42]. In Fig. 7 we plot the variance calculated from numerical simulation along with this analytical prediction. Equation (A5) also establishes that the relative fluctuation in log (T N ) is var log (T N ) This decay of the relative fluctuation as the number of interfaces grows implies self-averaging, and confirms that log (T N ) is normally distributed in the limit of sufficiently large N. However, since T N is log-normal distributed, T N itself is not selfaveraging and the difference between T N and e log(T N ) becomes exponentially different at large N. One can extract the distribution for T N ; since the derivation is cumbersome, here we only present the results. We refer interested readers to the calculations in the lecture notes by Müller and Delande [37]. If we define ∆T N = T N+1 − T N , then to leading order in R 1 Following the previous averaging process we obtain