Intrinsic non-magnetic $\phi_0$ Josephson junctions in twisted bilayer graphene

Recent experiments have demonstrated the possibility to design highly controllable junctions on magic angle twisted bilayer graphene, enabling the test of its superconducting transport properties. We show that the presence of chiral pairing in such devices manifests in the appearance of an anomalous Josephson effect ($\phi_0$ behavior) even in the case of symmetric junctions and without requiring any magnetic materials or fields. Such behavior arises from the combination of chiral pairing and nontrivial topology of the twisted bilayer graphene band structure that can effectively break inversion symmetry. Moreover, we show that the $\phi_0$ effect could be experimentally enhanced and controlled by electrostatic tuning of the junction transmission properties.

Introduction.-Since its discovery in 1962 [1], the Josephson effect has emerged as one of the most remarkable manifestations of quantum coherence at the macroscopic scale with multiple technological applications [2].
Material platforms on which Josephson junctions (JJs) can be fabricated do not cease to grow, ranging from conventional or unconventional superconductor tunnel junctions [3] to hybrid nanostructures including exotic materials [4]. The discovery of superconducting phases in magic angle twisted bilayer graphene (MATBG) [5] provides an additional playground for the study of the Josephson effect.
In fact, the possibility to produce tunable junctions on this material by electrostatic gating has been recently demonstrated [6,7], with enabled phase control in ring shaped configurations [8], and JJs exhibiting unconventional Fraunhofer patterns have been found [9].
Many theoretical studies point to the possibility of chiral (d + id or p + ip) pairing symmetry arising as a result of electron correlations and the peculiar topology of flat bands in MATBG [17][18][19][20][21][22]. Furthermore, robust nematic behaviour was observed for a variety of twist angles in the superconducting phase of twisted bilayer systems [23][24][25][26], a phenomenon that has attracted recent theoretical interest [27][28][29][30][31][32]. The next key step in the field is thus to find reliable signatures that would help us to elucidate the pairing mechanism in MATBG.
In this Letter we address this issue showing that chiral pairing symmetry would manifest in an anomalous ϕ 0 behavior (a nonzero supercurrent in the absence of superconducting phase difference) in monolithic MATBG JJs. We demonstrate that this effect is a consequence of the nontrivial topology of the normal MATBG bands and thus cannot be captured by trivial models. We further show that the effect is enhanced for extended junctions and can be controlled by electrostatic gating of a middle normal region.
Before entering into the peculiarities of MATBG JJs, let us give a more general context on the topic of ϕ 0 junctions. In general, ϕ 0 behavior requires breaking time reversal symmetry (TRS) and inversion symmetry, having been predicted to appear in different systems by the combined effect of spin-orbit interactions and magnetic fields [33][34][35][36][37][38][39][40][41] or in junctions through magnetic metals with broken inversion symmetry [42][43][44][45]. In all these proposals the broken symmetries are linked to the spin degree of freedom and we could refer to them as "magnetic"-ϕ 0 junctions. In the case of MATBG the valley degree of freedom comes into play, providing new possibilities. In Ref. 46 it has been suggested that ϕ 0 behavior could appear in MATBG JJs with conventional s-wave pairing provided that the junction is established through a "valley polarized" region.
This mechanism has also been proposed to explain recent observations of anomalous Josephson effect in twisted trilayer JJs [47]. By contrast, we here analyze direct junctions between two MATBG regions with chiral pairing. An important aspect of MATBG which guides our study is that the large size of its moiré pattern (>10 nm) would allow transport experiments on junctions along well defined directions on the moiré lattice. We show that, only when the nontrivial topology of the MATBG is taken into account, ϕ 0 behavior can appear spontaneously when the junction is defined along certain directions on the moiré pattern that break translation invariance, or along any orientation when graphene's sublattice symmetry is broken. We argue that these effects would allow one to distinguish among possible pairing symmetries in MATBG.
Qualitative description.-We consider a junction defined on a bulk MATBG sample as schematically illustrated in Fig. 1(a). External gates can be used to tune the doping level on the two sides of the junction independently and each side is characterized by the presence of a generic chiral pairing. The barrier height, and therefore the junction transmission, can also be externally controlled. The junction is placed on a loop that enables control over the superconducting phase difference ϕ, determining the junction current-phase At band fillings close to the van Hove singularities (VHSs), where superconductivity is observed, MATBG is characterized by a normal Fermi surface exhibiting trigonal warping [48,49]. Typical shapes of these surfaces on the two valleys (denoted by K and K ′ ) are shown in Fig. 1(b). The main effect of trigonal warping is to break the intra-valley inversion symmetry along certain directions, i.e., E(k ⊥ ) ̸ = E(−k ⊥ ), where k ⊥ denotes the wavevectors perpendicular to the junction interface. The Fermi velocities v F at the Fermi points for a given parallel momentum k ∥ thus differ. This is the case, for instance, for the junctions defined parallel to the line joining the moiré minivalleysK −K ′ (i.e., k ∥ = k y ), which we call "zigzag" (ZZ) in analogy with pristine graphene. On the contrary, junctions defined perpendicular to theK −K ′ line (i.e., k ∥ = k x ), called "armchair" (AC), are not affected by trigonal warping in the normal state [50]. However, as we discuss below, such symmetry can be broken even for AC junctions due to the combined effect of normal state topology and pairing properties. As illustrated in Fig. 1(c), these junction orientations correspond in real space to lines joining nearest neighbors (ZZ) or next nearest neighbors (AC) "sites" in the moiré pattern, which can be associated with the charge centers.
From this definition we observe that an AC junction breaks the moiré translation invariance along the interface.
Periodicity along the junction interface (i.e., along k ∥ = k x ) is recovered by cell duplication, which corresponds to folding the Brillouin zone into a rectangular one as shown in Fig. 1(d). Consequently, in the superconducting state the Andreev bound states (ABSs) formed at such interface result from coupling scattering states of different type on both sides, which might differ by a topological phase.
While this asymmetry is a necessary condition, to obtain a net ϕ 0 behavior this phase difference should survive after k ∥ and valley integration.
We find that spontaneous valley supercurrents I K J (ϕ = 0) ̸ = 0 appear in the AC case due to the combined effect of chiral pairing and nontrivial band topology, and in the ZZ due to Fermi-surface warping. However, while the ZZ currents on opposite valleys cancel each other [i.e., I K J (0) = −I K ′ J (0), see Fig. 1(e)], for AC we have Consequently, the valley currents do not compensate leading to a net ϕ 0 -junction behavior, see Fig. 1(f). As we show below, breaking graphene's sublattice symmetry yields an anomalous CPR also for junctions defined on a ZZ direction only in models where nontrivial band topology is taken into account.
Our main result, which we expose in detail below, is thus that the presence of chiral pairing on a nontrivial MATBG band manifests as an anomalous Josephson effect. This is always the case in AC junctions, due to the naturally broken moiré translation invariance at the interface, and can be induced on ZZ junctions by breaking graphene's sublattice symmetry.
CPR calculations for MATBG junctions.-To give a quantitative estimation of the effect discussed above, we consider a six band model (6BM) for MATBG junctions which accounts for both the appropriate flat bands dispersion and their nontrivial (fragile) topology [51]. Our calculations are based on the recursive Green functions method introduced in Refs. 52-54, which we extend here to the superconducting case, see the Supplemental Material (SM) [55] where we also include the pairing implementation for the trivial two band model (2BM) [49]. To introduce chiral superconductivity in these models we project the pairing order parameter onto the directional p ± orbitals which provide the main contribution to the flat bands below the VHS [55]. For simplicity, we only consider here intervalley, spin-singlet chiral d-wave superconductivity with zero center-of-mass momentum. From the irreducible representations of the crystal symmetry group [56] in the triangular lattice we build up the form factors for the pairing where the chiral pairing is the superposition ∆ d+id ′ = ∆ d x 2 −y 2 +i∆ dxy , L x,y are the orthogonal lattice vectors in the doubled unit cell, and we set the value of ∆ as roughly one order of magnitude smaller than the bandwidth. Some other configurations (e.g., chiral p-wave with S z = 0 or superpositions of nematic chiral orders with local s-wave) could be implemented using the same approach, but have minor effects on the CPR. Notice that chiral superconductivity is a necessary ingredient to obtain ϕ 0 behavior as it breaks TRS. However, not all TRSbreaking pairings lead to ϕ 0 effect, e.g., nodal intravalley pairings never show uncompensated valley currents for any junction configuration studied in this work [55]. In this sense, ϕ 0 behaviour is intrinsically linked only to chiral superconductivity. Along with pairing, we consider the effect of a sublattice symmetry breaking perturbation (δ S ) acting over p ± that can be produced by alignment with the hBN substrate [55]. We note that the nontrivial topology of the bands in the 6BM is already reflected in an anisotropic response of the superfluid weight [55]. The relevant effect of remote bands for superconductivity in MATBG has been discussed in several theoretical works [27,[57][58][59][60][61].
The junction is defined as a smooth barrier at the moiré scale for which k ∥ is a good quantum number. Its effective transmission is controlled by a coefficient τ such that for a symmetric junction (SS) with the same chemical potential on both leads, µ L = µ R = µ 0 , we recover the bulk system at zero phase bias and τ = 1. We also consider asymmetric junctions (SS ′ ) in which µ L ̸ = µ R and symmetric junctions through a finite-sized, heavilydoped normal region (SN S) where |µ L | = |µ R | ≪ |µ c | and τ = 1. We first discuss the case of AC junctions. Figure 2(c) shows the evolution of the zero-phase supercurrent I J (0) for τ = 0.8 as a function of the chemical potential, both in the case of equal (µ L = µ R ) and unequal (µ L fixed, varying µ R ) doping levels for the 6BM. Since the critical currents for all cases are of the same order, we normalize I J (0) to the critical current I 0 c for µ L = µ R = −1 meV. The red dashed lines indicate the VHS positions.
While the most striking result is the presence of ϕ 0 behavior for a symmetric junction without external fields, the computed |I J (0)/I 0 c | values are relatively small (less than 10%) in the range of doping levels where superconductivity in MATBG is expected to appear. We observe similar values of I J (0) for symmetric junctions with p or n doping close to the VHSs, see Fig. 2(c). The ϕ 0 effect is generally larger for asymmetric junctions of n-p type, i.e., with µ L < 0 and µ R > 0. In contrast to the symmetric case, in this situation we observe no changes in the sign of I J (0) as a function of the doping level. Finally for both symmetric and asymmetric junctions we find large non-vanishing |I J (0)/I c | values close to charge neutrality µ = 0. We notice, however, that experimental samples show no robust superconducting dome at those fillings and that the superconducting state could have different properties close to charge neutrality [11,62].
Having demonstrated the anomalous Josephson effect of short symmetric and asymmetric junctions with chiral pairings, we now show that the ϕ 0 effect can be enhanced and tuned by a middle normal region in the SN S configuration. We thus assume symmetric junctions including a finite normal central region with large doping beyond the flat bands, so that the barrier transmission is limited by the chemical potential mismatch. We show in Fig. 3 that it is possible to tune the |I J (0)/I c | ratio by varying the doping levels (a) or the length (b) of the central region. Notice that doping the central region at the dispersive bands or increasing its length reduces the critical current by an order of magnitude with respect to the cases analyzed in Fig. 2. Finally, we analyze ZZ junctions in Fig. 4. As discussed above, trigonal warping leads to compensated valley currents at zero phase bias. However, when sublattice symmetry is broken, the nontrivial topology of MATBG captured by the 6BM leads to uncompensated currents and thus ϕ 0 behavior. This behavior is not observed in the trivial 2BM [55]. Similarly to AC junctions, we find larger |I J (0)/I c | values for asymmetric junctions of n-p type. We also observe a change in sign of I J (0) through the p-p' to p-n transition. The size of |I J (0)/I c | also increases as the sublattice symmetry breaking parameter δ S value is increased [inset of Fig. 4(c)]. It is worth mentioning that we observe asymmetries between the negative and positive critical currents (superconducting diode effect [38,39,63]) in both AC and ZZ orientations for asymmetric SS ′ and SN S junctions. This effect is observed when the ϕ 0 behaviour is enhanced and could thus be electrically controlled by tuning the chemical potential and effective transmission. Let us further comment that if, additionally, valley polarization is included, e.g., by means of TRS breaking in the parent state as in Ref. [46], one would observe ϕ 0 behavior for all orientations.
Minimal model.-A scattering theory that confirms and gives further insights to the previous results can be built by linearizing the bulk Hamiltonian with respect to k ⊥ around each Fermi point for a given k ∥ [55]. For the 6BM, Andreev reflection coefficients for electron to hole processes (and their time reversed) acquire different phases θ(k ∥ ) and θ ′ (k ∥ ). Applying matching conditions at the interface corresponding to a single channel junction with transmission T k = T (k ∥ ), leads to the following set of Andreev bound states [55] ϵ(k ∥ ) , in agreement with the full numerical results. By contrast, the warping distortion of the Fermi surface in the ZZ case results in θ(k ∥ ) ̸ = θ ′ (k ∥ ), even for the non-topological 2BM [55]. However, in this orientation we find that θ K (k ∥ ) = −θ K ′ (−k ∥ ) (same for θ ′ ) and, consequently, I K J (0) = −I K ′ J (0). The ZZ valley currents are thus compensated unless sublattice symmetry is broken and fragile topology taken into account.

Conclusions:
We have shown that chiral pairing in MATBG junctions would manifest in the appearance of a nonmagnetic ϕ 0 behavior, stemming from the nontrivial topology of the MATBG bands. This effect can then be used to distinguish chiral superconductivity from other mechanisms such as nodal pairing. Moreover, we illustrated how the effect is enhanced for extended junctions and controlled by electrostatic gating of a middle heavily doped normal region. Although we have focused on d + id case, other chiral symmetries like p + ip or combinations of these symmetries with s-wave pairing, would show similar behavior. The orientation sensitivity of the ϕ 0 effect could help to distinguish this different orbital character in an actual experiment. Furthermore, it would allow us to distinguishing whether the anomalous behavior is due to a broken valley symmetry parent state or due to configurations involving symmetric contributions from both valleys and chiral pairing.
We thank F.S. Bergeret, R. Egger, and W.J. Herrera for valuable comments and discussions. This work was supported by the Agencia Estatal de Investigación project No. PID2020-117992GA-I00 and No. PID2020-117671GB-I00, the Spanish CM " Talento  The implementation of these models for normal transport calculations has already been discussed by us in Ref. 53. Here, we describe the necessary extensions of these models to account for superconducting transport.
In order to obtain the transport properties in defined armchair (AC) and zigzag (ZZ) orthogonal boundaries we have to consider the cell doubling of the unit cell [54] with lattice vectors L x = √ 3L m and L y = L m , where L m = a/(2 sin θ/2) ≈ 13 nm is the moiré lattice length given in terms of the graphene lattice parameter a and the twist angle θ, see Concerning the implementation of superconductivity in MATBG, we have focused on configurations involving symmetric contributions from both valleys, which are the most likely ones according to many previous theoretical studies for the SC domes around filling ν = −2 (see Refs. [12,15,64]). Therefore, we exclude further effects coming from possible flavor resets when varying the band filling of the flat bands. Regarding the rest of degrees of freedom of the system, we adopt the rather natural spin projection S z = 0 intervalley chiral superconductivity with zero net momentum Cooper pairs [16]. Despite the fact that the orbital character of the pairing cannot be determined by existent experimental data, we will focus on d + id ′ wave as a typical example of time reversal symmetry (TRS) breaking topological superconductivity. In any case, p + ip ′ with S z = 0 spin projection shows similar properties to the chiral d-wave, featuring the aforementioned ϕ 0 behaviour (not shown in this work).
We thus characterize the bulk MATBG by a Bogoliubov-de Gennes Hamiltonian (BdG), where the electron (hole) HamiltonianĤ e ( ⃗ k) [Ĥ h ( ⃗ k)] and pairing∆( ⃗ k) can take into account the multiband, fragile The intravalley term is used to rule out nodal, TRS-breaking (non chiral) pairing mechanism as a source of ϕ 0 behaviour.

Superconducting 6BM Hamiltonian with unit cell doubling
The normal state of the 6BM is defined in a basis of p z and p ± -orbitals in a triangular lattice, (τ, p z ) and (τ, p ± ), respectively, and three s-orbitals in a Kagome lattice (κ, s), see Fig. S 1(a). The local fermion operators in the doubled unit cell are defined asΨ µ α = (τ pz,µτp+,µτp−,µκ The low energy physics associated to the flat bands in the 6BM are described fundamentally by the directional p ± -orbitals in the triangular lattice (τ, p ± ) [51]. Therefore, we project the pairing order parameter over these orbitals only. We restrict ourselves to the simplest case in which there is no sublattice structure. The pairing wavefunctions that preserve the global symmetry of the lattice are obtained from the irreducible representations of the crystal symmetry group [56] in the triangular lattice, Recent experiments [11] based on Andreev reflection measurements point towards ∆ ≈ 0.3 meV. In our calculations we set the value ∆ = 0.1 meV, one order of magnitude smaller than the bandwidth. We project this form factors in the (τ, p ± ) subspace expanded to take into account the moiré superlattice degree of freedom in the basisΨ ± = (Ψ 1 p−,σ ) T and σ =↑↓ being the implicit spin degree of freedom. We observe the following intraorbital contributions for the nodal terms d = d x 2 −y 2 and d ′ = d xy , where ϕ 10 = e −ikxLx and ϕ 01 = e −ikyLy , with ϕī j = ϕ † ij . If we explicitly show the total electron and hole contributions of the Hamiltonian in the moiré superlattice space we obtain Experiments until this date have shown an incompatibility between superconductivity and full hBN alignment, but symmetry breaking superconducting domes have been observed close to charge neutrality [11,62]. For the purpose of this work, graphene's sublattice symmetry breaking can act as an example of neither magnetic, nor valley, polarization mechanism to obtain uncompensated valley currents for the ZZ case.

Superconducting 2BM Hamiltonian with unit cell doubling
We use the topologically trivial 2BM as a test to compare the features associated primarily to the Fermi surface properties and the ones in which the fragile topology is involved. We use the same pairing wavefunction projected in the triangular lattice, as for the 6BM. However, in the hexagonal lattice of the 2BM, we only consider the next nearest neighbours (nnn) terms to expand the order parameter, neglecting the nearest neighbours (nn) ones, so we can directly compare the effects of the Fermi surface over the exact same pairing wavefunction. We implement the chiral d-wave using the following spinorsΨ = (Ψ 1 where h 0 = t 2 ϕ 10 + t * 2 ϕ1 0 , (S 8a) Finally we can break sublattice symmetry adding the perturbationV 11 =V 22 = δ Sτ AB z , whereτ is a Pauli matrix.
The pairing contributions for the nodal d and d ′ have the form∆ Finally the chiral order parameter can be expressed aŝ

Analysis of superfluid response
To observe the inherited topological effects in the pairing, we now analyze the superfluid weight. The superfluid weight gives the size of the superfluid current for a given phase gradient in the bulk, and it is connected to central phenomenology of the superconducting phase, like the Meissner effect.
Recent works relate the superfluid weight in multiband systems featuring superconducting flat bands with topological properties like the Berry curvature [27,[57][58][59][60][61]. Following Ref. 27, we study the superfluid weight in the doubled unit cell defined via the static Meissner effect for local and nonlocal interactions where µ, ν = {x, y} are the different components of the superfluid weight tensor marking the orientation of the partial derivatives of the momentum ∂ µ = ∂ kµ , V is the area of the sample in k-space, f (E) is the Fermi-Dirac distribution in the zero temperature limit, andĤ|Ψ n k ⟩ = E n k |Ψ n k ⟩ represent the eigenvalue problem for the n-th band of the BdG Hamiltonian at momentum k contained in the first Brillouin zone.
The D s xy component of the superfluid weight, Fig. S 2, shows nonzero values for the 6BM case. This feature is associated to spontaneous C 3 rotational symmetry breaking and thus unexpected nematic response for the 6BM, which is specially noticeable for the AC orientation. Notice that this effect is not compensated between valleys in contrast to the 2BM ZZ case. Previous microscopic tight-binding calculations of the superfluid weight showed anisotropies in conventional MATBG bulk when the pairing is nonlocal [27]. This phenomenon can be explained by the nontrivial multi-orbital character of the 6BM that is captured in the quantum geometric tensor and, consequently, in the superfluid weight [60]. The value of the anisotropic response is proportional to ∆ and equivalent for local and nonlocal pairing wavefunctions.
The anisotropic D xy component demonstrates the presence of inherited unconventional properties associated to the pairing in the 6BM bulk when the doubled unit cell is considered. This anomalous response is specially guaranteed for the AC orientation and exemplifies the different properties of the 2BM and 6BM even in the superconducting phase, despite the fact that both models have the same Chern number when chiral superconductivity is taken into account.

TRANSPORT CALCULATIONS
The phase-biased equilibrium Josephson current can be obtained using Keldysh Green functions (GFs) [52] as where Ω k ∥ = 2π/L ∥ accounts for the integration limits, the trace tr N and Pauli matrixτ z act over Nambu space, Σ LR is the coupling between the leads andĜ +− are the Keldish components of the junction GF. At equilibrium we havê ) where j, j ′ = L, R are written in terms of the advanced GF (from now on we omit the superindex A). We obtain the edge GFsĜ jj ′ from the full junction GF whereĜ L/R is the boundary GF of the left (right) side associated to the superconducting MATBG lead.
To compute these boundary GFs, we use the recursive method described in Ref. 54 (notice that in the doubled cell space the 6BM is described up to nn and the 2BM up to nnn). Experimental setups typically consist on a TBG slab with several electrodes but no constrictions or abrupt junctions (e.g., other metal oxides, etc.). Our description must recover a perfect bulk when no phase bias is applied and the chemical potential along the sample is homogeneous. To do so, we use the following general self-energy structure for the junctionŝ where τ is an effective transmission coefficient of the junction andT LR (k ∥ ) are the nonlocal contributions of the BdG Hamiltonian such that for a symmetric junction and τ = 1 we recover a pristine bulk. With this general structure, and taking advantage of the recursive Green function method, we define three different types of junction. First, a symmetric SS junction, in whicĥ H L =Ĥ R , defined by a phase bias and a non-perfect transmission τ < 1. Equivalently, but varying the filling on one of the leads with respect to the other, we obtain an asymmetric SS ′ junction where µ L ̸ = µ R . Finally, we can compute the more realistic junction in which τ = 1 and the phase bias is applied along a finite normal region of a few moiré wavelengths. This corresponds to a SN S junction in which the normal region is heavily doped |µ L | = |µ R | ≪ |µ c |. In Fig. S 3 we show results for the subgap Andreev spectrum ρ(E, ϕ) obtained for the two different models for the case of a SS junction. The subgap spectrum is dramatically different for the two models: it satisfies ρ(E, ϕ) = ρ(E, π − ϕ) for the 2BM while this symmetry is broken in the 6BM. In particular, the phase shift (ϕ 0 ) for Andreev bound states (ABSs) in the 6BM varies with k ∥ , featuring maximum phase displacement atX/2. The ABSs satisfy that ϵ K (k ∥ , ϕ) = ϵ K ′ (−k ∥ , ϕ), enabling the appearance of a net ϕ 0 behaviour despite the different contributions from each transport channel at each valley.

Linearized scattering approach
To obtain the CPR from a scattering perspective, one needs to determine the Andreev reflection coefficients at an ideal interface between a normal and a superconducting MATBG region, assuming conservation of the momentum parallel to the junction interface k ∥ . We start by linearizing Eq. (S 2) around the i-th Fermi point with respect to the incident momentum k i ⊥ , i.e., The scattering states |Ψ i E ⟩ would then satisfy the equation for the correction of the momentum with respect to the Fermi point, δk i ⊥ , whereÂ is the adjoint matrix ofĤ i 1 . We can further simplify the problem by projecting Eq. (S 17) into the subspace spanned by the states |Ψ i ν;α ⟩, corresponding to the dominant electron (ν = e) and hole (ν = h) states at the i-th Fermi point. The label α denotes additional degrees of freedom which could be coupled by scattering at the junction interface. Finally, we project the problem into a 2×2 Nambu space, namely, where E is the incident energy of the quasiparticle and ∂ k = ∂ k ⊥ . The equation for δk i ⊥ acquires the eigenvalue problem structure with δk i ⊥ (E) = λ/ detÂ i . The eigenvalues of δk i ⊥ (E), that is, the corrections to the Fermi momentum in the direction perpendicular to the junction, are The eigenvalues δk i ± (E) (with associated eigenstates |Ψ i ± (E)⟩) become complex at subgap energies. Consequently, we can extract two Andreev reflection coefficients for each Fermi point, While for a conventional single-band superconductor the Andreev coefficients on opposite Fermi points exhibit exactly opposite phases, for general MATBG models with chiral pairing this inversion symmetry can be broken, even when it is preserved in the normal state as in the case of the AC orientation. The 2BM model satisfies that the eigenvalues δk i ⊥ (E) for the chiral d-wave come in complex conjugated pairs and thus we find that ∂ k ∆∆ † = (∂ k ∆ † ∆) † for E = 0. The nontrivial behavior of the 6BM is revealed in the phase of the Andreev reflection coefficients obtained from the eigenvectors z i In general, the z i ± coefficients on opposite Fermi points satisfy the following approximate ansatz for AC junctions based on the previous numerical results, with where∆ k is an effective gap which depends on k ∥ and ρ ≳ 1 a dimensionless parameter. In the AC configuration we have two symmetric Fermi points k 1,2 ⊥ ≡ ∓k F . The boundary modes take the form Notice that right and left moving solutions have a different phase of the Andreev coefficients for each Fermi point, θ k = θ(k ∥ ) and θ ′ k = θ ′ (k ∥ ), thus breaking inversion symmetry.
The boundary condition at the interface, imposing a single channel matching with effective transmission T (k ∥ ) = T k and fulfilling R k + T k = 1, takes the form where τ k = √ T k , r k = √ R k , and σ z is a Pauli matrix acting in Nambu space. The equation for the Andreev bound states is thus obtained imposing that the determinant of the matrix of coefficients (a, b, c, d) is zero. In general, the solutions of this equation, E 0 , are complex numbers with both ϵ(k ∥ ) and Γ(k ∥ ) real. The energy (real part) of the bound state equation adopts the simple form (S 27) withθ k = θ k + θ ′ k and δθ k = θ k − θ ′ k . But the imaginary part Γ(k ∥ ) is, in general, nonzero. For T k < 1, however, we find that the imaginary part vanishes and we have real Andreev bound states in the limit ρ → 1. Equation (S 27) thus shows that the ϕ 0 effect originates from the symmetry breaking induced by the different phases obtained from the eigenvalue problem in the 6BM case.
For the AC junction, the 2BM is recovered when θ ′ k = θ k so the phase shift vanishes, δθ k = 0. In this case, the trivial phases accumulated only produce a renormalization of the effective transmission of the junction. In the AC 6BM, we observe that θ K (k ∥ ) = θ K ′ (−k ∥ ) and, equivalently, for θ ′ (see Fig. S 4). Consequently, we observe different ϕ 0 contributions for each k ∥ that integrate to the same ϕ 0 -shifted total valley currents with equal magnitude and sign.
For ZZ junctions, in both models we observe compensated valley currents I K J (0) = −I K ′ J (0) for E ̸ = 0, in agreement with the full numerical results. Only in the 6BM, and when sublattice symmetry is broken, we expect to observe uncompensated valley currents due to asymmetries between valleys induced in the eigenvalues of the scattering problem.
Finally, when intravalley pairing is considered, we only observe ϕ 0 behaviour for chiral symmetries and under the same conditions as studied in the intervalley case, in agreement with the full numerical results. By contrast, we observe compensated valley currents for asymmetric SS ′ junctions when nodal pairing is studied, as predicted by the accumulated phases in the Andreev coefficients analysis. These valley currents for both AC and ZZ orientations are always compensated I K J (0) = −I K ′ J (0), even when graphene's sublattice symmetry is broken. Consequently, in the case of nodal pairing other valley polarizing mechanism would be required to lift the symmetry and observe the ϕ 0 effect. Therefore, both the scattering and the full numerical analysis link the anomalous Josephson effect only to chiral pairing in MATBG. No other TRS-breaking symmetry, like nodal intravalley, would exhibit ϕ 0 behaviour in direct Josephson junctions. As we explain in the main text, this fact provides a direct way to distinguish between nodal and full-gapped chiral pairing in MATBG.