Incoherent tunneling and topological superconductivity in twisted cuprate bilayers

Twisting two monolayers of a high-$T_c$ cuprate superconductor can engender a chiral topological state with spontaneously broken time reversal symmetry $\mathcal{T}$. A crucial ingredient required for the emergence of a gapped topological phase is electron tunneling between the CuO$_2$ planes, whose explicit form (in an ideal clean sample) is dictated by the symmetry of the atomic orbitals. However, a large body of work on the interlayer transport in cuprates indicates importance of disorder-mediated incoherent tunneling which evades the symmetry constraints present in an idealized crystal. The latter arises even in the cleanest single-crystal samples through oxygen vacancies in layers separating the CuO$_2$ planes, introduced to achieve the hole doping necessary for superconductivity. Here we assess the influence of incoherent tunneling on the phase diagram of a twisted bilayer. We show that the model continues to support a fully gapped topological phase with broken $\mathcal{T}$, even in the limit of disorder-mediated interlayer coupling. Compared to the model with a constant, momentum conserving interlayer coupling, the extent of the topological phase around the 45$^\circ$ twist decreases with increasing incoherence, but remains robustly present for parameters likely relevant to Bi$_2$Sr$_2$CaCu$_2$O$_{8+\delta}$.


I. INTRODUCTION
Twisted 2D van der Waals materials have emerged as an elegant platform to engineer and study correlated quantum phases with unprecedented experimental control [1][2][3][4][5][6]. At certain magic angles, the electronic structure is deformed into exceedingly narrow bands in a moiré Brillouin zone. The small bandwidth yields pronounced correlation effects and generally makes multiple neighboring quantum phases experimentally accessible in a single fabricated device through application of electrostatic gating [3].
Recently, a related paradigm of twisted bilayer structures was introduced that does not rely on engineering of correlated flat bands, but can produce interesting new phases by combining known non-trivial properties of constituent monolayers [7,8]. It was shown that two cuprate monolayers, stacked at a twist angle θ, can give rise to a spontaneously time-reversal symmetry broken state by virtue of simple electron tunneling between the layers. Most notably, in a finite range of angles around the critical twist θ c = 45 • , the ground state of an otherwise nodal d-wave superconductor becomes fully gapped and acquires a finite Chern number. At exactly 45 • the gapped phase persists up to the native critical temperature of the cuprate monolayer, thus furnishing the first known example of a high-T c topological superconductor.
Pioneering experimental work on very thin twisted Bi 2 Sr 2 CaCu 2 O 8+δ (BSCCO) flakes succeeded in fabricating bilayers at various twist angles [9]. Measurements of the interlayer Josephson current, Fraunhofer interference patterns and half-integer Shapiro steps in samples close to the 45 • are suggestive of a T -broken phase [10]. Strong twist angle dependence of the critical current has been reported elsewhere [11].
It was later noted by Song, Zhang and Vishwanath [12] that twist angle and momentum dependence of the interlayer tunneling matrix element g k , arising from the symmetry of the copper active orbitals, can play an important role in the emergence of the T -broken phase. As originally argued by Xiang and Wheatley [13], the matrix element has the form g k = g 0 cos 2θ + g 1 µ k (θ/2)µ k (−θ/2), which we generalized here to a twisted bilayer geometry following [12]. The g 0 term represents the direct tunneling between copper d x 2 −y 2 orbitals while the g 1 term describes the 'oxygen-assisted' tunneling process with form factor µ k (θ) = [cos (R θ k x ) − cos (R θ k y )]/2 and R θ the rotation matrix. Crucially, the part of g k that remains non-vanishing at θ = 45 • contains the form factor µ k that suppresses tunneling at the nodes of the d-wave order parameter. While spontaneous T -breaking is still found to occur in this case, the ground state remains gapless and hence does not support the gapped topological phase with non-zero Chern number predicted in Ref. [7].
In the present work we consider the twisted bilayer problem within a family of incoherent tunneling models [14][15][16][17][18] in which the transfer of electrons between two adjacent CuO 2 layers is mediated by impurities that are inherently present in the otherwise inert 'spacer' layers. Such incoherent tunneling models have been shown to yield better agreement with experimentally measured caxis transport properties of nominally clean single-crystal cuprates than models where momentum is strictly conserved [19,20]. Because random impurities break all spatial symmetries of the system, the form of the interlayer coupling is required to respect the crystal symmetry constraints only on average. One may thus expect that incoherent tunneling models will evade the difficulties noted above and produce a fully gapped topological phase near θ = 45 • .
Based on a perturbative diagrammatic treatment within a simplified continuum model we show that the incoherent tunneling model indeed delivers the same phenomenology as the coherent model of Ref. [7] while respecting with the point group symmetries of the physical system. Importantly, we show that for sufficiently slowly varying disorder the ground state near θ = 45 • is gapped in the T -breaking region and topologically non-trivial. These results are then confirmed in a more realistic setting through a full numerical diagonalization of a lattice model with parameters chosen to reproduce the actual cuprate band structure in the vicinity of the Fermi level. The effect of interface inhomogeneity on Josephson effects in twisted bilayers was recently considered in Ref. [8] where it was found that sufficiently strong disorder can leave the system in a topologically trivial state around 45 • . This is consistent with our deductions.
The paper is structured as follows. After summarizing the origin of the T -broken phase in the language of group theory (Sec. II), in Sec. III we introduce a model of incoherent interlayer tunneling within a continuum framework and show that a fully-gapped T -broken phase emerges in the vicinity of θ = 45 • . The phase is topological as evidenced by chiral edge modes traversing the gap that appear in the disorder-averaged boundary spectral function. In Sec. IV, we supplement our continuum results with a lattice model that simultaneously captures the characteristic cuprate Fermi surface geometry, the moiré effects and disorder in interlayer coupling. Concluding remarks appear in Sec. V.

II. GROUP THEORETICAL DISCUSSION OF T -BREAKING IN TWISTED CUPRATES
The phenomenology of T -breaking in twisted cuprates can be captured by a two-component Landau-Ginzburg theory with complex order parameters Ψ 1 and Ψ 2 with a relative phase ϕ. The ϕ-dependent part of free energy is of the general form Time reversal symmetry will be spontaneously broken whenever the free energy develops two minima that are related by T : ϕ → −ϕ. The Josephson coupling term, proportional to cos ϕ, has only a single minimum at ϕ = 0 or π, depending on the sign of B. Presence the fourth order term proportional to C cos 2ϕ is therefore necessary to break T . Additionally, one must have C > 0, since otherwise the two minima of cos 2ϕ occur at ϕ = 0, π which map to themselves under T . In Ref. [7] it was argued that C is indeed positive based on microscopic mean-field calculations. We confirm that C remains positive in the case of incoherent interlayer coupling in Sec. III B.
Given C > 0, the fourth order term is minimized at ϕ = ±π/2. Then, T -breaking occurs as a consequence of the competition among the two terms in Eq. (2). Specifically, T will be broken when A special situation clearly arises if symmetry requires B to vanish; then Eq. (3) is guaranteed to be satisfied for any C > 0. Next we describe a set of symmetry requirements under which the coefficient B vanishes and the system is forced into the T -broken phase. The order parameters Ψ 1 , Ψ 2 transform according to irreducible representations (irreps) of the point group of the crystal. Two cases must be distinguished: (a) Ψ 1 and Ψ 2 transform under two different 1D irreps or (b) (Ψ 1 , Ψ 2 ) transform under a 2D irrep [21,22]. The latter case is considered a generic pathway to T -breaking that occurs immediately upon entering the SC phase. The former case generically yields two successive phase transitions with distinct critical temperatures, T c and T c , with T -breaking setting at the lower one T c . Note that T c can be zero or negative, in which case the T -broken phase is physically not accessible [21,23].
The point symmetries of an untwisted cuprate bilayer form the point group D 4h . Here, the d x 2 −y 2 and d xy order parameters transform according to the 1D irreps B 1g and B 2g , respectively. At arbitrary twist angle, inversion and mirror symmetries are broken and the point group reduces to D 4 with d-wave irreps B 1 and B 2 . Thus, given the d-wave nature of the order parameter in cuprates, only pathway (b) to T -breaking is possible and no definite symmetry-based arguments can be made.
Precisely at 45 • , however, the symmetry group is enlarged to the non-crystallographic point group D 4d which contains an additional 8-fold improper rotation S 8 of the quasicrystalline lattice. Most notably, among the irreps of D 4d only the 2D E 2 -irrep supports d-wave order. Thus, the Josephson coupling term −B |Ψ 1 Ψ 2 | cos ϕ must necessarily be absent at 45 • . This is so because it descends from the −B(Ψ 1 Ψ * 2 + c.c.) term in the free energy which is not invariant under S 8 : (Ψ 1 , Ψ 2 ) → (Ψ 2 , −Ψ 1 ). Thus, T -breaking can be viewed as a fundamental consequence of the point group at θ = 45 • . We summarize our key arguments as follows: Two-component order parameters that transform under a 2D irrep naturally break T . At 45 • twist angle, because the point group of the bilayer is D 4d , any d-wave order parameter must necessarily transform under a 2D irrep. Therefore, the superconducting state breaks T right below T c .
The phase diagram of twisted bilayer cuprates derived in Ref. [7] can then be understood from continuity arguments. It is expected that the T -breaking phase will not be limited to the exact 45 • -twist, but will extend to a range of twist angles in its vicinity. Since at twists slightly away from 45 • the order parameters transform under two 1D-irreps, two distinct critical temperatures are permitted and T -breaking will no longer coincide with the critical temperature T c of the spontaneous U (1)-symmetry breaking. This naturally leads to the wedgeshaped T -broken domain in the phase diagram explicitly computed in Ref. [10]. Our symmetry arguments will be manifest in a microscopic description of the bilayer system.

III. INCOHERENT TUNNELING
A. Background and model definition Experimental measurements of the c-axis transport in bulk crystals of BSCCO and other cuprates, summarized for example in Ref. [19], have been interpreted as evidence of interlayer tunneling dominated by disorder-mediated, incoherent processes. The c-axis superfluid stiffness, accessible through the measurement of the c-axis London penetration depth [24,25], provides a particularly clear evidence. Experimentally, the temperature dependence of the c-axis superfluid stiffness in clean single crystals was observed to follow an approximate power-law behavior ρ c = a−bT α with α 2 at low temperatures, whereas the in-plane stiffness showed a T -linear dependence [26]. The latter is the canonical behavior expected of a clean d-wave superconductor, reflecting the presence of lowenergy excitations with a Dirac spectrum [27]. Models with coherent tunneling between CuO 2 predict the same linear T -dependence for the c-axis stiffness [28], in clear disagreement with experimental data. If the interlayer tunneling were dominated by the oxygen-assisted processes (the g 1 term in Eq. (1)) then theory predicts ρ c = a − bT 5 [13], again at variance with experiment.
As demonstrated in Refs. [15][16][17]20] a description that captures the correct ∼ T 2 scaling along the c-axis (while preserving the T -linear behavior in the ab plane) can be given using the incoherent c-axis tunneling approach. In the following we shall review the relevant model and then apply it to the problem of a twisted bilayer.
A minimal model of the uncoupled bilayer system consists of the second-quantized Hamiltonian where l = 1, 2 denotes the layer index, Nambu-Gorkov spinors Ψ kl = (c kl↑ , c −kl↓ ) T . In the BCS approximation, we have with Pauli matrices σ j acting in the Nambu space and ∆ kl , ∆ kl denoting real and imaginary parts of the superconducting gap function. We adopt units such that = e = k B = m e = a 0 = 1, where mass is measured in units of electron mass m e and length scales in units of lattice constant a 0 . To make the model tractable we assume a simple parabolic band dispersion given by ξ k = k 2 /2m − µ in each layer. (In Sec. IV we consider a more realistic band structure and show that it leads to similar results.) The two superconducting d-wave order parameters are where ϕ is the phase difference and α k denotes the polar angle of k. The layers are coupled by the term and H = H 0 + H constitutes the full model. The lack of momentum conservation in Eq. (7) is the defining feature of the incoherent tunneling models and originates, physically, from the disorder present in the spacer layers separating the copper-oxygen planes. The disorder is captured via a set of Gaussian-distributed random variables γ q of average γ q = 0 and variance given by The scale Λ defines the characteristic momentum change that an electron undergoes when tunneling between the two layers. The factor 1/3 is chosen to reproduce the phase diagram of the coherent model of Ref. [7] in the limit Λ → 0 for the same value of g. For simplicity we have neglected any θ-dependence of the interlayer coupling although we expect the randomness to be stronger in twisted samples due to the increase in interface roughness, added strain, and moiré lattice modulations. The above form of incoherent interlayer tunneling is consistent with all lattice symmetries because γ q vanishes on average. This constitutes the key difference to a coherent coupling of the form k (g c † k,1 c k,2 + h.c.). As was pointed out in Ref. [12], at 45 • twist the two participating Cu d-orbitals transform under a 2D representation of D 4d and a coherent tunneling term is therefore not invariant under S 8 : (c k,1 , c k,2 ) → (c k,2 , −c k,1 ). It thus vanishes by virtue of the same argument as the Josephson coupling in Eq. (3). This is indeed owed to the coincidence of the atomic Cu orbitals transforming under the same representation as the superconducting order parameters.
It is instructive to consider the limit Λ → 0 of the incoherent tunneling in Eq. 8. Here, one has γ * q γ q+p = g 2 δ q,0 δ p,0 /3 and momentum is conserved in the interlayer tunneling process. Yet, Λ → 0 is not the clean limit in the sense that, in real space, it corresponds to the case where macroscopic regions are correlated with the same random value of interlayer tunneling g/ √ 3. From the viewpoint of disorder-induced incoherence, the random values of g should only be correlated in the vicinity of an impurity which sets the appropriate scale for 1/Λ. In the discussion below, the Λ → 0 thus serves as an abstract but convenient reference point that connects the present model to the calculations in the original work [7]. We will refer to it as the coherent limit.

B. Free energy and phase diagram
The physics of T -breaking is captured by the ϕdependence of the free energy. To determine the free energy we begin by implementing the global gauge transformation (c k1 , c k2 ) → c k1 e iϕ/4 , c k2 e −iϕ/4 which moves the superconducting phase difference from the order parameters in Eq. (6) to H according to In this gauge, the disorder-averaged interlayer current is given by where G(k, k , τ ) = T τ c k (τ )c † k (0) is the full imaginary time ordered Green's function of the disordered system and the current vertex is Note that the trace is to be performed over all momenta k, q and Matsubara frequencies ω n = (2n + 1)π/β, in addition to interlayer and Nambu indices. From the Josephson relation J(ϕ) = 2∂F(ϕ)/∂ϕ one then obtains the functional dependence of the free energy on the interlayer phase difference ϕ by simple integration. We expand Eq. (10) up to fourth order in g while treating H as a perturbation. Three different terms arise, which are diagrammatically represented in Fig. 1. Panel (a) corresponds to the term which is quadratic in g. Here, is the unperturbed, translationally invariant Green's function with H k = diag(H k1 , H k2 ) and is the impurity vertex. The disorder average acts on γ q factors and is performed according to Eq. (8). The diagrams in Fig. 1(b-c) represent terms of order g 4 : where impurity averaging is assumed but not explicitly shown for clarity of notation. Evaluating the traces, we obtain the current of the form with coefficients, to lowest order of g, Here, we have defined and ( * ) denotes a convolution integral, a k * b k = q a q b k−q . The quasiparticle dispersion of the unperturbed bands is given by E kl = ξ 2 k + ∆ 2 kl . From Eq. (14) one obtains the free energy The T -breaking phase transition occurs as a consequence of competition between cos ϕ and cos 2ϕ terms. Clearly, J c2 > 0 and the ground state acquires a finite phase difference for where it spontaneously breaks T . From our discussion in Sec. II it follows that that J c1 must vanish at twist of θ = π/4. Explicitly, one can see this result as follows. The functions E k,l , |γ k | 2 transform under the A 1g irrep of D 4h whereas the f k,1 , f k,2 transform under B 1g and B 2g , respectively. We note that convolution with the A 1g -symmetric impurity distribution |γ k | 2 does not change the symmetry of the convolution integral. Hence, j k transforms under B 1g ⊗ B 2g = A 2g and all terms in Eq. (15) average to zero at 45 • -twist. However, j 2 k is A 1g -symmetric and J c2 will consequently be finite and positive. Thus it is clear that condition (20) is generally satisfied at θ = 45 • and T will always be broken as soon as the system enters the SC state below T c . We conclude that impurity-mediated tunneling must not qualitatively change the T -breaking phase diagram relative to the model of Ref. [7]. The incoherent tunneling, however, shifts the phase boundaries. As shown in Appendix A, J c1 ∼ 1/Λ and J c2 ∼ 1/Λ 2 . Since J c1 is only weakly dependent on θ, and J c1 vanishes linearly around 45 • twist, it follows from Eq. (20) that the width of the T -breaking phase space is proportional to J c2 (θ = 0)/J c1 (θ = 0) ∼ 1/Λ. In the perfectly incoherent limit, Λ → ∞, the free energy becomes independent of ϕ and the T -breaking phase disappears.
To quantitatively ascertain the effect of incoherent tunneling on the phase diagram, we numerically evaluate the coefficients J ci . In principle, all Matsubara sums can be evaluated analytically, at the cost of removing the simple convolution structure in Eq. (17). This leaves three remaining momentum integrals to be numerically evaluated at complexity O(N 3 ) where N is the number of k-points of the 2D mesh used to perform the integrals. A more efficient approach is to exploit the convolution structure of Eq. (17) using the fast Fourier transform (fft) algorithm and numerically evaluate M Matsubara frequencies, affording evaluation of diagrams Fig. 1(a-b) at order O(M N log N ). The crossed diagram Fig. 1(c) does not possess a convolution structure. As we show in Appendix B, it can be evaluated at a cost of O(M N 2 ).
The resulting phase diagram is shown in Fig. 2 for coupling strength g = 10.5 meV and several values of Λ. We see that the T -breaking phase space is largest in the 'clean' limit Λ → 0 where it extends between (45 ± 6) • at T = 0. Increasing Λ gradually reduces the extent of the T -broken phase which eventually vanishes in the perfectly incoherent limit when Λ ∼ k F , i.e. when impu-rity correlations are on the scale of the lattice constant. Physically, this occurs because at this level of incoherence the Cooper pair essentially looses all memory of its momentum structure in the process of tunneling between layers.

C. Spectral gap and topological superconductivity
Having discerned the fate of the T -breaking phase in the presence of impurity-mediated tunneling we proceed to examine the topological properties of the resulting ground state. In the clean limit, T -breaking establishes a topological phase with Chern number C = 4 [7]. Since the disordered model is connected to the clean case by taking the limit Λ → 0, it is reasonable to expect the same C = 4 phase as long as the quasiparticle gap does not close.
Here, we show that these expectations are indeed met. To this end, we evaluate the Green's function in the Born approximation where with layer-indices τ = ±1. Here, we regularized the continuum model on a square lattice using + τ sin k x sin k y e −iϕ sin θ] with parameters chosen to match the continuum model. Following the method introduced in Ref. [29], we compute a spatially resolved Green function in the presence of a strong repulsive potential at x = 0 which simulates an edge and thus allows us to inspect the edge modes of the disordered system. We outline the method and give a derivation of Eq. (24) in Appendix C. In Fig. 3 we plot the analytically continued boundary spectral function at the edge (x = 1) as well as the bulk spectral function. We clearly observe two chiral edge modes traversing the bulk gap thus confirming the non-trivial topology of the system. The edge modes display a degeneracy in the layer degree of freedom, suggesting that the model is in a topological phase with Chern number C = 4. The bulk gap is reduced but remains finite as Λ increases.

IV. LATTICE MODEL
So far we have looked at the role of disorder in a continuum formulation of a twisted bilayer. The two-site unit cell of the regularized model allowed for analytical expressions for the layer Green's function to which we have systematically added incoherent interlayer tunneling and calculated the free energy up to fourth order in g. Another approach to tackle the problem and corroborate the results in a more general setting is to perform a BCS mean field theory on two twisted square lattices that represent the two CuO planes. In this case one can incorporate more realistic band structures but it is also harder to obtain the Green's functions analytically using Feynman diagrams. The reason is twofold: For one, at an arbitrary twist the lattice model is not commensurate. Secondly, at commensurate twist angles close to 45 • , the moiré unit cell contains many sites. To get around this, we perform a brute force disorder average wherein several disorder realizations with the same microscopic parameters are taken into account. While it limits us to real space, such a treatment is exact because all orders in perturbation theory are implicitly accounted for.
For each layer, we consider a square lattice Hubbard model with nearest neighbor density-density interactions such that a mean-field decoupling produces a d-wave order parameter. Including the interlayer tunneling processes with amplitudes g ij , the bilayer is described by where l is a layer index, t (t ) is the (next-)nearestneighbor hopping amplitude, µ is the chemical poten-tial that controls on-site particle density n iσl and ∆ ij,l denotes the complex order parameter on the bond connecting sites i and j on layer l. Considering a fully coherent interlayer tunneling, Ref. [7] employs a circularly symmetric, exponentially decaying form g ij = e −(rij −c)/ρ which connects sites i and j separated by r ij . Therein c in the interlayer separation and ρ is defined by the radial extent of the participating orbitals. The twist angle θ between the layers determines connectivity and the strength of the interlayer tunnelings. The free energy of this model shows a double-well structure for twist angles around 45 • [7].
To incorporate incoherent processes, we introduce a random tunneling factor that vanishes on average, but encodes the correlation between different processes depending on spatial separation. That is, where R = (r i +r j )/2 denotes the center of mass location of the hopping and g R = 0, Analogous to the parameter Λ in the continuum model, Λ sets the length scale for the correlation between different tunneling amplitudes and is indicative of disorder strength. We distinguish the two simply because of the slightly differing definitions. To simulate the Fermi surface of optimally doped BSCCO with hole pocket around (π, π), we set t = 153meV, t = −0.45t and µ = −1.35t [30]. Further, we choose c = 2.2 and ρ = 0.4 (in units of the lattice constant) to set interlayer distances. The dwave order parameters in cuprates originates in the CuO planes and the interlayer coupling is a minor perturbation that does not influence the order parameter magnitude. In other words, temperature dependence of the gap in each layer is independent of twist and coupling strength strength g, which we peg at 20meV. Therefore, we use a ∆ calculated self-consistently in a monolayer, which has a maximum of ∼ 40meV at 0K in accordance with experimental findings in cuprates [31,32].
To look for T -breaking we examine the free energy of the system, which can be calculated from the eigenvalues E i of the BdG Hamiltonian (26): In particular, for a given twist θ and disorder parameter Λ, we draw from the distribution (27) and average the free energy over 50 independent realizations. We choose a square bilayer sample as shown in Fig. 4(a), but the results are independent of the shape. Further, the exact number of sites in the system depends on the cut and the twist angle, but the free energy does not show an appreciable change beyond ∼ 900 sites per layer. In agreement < l a t e x i t s h a 1 _ b a s e 6 4 = " E r p m s 8 h w D i 2 k / 7 t S n T T g A U R Q V y s = " > A A A B / 3 i c b V D L S g N B E O z 1 G e M r 6 t H L Y B D i J e x K Q L 0 F v X i M a B 6 Q L G F 2 M p s M m Z 0 d Z m a F Z c n B D / C q n + B N v P o p f o G / 4 S T Z g y Y W N B R V 3 X R 3 B Z I z b V z 3 y 1 l Z X V v f 2 C x s F b d 3 d v f 2 S w e H L R 0 n i t A m i X m s O g H W l D N B m 4 Y Z T j t S U R w F n L a D 8 c 3 U b z 9 S p V k s H k w q q R / h o W A h I 9 h Y 6 b 6 C z / q l s l t 1 Z 0 D L x M t J G X I 0 + q X v 3 i A m S U S F I R x r 3 f V c a f w M K 8 M I p 5 N i L 9 F U Y j L G Q 9 q 1 V O C I a j + b n T p B p 1 Y Z o D B W t o R B M / X 3 R I Y j r d M o s J 0 R N i O 9 6 E 3 F f z 0 5 S j U j e m G 9 C S / 9 j A m Z G C r I f H u Y c G R i N A 0 D D Z i i x P D U E k w U s w 8 g M s I K E 2 M j K 9 p k v M U c l k n r v O r V q l d 3 t X L 9 O s + o A M d w A h X w 4 A L q c A s N a A K B I T z D C 7 w 6 T 8 6 b 8 + 5 8 z F t X n H z m C P 7 A + f w B k t G W T g = = < / l a t e x i t > (b) < l a t e x i t s h a 1 _ b a s e 6 4 = " n k r J E 2 T H g O K 9 2 s U N 8 D n x l b v V 2 3 w = " > A A A B / 3 i c b V D L S g N B E O z 1 G e M r 6 t H L Y B D i J e x K Q L 0 F v X i M a B 6 Q L G F 2 M p s M m Z 0 d Z m a F Z c n B D / C q n + B N v P o p f o G / 4 S T Z g y Y W N B R V 3 X R 3 B Z I z b V z 3 y 1 l Z X V v f 2 C x s F b d 3 d v f 2 S w e H L R 0 n i t A m i X m s O g H W l D N B m 4 Y Z T j t S U R w F n L a D 8 c 3 U b z 9 S p V k s H k w q q R / h o W A h I 9 h Y 6 b 4 S n P V L Z b f q z o C W i Z e T M u R o 9 E v f v U F M k o g K Q z j W u u u 5 0 v g Z V o Y R T i f F X q K p x G S M h 7 R r q c A R 1 X 4 2 O 3 W C T q 0 y Q G G s b A m D Z u r v i Q x H W q d R Y D s j b E Z 6 0 Z u K / 3 p y l G p G 9 M J 6 E 1 7 6 G R M y M V S Q + f Y w 4 c j E a B o G G j B F i e G p J Z g o Z h 9 A Z I Q V J s Z G V r T J e I s 5 L J P W e d W r V a / u a u X 6 d Z 5 R A Y 7 h B C r g w Q X U 4 R Y a 0 A Q C Q 3 i G F 3 h 1 n p w 3 5 9 3 5 m L e u O P n M E f y B 8 / k D l G u W T w = = < / l a t e x i t > (c) < l a t e x i t s h a 1 _ b a s e 6 4 = " 6 s T x j + c Q 3 s P X Z d g D T I J y n u 1 i 9 5 o = " > A A A B / 3 i c b V D L S g N B E O z 1 G e M r 6 t H L Y B D i J e x K Q L 0 F v X i M a B 6 Q L G F 2 M p s M m Z 0 d Z m a F Z c n B D / C q n + B N v P o p f o G / 4 S T Z g y Y W N B R V 3 X R 3 B Z I z b V z 3 y 1 l Z X V v f 2 C x s F b d 3 d v f 2 S w e H L R 0 n i t A m i X m s O g H W l D N B m 4 Y Z T j t S U R w F n L a D 8 c 3 U b z 9 S p V k s H k w q q R / h o W A h I 9 h Y 6 b 5 C z v q l s l t 1 Z 0 D L x M t J G X I 0 + q X v 3 i A m S U S F I R x r 3 f V c a f w M K 8 M I p 5 N i L 9 F U Y j L G Q 9 q 1 V O C I a j + b n T p B p 1 Y Z o D B W t o R B M / X 3 R I Y j r d M o s J 0 R N i O 9 6 E 3 F f z 0 5 S j U j e m G 9 C S / 9 j A m Z G C r I f H u Y c G R i N A 0 D D Z i i x P D U E k w U s w 8 g M s I K E 2 M j K 9 p k v M U c l k n r v O r V q l d 3 t X L 9 O s + o A M d w A h X w 4 A L q c A s N a A K B I T z D C 7 w 6 T 8 6 b 8 + 5 8 z F t X n H z m C P 7 A + f w B l g W W U A = = < / l a t e x i t > with the continuum model, Fig. 4(b) shows that the presence of T -breaking free energy minima is controlled byΛ. Namely, small values ofΛ support the T -broken ground state while larger values do not.

V. CONCLUSIONS
Twisted bilayers of high-T c cuprates hold the potential for realizing topological superconductivity, wherein a topological gap is spontaneously induced. As per the symmetry informed momentum space form factors, which determine the electron hopping between interlayer Cu atoms, the tunneling amplitude vanishes along the nodal directions and a spectral gap may not appear. In this work we highlight that an important aspect to consider in such an analysis is the disorder mediated tunneling. Not only does disorder appear naturally due to oxygen doping and interfacial defects, but incorporating momentum non-conservation has been shown to better represent experimental data in clean single crystals.
Using perturbative diagrammatic calculations and disorder averaging on the lattice, we find that an experimentally motivated incoherent tunneling model that respects all point group symmetries of the physical system gives rise to a qualitatively similar phase diagram as obtained in Ref. [7]. Specifically we find a substantial range of twist angles around 45 • and temperatures where spontaneous T -breaking occurs and produces a fully gapped topological phase with non-zero Chern number. The angular extent of the T -broken phase depends on the disorder length scale Λ −1 where the coherent limit Λ → 0 recovers the phase diagram of Ref. [7] and increasing Λ corresponds to a shrinking extent of the topological phase. Only when the incoherence length scale is comparable to the Fermi momentum, the twist angle for spontaneous T breaking is reduced to exactly 45 • .
From an experimental point of view, the inhomogeneity due to oxygen doping the BiO planes of untwisted BSCCO was found to be correlated over ≈ 14 A [33]. Since the CuO plane lattice constant is ≈ 5 A, that amounts to a correlations over 3 unit cells, i.e.,Λ ≈ 0.3. In a twisted geometry, one may expect the characteristic length scale to decrease and, hence, the estimate forΛ could shift up. That said, the role of complex atomic arrangements, moiré length scales and strong correlations are difficult to incorporate into such a heuristic reasoning. One would probably have to await data from complementary experimental probes, such as transport and optical response, to discern the nature of the superconducting state around 45 • .
It was noted in Ref. [12] that the measured critical current density J c in both twisted and untwisted Bi2212 is about factor of 500 smaller than the theory prediction based on the slave-boson mean field theory of a t-J model used in that study. We checked that a similar discrepancy occurs in the calculation using BCS mean field theory of Ref. [7]. As indicated in Fig. 5 the discrepancy is somewhat reduced in the incoherent tunneling model (by about one order of magnitude at large Λ) but nevertheless significant disagreement with experiment persists. As noted in Ref. [20] this is a known problem that affects superconductors in the cuprate family and becomes increasingly severe in the underdoped part of their phase diagram. A phenomenological fix can be implemented [20] by restricting the momentum sums in the expression for J c to patches of linear size ∼ x (the hole doping) around the nodal points of the d-wave order parameter. This modification leaves the temperature dependence of ρ ab (T ) and ρ c (T ) unchanged, but reduces their T = 0 magnitude to experimentally observed values. It similarly fixes the problem with J c . As with many aspects of cuprates a truly microscopic understanding of this phenomenon remains a challenge to the theory community. With regards to twisted cuprate bilayers it would be interesting to explore the effect of the phenomenological fix outlined above on the phase diagram but we leave this to future work. be reduced to O(M N 2 ) by exploiting the Gaussian form of |γ q | 2 ∼ e −q 2 /Λ 2 . For ease of notation, we will omit prefactors and Matsubara indices n and set Λ = 1.
We perform the substitution q → q − k and q → q − k which modifies Eq. (B1) to qq e −(q−q ) 2 /2 f q,2 f q ,2 (B2) Choosing the center of mass frame Q = q + q , P = q − q , the term is expressed as Them sum over k requires N 2 steps (N steps to perform the k-summation times N steps to resolve the residual Qdependence). Analogously, the Q-sum can be computed in N 2 steps and the remaining P-sum is performed in N steps. Taking into account the Matsubara summation, we arrive at a total complexity of O(M N 2 ).  FIG. 9. Layer-resolved bulk (a-d) and boundary (e-h) spectra for incoherently coupled cuprate bilayers with Λ = 0 (left two columns) and Λ/kF = 0.08 (right two columns) at 45 • twist angle for superconducting gaps plotted in Fig. 8(c-d). In each layer, the spectral functions display two distinct chiral edge modes, indicating a total Chern number C = 4.
The two gap functions are plotted in the bottom panels of Fig. 8. Here, the two gaps are no longer related to each other through a mirror symmetry. The corresponding layer-resolved spectral functions hence differ in the two layers of the bilayer structure, as seen in Fig. 9. In the absence of any interlayer coupling, g = 0, each layer possess four Dirac cones. Finite interlayer coupling g has two effects: (a) it gaps the Dirac cones and (b) induces the four gapped cones of layer 1 onto layer 2 and vice versa. The induced gapped Dirac cones are characterized by light spectral weight in Fig. 9. The positions of gapped Dirac cones correspond to intersections of the Fermi surface in Fig. 8(c-d) with the gap nodes, projected onto the k y axis. In each layer, the gapped Dirac cones are traversed by two Chiral edge modes, indicating a total topological Chern invariant of C = 4. This is still the case in the presence of impurities, Λ/k F = 0.08. Here, the incoherent nature of the interlayer tunneling causes a broadening of the induced gapped Dirac cones.