Science, Faculty of
Physics and Astronomy, Department of
DSpace
UBCV
Tummuru, Tarun Reddy
Doctor of Philosophy - PhD
University of British Columbia
The notion of twisting and stacking two-dimensional van der Waals materials has emerged as a paradigm for realizing novel electronic states. With the goal of engineering topological superconductivity, we go beyond the archetypal example of twisted bilayer graphene and consider structures composed of proximitized quantum wires and high-Tc cuprate superconductors. On the basis of symmetry arguments, backed by microscopic modeling, these setups are predicted to spontaneously break time-reversal symmetry and give rise to topological p- and d-wave phases with chiral Majorana edge modes. Given the experimental pertinence, we investigate a suite of Josephson transport phenomena in twisted cuprate bilayers and identify signatures of the topological phase. Further, the analysis is extended to twisted multilayers where the occurrence of higher-order band crossings makes the groundstate susceptible to secondary instabilities and topological superconductivity appears for an extended range of twists.
Topological superconductivityin twisted 2D structuresbyTarun Reddy TummuruB.Sc., Shiv Nadar University, 2017M.Sc., The University of British Columbia, 2019A THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFDOCTOR OF PHILOSOPHYinThe Faculty of Graduate and Postdoctoral Studies(Physics)THE UNIVERSITY OF BRITISH COLUMBIA(Vancouver)February 2023© Tarun Reddy Tummuru 2023The following individuals certify that they have read, and recommend to the Faculty ofGraduate and Postdoctoral Studies for acceptance, the thesis entitled:Topological superconductivity in twisted 2D structuressubmitted by Tarun Reddy Tummuru in partial fulfillment of the requirements for thedegree of Doctor of Philosophy in Physics.Examining Committee:Marcel Franz, Professor, Physics and Astronomy, UBCSupervisorDouglas Bonn, Professor, Physics and Astronomy, UBCUniversity ExaminerRoman Krems, Professor, Chemistry, UBCUniversity ExaminerMoshe Rozali, Professor, Physics and Astronomy, UBCSupervisory Committee MemberAdditional Supervisory Committee:Ziliang Ye, Assistant Professor, Physics and Astronomy, UBCSupervisory Committee MemberAndrew Potter, Assistant Professor, Physics and Astronomy, UBCSupervisory Committee MemberiiAbstractThe notion of twisting and stacking two-dimensional van der Waals materials has emergedas a paradigm for realizing novel electronic states. With the goal of engineering topologicalsuperconductivity, we go beyond the archetypal example of twisted bilayer grapheneand consider structures composed of proximitized quantum wires and high-Tc cupratesuperconductors. On the basis of symmetry arguments, backed by microscopic modeling,these setups are predicted to spontaneously break time-reversal symmetry and give rise totopological p- and d-wave phases with chiral Majorana edge modes. Given the experimentalpertinence, we investigate a suite of Josephson transport phenomena in twisted cupratebilayers and identify signatures of the topological phase. Further, the analysis is extendedto twisted multilayers where the occurrence of higher-order band crossings makes thegroundstate susceptible to secondary instabilities and topological superconductivity appearsfor an extended range of twists.iiiLay SummaryA superconductor is a phase of matter that allows electric current to flow through withoutresistance and expels magnetic fields. Since their discovery in 1911, several materials,including some commonplace metals like aluminum and lead, have been shown to becomesuperconducting when cooled to very low temperatures. Over the past two decades, a specialclass of superconductors called topological superconductors have garnered academic interestbecause of their relevance for quantum computation. However, material characterizationhas been tricky and experimental consensus is still lacking. Building on recent experimentalinnovations, we propose conceptually new ways to engineer topological states: assemblingindividual constituents in a specific manner, like in a Lego.ivPrefaceContents of this thesis have been adapted from published works as follows:• Chapter 2 : T Tummuru, O Can and M Franz, ‘Chiral p-wave superconductivity in atwisted array of proximitized quantum wires’ Phys Rev B 103, L100501 (2021).While we were looking at twisted cuprates, Marcel brought up the possibility ofengineering a chiral p-wave state. I set up the model and determined the phasediagram, while Marcel simulated the vortex and Oguz looked at an alternative modelfor the quantum wires that takes into account the spin-orbit coupling and magneticfields. I drafted the manuscript in collaboration with Marcel and Oguz.• Chapter 3 : O Can, T Tummuru, R Day, I Elfimov, A Damascelli and M Franz, ‘High-temperature topological superconductivity in twisted double-layer copper oxides’ NatPhys 17, 519 (2021).The idea to twist cuprates originated from Marcel, who also formulated the Ginzburg-Landau theory and the continuum model. Together with Oguz, we developed thelattice model and performed mean-field analysis. Oguz carried out calculations ofthe phase diagram, while I simulated strip geometries and Josephson effects. RyanDay and Ilya Elfimov performed DFT calculations. The manuscript was compiled byMarcel with input from all authors.• Chapter 4 : T Tummuru, S Plugge and M Franz, ‘Josephson effects in twisted cupratebilayers’ Phys Rev B 105, 064501 (2022).Our interest in the Josephson physics of twisted cuprates was sparked by experimentaldata from Philip Kim’s group at Harvard University. Marcel analyzed the temperaturedependence of current and the gap in thicker flakes using continuum models. TheFraunhofer pattern calculation was performed by Stepan Plugge. In addition tothe Shapiro steps, I helped determine the temperature and interlayer tunnelingdependence of the critical current using a lattice formulation of the problem. Stepanand I worked on the manuscript, originally drafted by Marcel.• Chapter 5 : T Tummuru, E´ Lantagne-Hurtubise and M Franz, ‘Twisted multilayernodal superconductors’ Phys Rev B 106, 014520 (2022).The possibility of cubic band touchings in trilayers was pointed out by Marcel.Working with E´tienne Lantagne-Hurtubise, we extended that idea and looked atphase diagrams of multilayers in different twist configurations and performed self-consistent mean-field calculations. We drafted the work with input from Marcel.vTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiLay Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiiList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ixAcknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . x1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 Chiral p-wave superconductivity with proximitized quantum wires . . 42.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42.2 Ginzburg-Landau theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42.3 Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52.4 Majorana zero mode in a vortex . . . . . . . . . . . . . . . . . . . . . . . . 82.5 Outlook . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93 Chiral d-wave superconductivity in twisted bilayer cuprates . . . . . . . 113.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113.2 Ginzburg-Landau theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123.3 Microscopic models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 143.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 184 Josephson effects in twisted cuprate bilayers . . . . . . . . . . . . . . . . 194.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 194.2 Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 224.3 Temperature dependence of critical current . . . . . . . . . . . . . . . . . . 244.4 Estimate of interlayer coupling . . . . . . . . . . . . . . . . . . . . . . . . . 264.5 Interlayer coupling dependence of critical current . . . . . . . . . . . . . . . 274.6 Magic angles and secondary order parameters . . . . . . . . . . . . . . . . 284.7 Fraunhofer pattern . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 294.8 Shapiro steps . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31viTable of Contents4.9 Thick flakes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 334.10 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 345 Twisted multilayer nodal superconductors . . . . . . . . . . . . . . . . . . 375.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 375.2 Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 385.3 Trilayer: Alternating twist . . . . . . . . . . . . . . . . . . . . . . . . . . . 395.3.1 Spontaneous T -breaking . . . . . . . . . . . . . . . . . . . . . . . . 415.3.2 Quadratic band touching . . . . . . . . . . . . . . . . . . . . . . . . 415.4 Trilayer: Chiral twist . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 455.4.1 Spontaneous T -breaking . . . . . . . . . . . . . . . . . . . . . . . . 455.4.2 Cubic band touching . . . . . . . . . . . . . . . . . . . . . . . . . . 455.5 Generalization to multilayers . . . . . . . . . . . . . . . . . . . . . . . . . . 475.5.1 Alternating twist . . . . . . . . . . . . . . . . . . . . . . . . . . . . 475.5.2 Chiral twist . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 485.6 Outlook . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 496 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55AppendicesA Appendices for Chapter 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64A.1 Chern number . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64A.2 Construction from nanowires . . . . . . . . . . . . . . . . . . . . . . . . . . 64A.3 Details of the vortex calculation . . . . . . . . . . . . . . . . . . . . . . . . 67B Appendices for Chapter 3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68B.1 Microscopic derivation of GL coefficients . . . . . . . . . . . . . . . . . . . 68B.2 Self-consistent lattice mean-field theory . . . . . . . . . . . . . . . . . . . . 69C Appendices for Chapter 4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72C.1 Interlayer supercurrent analysis . . . . . . . . . . . . . . . . . . . . . . . . 72C.2 Expansion in powers of g . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73C.3 In-plane superfluid stiffness . . . . . . . . . . . . . . . . . . . . . . . . . . . 73D Appendices for Chapter 5 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76D.1 Layer decoupling unitary transformation . . . . . . . . . . . . . . . . . . . 76D.2 Low energy theory at magic angles . . . . . . . . . . . . . . . . . . . . . . . 77D.3 Secondary instabilities . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79D.4 Density of states and gap scaling . . . . . . . . . . . . . . . . . . . . . . . . 79D.5 Lattice model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82viiList of Tables1.1 Periodic table of topological insulators and superconductors . . . . . . . . . 25.1 Low energy spectra at magic angles . . . . . . . . . . . . . . . . . . . . . . . 48viiiList of Figures2.1 Geometry of the proposed device . . . . . . . . . . . . . . . . . . . . . . . . 52.2 Kitaev wire array model and its phases . . . . . . . . . . . . . . . . . . . . . 62.3 Majorana zero mode in a vortex . . . . . . . . . . . . . . . . . . . . . . . . 83.1 Time-reversal symmetry breaking in twisted bilayers . . . . . . . . . . . . . 133.2 Lattice model results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 154.1 Schematic of the setup and phase diagram . . . . . . . . . . . . . . . . . . . 204.2 Temperature dependence of supercurrent . . . . . . . . . . . . . . . . . . . . 234.3 Critical current and superfluid stiffness . . . . . . . . . . . . . . . . . . . . . 254.4 Lattice model: critical current and g dependence . . . . . . . . . . . . . . . 274.5 Fraunhofer patterns . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 304.6 Shapiro steps . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 314.7 Critical currents in thick flakes . . . . . . . . . . . . . . . . . . . . . . . . . 334.8 Spectral gap in thick flakes . . . . . . . . . . . . . . . . . . . . . . . . . . . 355.1 Time-reversal symmetry breaking in twisted trilayers . . . . . . . . . . . . . 405.2 Chiral edge modes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 425.3 Magic angles and density of states . . . . . . . . . . . . . . . . . . . . . . . 435.4 Phase diagrams of twisted trilayers . . . . . . . . . . . . . . . . . . . . . . . 445.5 Mean-field results on quadrilayers . . . . . . . . . . . . . . . . . . . . . . . . 465.6 Phase diagram comparison of chiral twist multilayers . . . . . . . . . . . . . 475.7 Low energy dispersion of chiral quadrilayer . . . . . . . . . . . . . . . . . . 49A.1 Energy spectra with SOC dispersion along the wires . . . . . . . . . . . . . 65A.2 Energy spectra with generic SOC dispersion . . . . . . . . . . . . . . . . . . 66B.1 Analysis of energy integrands . . . . . . . . . . . . . . . . . . . . . . . . . . 69B.2 Commensurate twist angles . . . . . . . . . . . . . . . . . . . . . . . . . . . 70C.1 Superconducting gap and in-plane superfluid stiffness . . . . . . . . . . . . . 74D.1 Secondary s-wave order in trilayers . . . . . . . . . . . . . . . . . . . . . . . 80ixAcknowledgementsLooking back, the time as a graduate student has been one of the most enriching anddynamic periods of my life. Numerous individuals have molded my experiences and I wouldlike to acknowledge them all, including those I inadvertently missed below.Foremost, I would like to thank my advisor Marcel Franz whose guidance shaped theform and substance of the work described here. Marcel’s ability to distill complex problemsinto succinct questions rooted in physical intuition has often provided me with clarityof thought and a sense of direction. His constant encouragement to pursue interestingquestions has been crucial in learning to think as a researcher. Marcel’s structured courseswith an overarching view of the field of topology and other topics in condensed matterwere music to the uninitiated ears. The group outings and get-togethers hosted by Noraand Marcel will always remain pleasant memories.Research has been fun in good part due to conversations and collaborations with peersat the department. Thanks to: Oguz Can for exploring together the exciting physics oftwisted cuprates and initiating discussions on myriad physics concepts; Ryan Day forexplaining details of DFT simulations; Stephan Plugge for sharing his expertise withtransport calculations; E´tienne Lantagne-Hurtubise for identifying important aspects inthe multilayer story and innumerable conversations about physics and beyond; RafaelHaenel who helped me grasp the beauty of group theory (among other topics) and hasbeen a companion during several late nights spent on assignments and exams; ChengshuLi, who was ever willing to discuss and share his knowledge of numerics. Many othershave enlivened the meetings and everyday work: Sayak Dasgupta, Shannon Egan, StepanFomichev, Pedro Lopes, Alberto Nocera, Vedangi Pathak, Sharmistha Sahoo, Kyle Wamer,Wang Yang, Xiao-Xiao Zhang and Benjamin Zhou.Having come to UBC as a master’s student, I owe a lot to Ian Affleck for his supportand guidance during the initial years. Ian’s dedication to the craft, combined with his workethic and humility, has made a lasting impact. I take the opportunity to express gratitudeto Mona Berciu whose exceptional solid-state physics course formed the foundation formany graduate students. Besides, over multiple conversations, Mona freely shared herexperiences and provided invaluable advice regarding different aspects of life within andwithout academia. I would also like to thank the members of my supervisory committeeAndrew Potter, Moshe Rozali and Ziliang Ye for their valuable feedback.Over the last year, I’ve had the privilege of spending time at the University of Zurich.I am indebted to Titus Neupert for generously hosting and supporting me. WitnessingTitus’ grasp on the subject, the constant stream of ideas and the breadth of his knowledgehas been inspiring. I have also benefited from the mentorship of Toma´sˇ Bzdusˇek, whointroduced me to a plethora of mathematical concepts. The talks listened to, visitors metand schools attended have helped broaden my scope of thinking. I am also grateful forxAcknowledgementsconversations on an eclectic mix of topics with Nikita Astrakhantsev, Michael Denner,Mark Fischer, Bastien Lapierre, Patrick Lenggenhager, Aleksandra Nelson, AnastasiiaSkurativska, Martina Soldini, Glenn Wagner and Songbo Zhang. You all have beenextremely welcoming and filled my stay with warmth.Life outside of work has been adventurous and ever-interesting thanks to E´tienne,Fiana, Irene, Ji Eun, Jihyun, Khalid, Kyle, Kyoko, Oliver, Rafael, Roger, Samir, Samuel,Yuki, Zoe and several others at St. John’s College. If I could call Vancouver a home, it isbecause of you all.Thank you Hyunju for being an anchor through all times, thick and thin. Three yearsis a long time to be apart, but the conversations, adventures and laughs have sustainedme. Looking forward to the future!The journey over the past few years was possible because of the trust my family placedin me, even though these undertakings lead me to an uncertain path on the other side ofthe world. Thank you for the love and support.xiChapter 1IntroductionTopology studies the properties of spaces that are preserved under continuous deformations.Two geometric objects are said to be equivalent if they can be deformed into each otherwithout any tearing or gluing of parts. Mathematically, the letter D is homeomorphicto O, but not B. In a similar vein, condensed matter systems can be classified intoequivalence classes based on their Hamiltonians. Topological materials constitute a uniquestate of matter where topology of quantum states residing in the bulk manifests itself asunconventional properties on the boundary.The field has its origins in the experimental discovery of quantum Hall effect, wherequantized Hall conductance was measured in a 2D electron gas subjected to a magneticfield [1]. This phenomenon, which was eventually traced to the existence of a topologicalinvariant [2], inspired Haldane’s seminal work: quantized response could be elicited evenin the absence of an external magnetic field [3]. Christened as quantum anomalous Hallinsulator (for the lack of magnetic field), when time reversal symmetry T is broken, themodel develops a topological character that is revealed in conducting quantum stateslocalized to the edges of an otherwise insulating system. The topology is quantified interms of the so-called Chern number, which is necessarily an integer – with zero signalinga trivial topology. By combining two time-reversed copies of Haldane’s model, Kane andMele [4] provided the proof of principle for a T invariant topological state. Also knownas the quantum spin Hall insulator, the phase is labelled by a Z2 index and exhibitscounter-propagating edge states.Without the need for T breaking, it thus became apparent that non-trivial topologymight be realized in actual materials. Following extensive theoretical predictions, severalindicators of topology in normal sate band structures have been experimentally discovered.To name a few, signatures of the quantum spin Hall effect were observed in HgTe [5] andInAs/GaSb quantum wells [6], Chern insulating states have been reported in magneticallydoped topological insulator thin films [7] and, more recently, in graphene-based moire´materials [8, 9].Topological materials can be framed in an elegant mathematical framework thatencompasses different topological invariants and material realizations cited above [10, 11].The classification can be understood in terms of symmetries of the Bloch Hamiltonian hkin d spatial dimensions. The symmetries in consideration, time reversal T , particle-hole Pand chiral C, are very generic and their action can be expressed asT †hkT = h−k, P†hkP = −h−k and C†hkC = −hk. (1.1)The momentum sign change in the first two relations is a consequence of the anti-unitarynature of the symmetries. An anti-unitary operator can be written as a combination of aunitary operator and complex conjugation K. T and P come in two flavours, depending1Chapter 1. Introductionclass T P C d = 0 1 2 3A 0 0 0 Z 0 Z 0AIII 0 0 1 0 Z 0 ZAI + 0 0 Z 0 0 0BDI + + 1 Z2 Z 0 0D 0 + 0 Z2 Z2 ZDIII − + 1 0 Z2 Z2 ZAII − 0 0 2Z 0 Z2 Z2CII − − 1 0 2Z 0 Z2C 0 − 0 0 0 2Z 0CI + − 1 0 0 0 2ZTable 1.1: Periodic table of topological insulators and superconductors. The ten symmetryclasses in the first column are characterized by the presence (absence) of time-reversal T ,particle-hole P, and chiral C symmetry, as denoted by ±1 (0). d is the spatial dimensionand the entries Z, Z2, 2Z, represent possible values the topological invariant can take,while 0 implies that the system does not have a topological phase.on whether they square to identity or negative identity. Specifically, in the case of super-conductors, the Bogoliubov de Gennes (BdG) framework imposes particle-hole symmetrythat can be traced back to the Fermi statistics of the electronic degrees of freedom. Chiralsymmetry is defined as the combination of the other two: C = T P. Therefore, whenevertwo of them are present, the third is a symmetry as well. However, if both T and P areabsent, then C may or may not be a symmetry. Together, these considerations give us tensymmetry classes that are summarized along with the corresponding topological invariantsin Table 1.1.When the symmetries of a model are known, its groundstate may be labelled astopologically trivial or non-trivial according to the corresponding topological invariant,which is either an integer (Chern/winding number) or a Z2 quantity. The invariant signifiesthat the system remains in the same phase even when the Hamiltonian is continuouslydeformed provided that its energy eigenvalues do not become degenerate and the symmetriesare preserved. Importantly, each non-trivial entry in the table is expected, via the bulk-boundary correspondence, to have gapless states on the boundary.Among the different kinds of topological states, emergent quasiparticles in topologicalsuperconductors hold immense potential for technological advancement. Originally pro-posed as a solution of the Dirac equation, a Majorana fermion has the unusual property ofbeing its own anti-particle. In the condensed matter setting, Majorana zero modes (MZMs)are 0-dimensional charge neutral bound states that exist strictly at zero energy and obeynon-Abelian statistics [12, 13]. MZMs are ideal candidates for braiding based topologicalquantum computation that minimizes the influence of error inducing environmental per-turbations [14, 15]. Prominent examples of candidate topological materials with effectiveodd-parity Cooper pairing that are believed to host zero modes in vortex defects includeUPt3 [16], CuxBi2Se3 [17, 18] and FeTe0.55Se0.45 [19, 20]. CuxBi2Se3 is an interesting casebecause the normal state corresponds to a doped 3D topological insulator (Bi2Se3). In an-2Chapter 1. Introductionother manifestation, Majorana fermions appear as propagating states localized to the edgeof a chiral superconductor [21, 22]. Schemes to realize chiral spin-singlet superconductivityinclude c-axis twist Josephson junctions [23] and doped graphene [24].Despite the possibilities, proposals and interesting leads, conclusive evidence for topo-logical superconductivity and Majorana modes remains elusive. The dearth of naturallyoccurring materials prompted the notion of hybrid structures where a topological groundstate appears by design. A well-known example is a semiconductor nanowire proximitizedwith a regular s-wave superconductor that behaves like a spinless 1D p-wave superconductor[25–27]. Another platform that effectively simulates a chiral p-wave superconductor isthe Fu-Kane setup [28], where a 3D topological insulator is interfaced with a conven-tional superconductor. The existence of an MZM could be demonstrated spectroscopicallyas a quantized conductance peak, however, experimental consensus is still awaited [29].Similarly, the observation of 1D propagating Majorana modes at the edges of 2D chiralsuperconductors has been contested [30]. Given this context, we explore other ways toengineer 2D topological superconductivity with localized Majorana zero modes and/orpropagating Majorana edge modes. In particular, we focus on an experimental innovationthat has risen to prominence in the past few years: twisting and stacking 2D materials.Twistronics, from twist + electronics, produces emergent electronic states that areoriginally absent in the building blocks. The archetypal example, which was also the firstto be experimentally realized, is twisted bilayer graphene. A single layer of graphene isa semi-metal with gapless Dirac cones characterizing its dispersion. However, when twographene layers are twisted to 1.1◦, the spectrum is characterized by flat bands [31, 32].And because the kinetic energy of the electrons (given by the curvature of the bands)is quenched, interaction effects become prominent. It has now been established thattwisting can pave the way to an understanding of strong correlations and bring abouta slew of interesting insulating, topological and superconducting phases. Motivated bythese discoveries, twist angle engineering is now being applied to a range of van der Waalsmaterials [33] including more elaborate graphene-based structures [34–40], transition metaldichalcogenides [41–44], 2D magnets [45, 46] and thin films of high-Tc superconductors[47–51].In what follows, we propose an expansion of the programme of twistronics to nodalsuperconductors and show that it is indeed possible to create topological superconductors.We discuss potential routes to chiral p- and d-wave states using proximitized quantumnanowires and high-Tc cuprate superconductors, respectively. A monolayer d-wave su-perconductor belongs to class CI, but by combining two such layers with a twist onespontaneously breaks T and ends up in class C that is characterized by a 2Z invariant.A spinless 2D chiral p-wave superconductor, on the other hand, falls in class D and isdefined by a Z invariant. Here vortex defects are shown to host Majorana zero modes.Since twisting cuprates is within the reach of current fabrication techniques, we identifypertinent questions from the perspective of quantum transport and highlight signaturesof the topological state. Finally, as an extension of the bilayer setup, we model twistedmultilayers of cuprates which, in addition to chiral superconductivity, show interestingfeatures such as higher-order band crossings and extended topological phases.3Chapter 2Chiral p-wave superconductivitywith proximitized quantum wires2.1 IntroductionA key property of Majorana modes in condensed matter systems, which makes themmuch sought after, is their exchange statistics. Non-trivial statistics form the basis forfault-tolerant topological quantum computation [15, 52]. Moreover, a Majorana zero mode(MZM) is a quantum state that constitutes one-half of a usual complex fermion. By theway of enabling spatial separation of information encoded in a qubit, MZMs can serve asquantum memories that are immune to decoherence [14].Prominent examples of systems believed to host MZMs include proximitized semicon-ductor quantum wires [25–27, 53–56] that realize the Kitaev chain model [14] and theso-called Fu-Kane superconductor [28], which emerges at the interface of a 3D topologicalinsulator and a conventional superconductor. The latter may be realized intrinsically inFeTexSe1−x [19], with some experiments showing evidence for MZMs in vortex cores [57–60].A comprehensive discussion of these platforms can be found in recent review articles [61–66].Another setup that exhibits MZMs is a spinless px + ipy superconductor [12, 67] where,similar to the Fu-Kane model, zero modes are localized in vortex cores. However, px + ipysuperconductors have been notoriously difficult to find. The leading candidate, Sr2RuO4[68], was recently shown to be incompatible with the required spin-triplet pairing [69].In this chapter, we show that chiral px + ipy superconductivity can be actualized usingproximitized semiconductor nanowires. Quantum nanowires have been extensively studiedover the past decade and are known to behave as 1D p-wave superconductors [66]. Ouridea relies on the observation that a 2D superconductor with p-wave symmetry can bebuilt by assembling an array of quantum wires. Remarkably, when two such layers arestacked on top of each other with a relative twist angle of 90◦ (Fig. 2.1), a non-zero phasedifference develops between the superconducting order parameters, spontaneously breakinga residual time reversal symmetry T ′. As a result, each monolayer effectively behaves as afully gapped chiral superconductor. In what follows, we flesh out the idea and demonstrateby model calculations that such synthetic px+ ipy superconductors exhibit protected chiralMajorana edge modes and harbor MZMs in the cores of Abrikosov vortices.2.2 Ginzburg-Landau theoryTo illustrate the basic mechanism, let us analyze the properties of a system composed oftwo p-wave superconductors, each in the form of a monolayer, using a phenomenological42.3. ModelBSC coatingquantumwireFigure 2.1: Geometry of the proposed device. Semiconductor quantum wires coated onone side with a thin layer of conventional superconductor, such as Al, are assembled intotwo arrays and stacked with twist angle close to 90◦ as shown.Ginzburg-Landau (GL) theory. A time-reversal symmetry breaking state requires twocomplex order parameters, η1 and η2, one for each layer. Specifically, we consider thecase where the two order parameters have symmetries η1 ∼ px and η2 ∼ py. Cooperpair tunneling amplitude between two superconductors is given by an overlap of the twowavefunctions. Because px and py orbitals are mutually orthogonal, single Cooper pairscannot tunnel between the two layers and the lowest order coupling ∼ (η1η∗2 + c.c.) isprohibited. Taking this into account, the GL free energy functional, which is a scalarunder the symmetry operations of the system, is given by F = F1 + F2 + F12 withFa = α|ηa|2 + β|ηa|4 andF12 = γ|η1|2|η2|2 + ϵ(η21η∗22 + c.c.). (2.1)Here α, β, γ and ϵ are real phenomenological parameters that encode physical informationof microscopic origin and a = 1, 2 labels the layers. We focus on the case of spatiallyuniform superconducting (SC) order and only include gradient free terms. Without loss ofgenerality, the two layers can be assumed to be identical, in which case the order parameteramplitudes can be denoted by η1 = η and η2 = ηeiφ, where φ is the SC phase difference.The free energy then becomesF = F0 + ϵη4 cos (2φ), (2.2)where F0 is independent of φ. For ϵ > 0, free energy minima occur when φmin = ±π/2,which indicates a time-reversal symmetry breaking px ± ipy order.2.3 ModelMoving towards a concrete model, we begin by noting that a 2D p-wave superconductor canbe engineered using an array of proximitized semiconductor quantum wires. For simplicity,we represent each wire as a Kitaev chain [14]. In Appendix A.2, we analyze a more realistic52.3. ModelFigure 2.2: Kitaev wire array model and its phases. (a) Single layer geometry withtunneling amplitudes and SC pairing indicated. (b) Phase diagram of the bilayer at theoptimal phase difference φ = π/2 with ∆ = 0.8 and t′ = 0.3. The color of the markersindicates the Chern number as shown in the legend, while their size shows the magnitudeof minimum gap in arbitrary units. Grey corresponds to points where C could not bedetermined reliably. (c) Fermi surface evolution with increasing chemical potential µ fort′ = 0.3 and g = 0.4. (d) The energy spectrum for an infinite strip geometry in C = 1 andC = 2 phases. The color scale indicates the normalized expectation value of the positionoperator along the finite direction of the strip. Note that only a section of the Brillouinzone around k = π, where the edge states cross the gap, is shown. Parameters used are∆ = 0.2, t′ = 0.3, and (g, µ) = (0.8,−2.4) and (0.5,−1.8) in the left and right panels,respectively.model which includes spin-orbit coupling and the Zeeman field effect. Both models showthe phenomenology outlined above.A Kitaev chain is characterized by an inter-site tunneling amplitude t, a nearestneighbor pairing field ∆ and an on-site chemical potential µ. Consider chains that arestacked next to each other, parallel to the x axis. If they are coupled via an inter-wiretunneling amplitude t′, one obtains a square lattice as shown in Fig. 2.2(a). It is easy tosee that the resulting 2D superconductor has a gap function with px symmetry.Similarly, one can imagine constructing another monolayer by stacking wires parallelto the y axis to produce a py superconductor. This layer may be aligned and stacked atopthe first layer as shown in Fig. 2.1. The coupling between the two layers can be modeledby an inter-layer tunneling g that connects the two sites in the square lattice unit cell. Inreal space, the lattice model of the bilayer readsH =∑r,a[−tc†r,acr+δ,a − t′c†r,acr+δ′,a +∆c†r,ac†r+δ,a −µ2c†r,acr,a − gc†r,1cr,2 + h.c.] (2.3)where δ is the unit vector parallel to the individual wires in layer a, while δ′ is perpendicular62.3. Modelto it. Henceforth, we will work in units where t = 1.An infinite periodic system can be described by defining a four-component Nambuspinor Ψk = (ck1, c†−k1, ck2, c†−k2)T , where c†ka creates a fermion with momentum k in layera. In this notation, we have H = 12∑kΨ†khkΨk + E0 with the Bogoliubov-de Gennes(BdG) Hamiltonianhk =ξk1 ∆k1 g 0∆∗k1 −ξk1 0 −gg 0 ξk2 ∆k20 −g ∆∗k2 −ξk2 , (2.4)and E0 =∑k(ξk1 + ξk2). Therein,ξk1 = −(2t cos kx + 2t′ cos ky + µ),ξk2 = −(2t cos ky + 2t′ cos kx + µ), (2.5)∆k1 = 2∆sin kx, ∆k2 = 2∆eiφ sin ky.While the physical time-reversal symmetry T that reverses the spin of the electron isbroken in quantum wires due to the applied magnetic field B, the Kitaev chain [14] andour BdG Hamiltonian (2.4) with φ = 0 obey an antiunitary symmetry T ′ that can beregarded as time reversal symmetry for spinless fermions (T ′2 = 1). It is generated byτ zh∗kτz = h−k, where Pauli matrices τα act in Nambu space. This is in addition to thecharge conjugation symmetry K generated by τxh∗kτx = −h−k. Formation of a px ± ipystate in the bilayer structure is marked by a spontaneous breaking of this residual timereversal T ′.Diagonalization of (2.4) produces two pairs of bands ±E1,2k and the associated Blocheigenstates |ψ1,2k ⟩. The ground state energy, computed as Eg = E0 −∑k(E1k + E2k), cannow be minimized with respect to the phase difference φ between order parameters ofthe two layers. We find that two degenerate minima exist at φmin = ±π/2 for the entirerange of system parameters considered, in accordance with the GL analysis with ϵ positive.When the spectrum is gapped, one can characterize the topology of the bulk bands usingthe Chern number C (see App. A.1). It can be conveniently calculated from the Blocheigenstates using the gauge-invariant formulation of the Berry curvature given in Ref. [70].The resulting phase diagram is shown in Fig. 2.2(b).We observe that the model shows five distinct phases with C ranging from −2 to+2. While the detailed position of the phase boundaries depends on model parameters,the phase diagram topology remains robust. The origin of each of these phases can beattributed to the nature of the occupied bands in the normal state, which can be determinedfrom the Fermi surface (FS). Starting on the electron doped side to the extreme left inthe phase diagram Fig. 2.2(b), at a finite interlayer coupling say g = 0.4, the bands areempty and, hence, C = 0. When the value of µ is above the lowest point of the bottomband, a single electron-like FS is formed as illustrated in Fig. 2.2(c). In the SC state, thesystem can be viewed as an effectively single-layer px + ipy SC, which is known to exhibitC = 1. As µ is tuned up, another electron-like FS develops and each layer behaves as apx + ipy superconductor with an aggregate C = 2. Further increase in µ gives rise to oneelectron-like FS plus one hole-like FS. The electron-like and the hole-like bands contribute72.4. Majorana zero mode in a vortex760 780 800 820 840eigenstate index-0.4-0.200.20.4Energy no vortexn = +1 vortexn = −1 antivortexmu=-1.7, del=0.5, g=0.8, tp=0.2Out[!]=Out[!]=bulk gapFigure 2.3: Energy levels of the lattice model in a square sample with 20×20 unit cells andopen boundary conditions. The C = 1 phase is shown with parameters g = 0.8, t′ = 0.2,µ = −2.4 and ∆ = 0.25. Results with n = ±1 vortex placed at the central plaquetteare compared with n = 0 vortex-free configuration. Insets show wavefunction amplitudesassociated with two near-zero modes of the n = +1 vortex solution.1 and −1 respectively in the SC spectrum, resulting in the C = 0 phase seen to occupythe central region of the phase diagram. Continuing this line of reasoning, one eventuallyhas two and one hole-like FSs with C = −2 and C = −1, respectively, for positive µ, andfinally when the FS vanishes again C = 0. As expected, when the topological invariant isnon-zero, in a strip geometry the spectrum exhibits protected edge modes that traversethe gap and connect the bulk bands as shown in Fig. 2.2(d).2.4 Majorana zero mode in a vortexIn addition to protected edge states, a px ± ipy superconductor is expected to host a singleMajorana zero mode at the core of an Abrikosov vortex [67]. As the prediction applies toa single-layer px ± ipy SC, one might anticipate that in a bilayer geometry MZMs fromthe two layers would hybridize and form a complex fermion at non-zero energy. For thisreason, we expect vortices to host MZMs only in C = ±1 phases of our model whichbehave, effectively, as single-layer px ± ipy superconductors.To test these conjectures we solve the problem defined by the Hamiltonian (2.3) on asquare domain with L× L sites, open boundary conditions and a single vortex placed atthe origin. The vortex is encoded by attaching a position dependent phase factor to theSC pair field according to∆rδ = ∆e−inθrδ . (2.6)Here ∆rδ is the order parameter living on the bond connecting lattice sites r and r + δ,while θrδ denotes the angle between the line connecting the midpoint of this bond to thevortex center and the x axis. Integer n denotes the vorticity. Details of this calculationand additional remarks are given in Appendix A.3.Fig. 2.3 shows our results in the C = +1 topological phase. Even without the vortex, we82.5. Outlookobserve a number of states with energies inside the bulk SC gap – these are the protectededge modes mandated by the bulk-boundary correspondence. In the n = 1 vortex case,two modes appear near zero energy. Inspection of the associated eigenstates reveals thatthey are composed of linear superpositions of an edge mode and a mode localized atthe vortex center (inset Fig. 2.3). Note that the model at hand only has C2 rotationsymmetry and hence the eigenstates also have this symmetry. The small energy splittingϵ0 can be attributed to the finite size effect, specifically the overlap between the zeromode wavefunction and the corresponding edge mode. As expected we find that ϵ0 decaysexponentially with increasing system size L. Therefore, we identify these states as theMajorana zero modes bound to the vortex and the system edge, respectively.Analytical solution to the simple model [67] indicates that for the px + ipy case MZMis present in the vortex (n = +1) but absent in the antivortex (n = −1). For px − ipy thesituation is reversed. In our numerical solution, n = −1 antivortex does not show anynear-zero modes, in accordance with the above expectation. Likewise, we find no zeromodes in |C| = 0, 2 phases.2.5 OutlookFollowing the original signatures of Majorana zero modes in proximitized InSb wires[27], great strides have been made in the preparation and characterization of high-qualityquantum wires and their assembly into complex structures [66]. As a result, it is nowpossible to fabricate T- and Y-junctions between individual wires as well as assemble wiresinto arrays with various geometries. In light of these developments, it might be possible tofabricate and probe the setup discussed here.A key step in this direction is the recently demonstrated ability [71] to apply the SCcoating on one side of the wire, as schematically indicated in Fig. 2.1. For our proposeddevice this feature is essential as we require single-electron tunneling between the layers todominate. While p-wave pairs cannot tunnel in the 90◦ twisted geometry, s-wave Cooperpair tunneling is allowed. Should the coating material on wires from different layers comeinto contact, such tunneling would lock the SC phases of the two layers together and drivethe system away from the desired px ± ipy state, which requires φ = ±π/2.In the absence of naturally occurring chiral px ± ipy superconductors, our constructionoutlines a new path towards artificially engineered structures that behave effectively assuch. The Fu-Kane model [28] is often quoted as a realization of the px ± ipy SC, butit is important to note that although both harbor MZMs in vortices, they represent twofundamentally different quantum states: A Fu-Kane SC respects T , while the chiral p-wavestate breaks it. Further, the former cannot be realized as a purely 2D system, in contrastto the latter. It follows that even if FeTexSe1−x is eventually confirmed as a realizationof the Fu-Kane paradigm, as seems likely based on existing evidence [57–60], search for agenuine T -broken px ± ipy superconductor in natural and artificially engineered systemswill remain a valuable pursuit.Our architecture relies on quantum wires, which themselves host MZMs. In the basicversion of such wires [27], MZMs are permanently bound to their endpoints and are thusimmobile. By contrast, in the 2D structure discussed here MZMs are attached to Abrikosov92.5. Outlookvortices that are intrinsically mobile objects. The mobility of vortices in two dimensionscould assist the long-term goal of verifying their predicted non-Abelian exchange statistics,which is the foundation for topologically protected, fault-tolerant quantum computation.10Chapter 3Chiral d-wave superconductivity intwisted bilayer cuprates3.1 IntroductionA key take away from the previous chapter is that when two monolayers of a nodalsuperconductor are stacked with a relative twist, spontaneous time-reversal symmetrybreaking could give rise to topological superconductivity. We had to artificially assemblea p-wave superconductor, but one form of 2D nodal superconductivity that does occurnaturally is the d-wave state. In a remarkable recent development Yu and co-workers [72]succeeded in isolating and probing monolayer samples of high-temperature superconductorBi2Sr2CaCu2O8+δ (Bi2212) with critical temperatures Tc essentially indistinguishable fromthe bulk Tc ≃ 90K of optimally doped Bi2212 crystals. This result delivers the long-awaitedproof that superconductivity in Bi2212 (and possibly all the cuprates) is an intrinsicproperty of the 2D layers. It also opens exciting possibilities for engineering structureswith novel behaviour and functionality, from individual layers. Motivated by the results ofRef. [72] and recent developments in van der Waals-bonded materials [32, 43, 73–76], weconsider here structures composed of two cuprate monolayers with a twist angle θ betweenthem.Twisted bilayer graphene famously exhibits a sequence of gate-voltage controlledcorrelated insulating and superconducting phases [32, 73–76] when θ is close to the magicangle 1.1o. Similar physics is reported [43] in twisted bilayer WSe2 in the range of anglesbetween 4.0o and 5.1o. In both cases, the observed behaviour is thought to reflect asubtle interplay between the band structure energetics and topology brought about bythe moire´ patterns, electron-electron interactions and disorder [31, 77–85]. As such, manyaspects of the reported experimental observations in twisted graphene and WSe2 still awaitfull theoretical understanding [86]. The problem of twisted Bi2212 bilayers is potentiallysimpler to analyze theoretically, at least at the level where each Bi2212 layer is treated as aBCS superconductor with a dx2−y2 order parameter. Results of such a theoretical analysisare, nevertheless, non-trivial: over a continuous range of twist angles θ in the vicinityof 45o, we find that below a critical temperature T˜c(θ) the bilayer system spontaneouslybreaks time-reversal symmetry T and forms a chiral topological phase that can be thoughtof as an emergent dx2−y2 ± idxy superconductor (d± id′ for short). Such a superconductorhas a fully gapped bulk and is characterized by Chern number C = ±2 per layer with a gapof dx2−y2 ± idxy symmetry, respectively. It exhibits two gapless, topologically protectedchiral edge modes per monolayer. These modes are coherent superpositions of electron andhole states and can be therefore classified as Majorana particles [65]. When θ ≈ 45o, our113.2. Ginzburg-Landau theoryresults indicate that T˜c ≃ Tc, suggesting that, remarkably, in carefully assembled Bi2212samples, the topological phase could set in well above the liquid nitrogen temperature.In momentum space, the dx2−y2 order parameter that characterizes each monolayercan be represented as ∆k = ∆0(kˆ2x − kˆ2y), where ∆0 denotes the amplitude and (kˆx, kˆy) =(kx, ky)/k are momentum components defined relative to the principal crystal axes of theunderlying copper-oxygen lattice. When one of the layers is rotated by angle θ, its orderparameter is transformed in the original coordinate frame to∆(θ)k = ∆0[cos(2θ)(kˆ2x − kˆ2y) + sin(2θ)2kˆxkˆy]. (3.1)From the vantage point of the unrotated monolayer, the second term in Eq. (3.1) hasthe functional form of a dxy order parameter. For coupled layers, we then expect theCooper pair tunneling to induce a subdominant dxy order parameter in the unrotatedlayer. By symmetry, the same should happen in the other layer. Furthermore, if the twoorder parameters combine with a complex phase, the result is a fully gapped topologicalsuperconductor with broken T . The gapped nature of the phase follows from the factthat the excitation gap is given by |∆k| which, for a d± id′ gap function, is non-vanishingeverywhere at the Fermi level, see Fig. 3.1(a,b). T -breaking is evident because time reversalmaps d+ id′ to d− id′. We demonstrate that for a range of twist angles θ between θ−c andθ+c , such a T -broken phase becomes energetically favorable below a critical temperatureT˜c(θ).T -broken d± id′ phases have been discussed theoretically in various contexts, includingbulk high-Tc cuprates in applied magnetic field [87–89], c-axis twist Josephson junctions [23,90–92], and highly doped single-layer graphene [24, 93]. However, conclusive experimentalevidence for these states has never been reported. By placing emphasis on monolayer-thinsamples of an established dx2−y2 superconductor, we argue that for ≃ 45o twist, theT -breaking d± id′ phase is nucleated robustly and reliably through the simple and wellunderstood mechanism of interlayer Cooper pair tunneling that is largely insensitive tomicroscopic details.3.2 Ginzburg-Landau theoryThe phenomenon outlined above is most easily quantified by examining the phenomeno-logical Ginzburg-Landau (GL) theory. For the sake of clarity, we begin by consideringmodels with a single CuO2 plane per monolayer. Though directly relevant to Bi2Sr2CuO6+δ(Bi2201), our results are broadly applicable to cuprates with multiple CuO2 plane structure,such as Bi2212, which has been exfoliated [72] to the monolayer limit.The GL free energy density for two coupled monolayers can be written asF [ψ1, ψ2] = f0[ψ1] + f0[ψ2] +A|ψ1|2|ψ2|2 +B(ψ1ψ∗2 + c.c.) + C(ψ21ψ∗22 + c.c.), (3.2)where ψa=1,2 are complex order parameters of the two layers and f0[ψ] = α|ψ|2 + 12β|ψ|4 isthe free energy of a monolayer. We are interested in spatially uniform solutions, so thegradient terms have been omitted. If the two layers are physically identical, then the order123.2. Ginzburg-Landau theoryπ/40 π/21.00.50.0π0-π1.00.50.0-0.5θφθ = 9.53o < θ−c θ = 36.87o > θ−cθ θ+cθ−cT/TcTSCϵk ϵkk⊥ k⊥k∥ k∥a) b) c)θ = 0θ = π /4θ = θ−cdx2−y2 + idxydx2−y2 + dxynormaldSC dSCℱ(φ)d)φ minf)π/40 π/20π/2π1020θgap[meV]TSCBdGGLe)θkxkyθ−c θ+cFigure 3.1: Time reversal symmetry breaking in twisted bilayers. Panel (a) illustratesthe lattice geometry at small twist angle θ < θ−c where the free energy is minimized forinterlayer phase difference φ = 0. This results in a gapless spectrum with weakly splitDirac points. For twist angles θ > θ−c in panel (b) the preferred state has φ ̸= 0. It breakstime reversal and is fully gapped as indicated by the massive Dirac spectrum. The spectraare obtained by diagonalizing the BdG Hamiltonian (3.10); k⊥ and k∥ denote momentumcomponents perpendicular and parallel to the Fermi surface near a nodal point. (c) showsthe GL free energy F(φ) for twist angles θ ranging from 0o to 45o in 3o increments forK = 0.125. When θ > θ−c , minima occur away from φ = 0. (d) Predicted phase diagrambased on GL theory (3.2), where T˜c(θ) is given by Eq. (3.6). TSC denotes the gappedtopological phase, dSC stands for gapless d-wave superconductor. Panel (e): Fermi surfaceand d-wave gap in the continuum formulation of the microscopic model. (f) displays therelative phase between the layers φmin that minimizes the GL and the BdG free energy asa function of twist angle θ. Orange symbols represent the spectral gap.parameter amplitudes must be the same and the most general solution (up to an overallphase) isψ1 = ψ, ψ2 = ψeiφ, (3.3)where we take ψ to be real and positive.Because ψa describe order parameters with d-wave symmetry, they transform asψa → −ψa under a 90o rotation. Suppose we increase the twist angle θ continuouslyfrom 0 to π/2. Since ψ2 changes sign under this operation, we conclude that for the freeenergy (3.2) to remain invariant, the parameter B must also change sign. The simplestdependence consistent with this constraint is B = −B0 cos(2θ), where we take B0 > 0. Thenegative sign is chosen to allow the phase difference between the layers to vanish as thetwist angle approaches zero. This also implies that B = 0 for θ = π/4, consistent with theintuition that a dx2−y2 Cooper pair cannot tunnel between the layers in this configuration.Parameters A and C are not required by symmetry to depend on θ, so we treat them asconstants henceforth. Although coefficient A does not have a bearing on the phase, it133.3. Microscopic modelsis expected to be either negative or a small positive quantity. To see this, consider anuntwisted multilayer stack: a large positive A would translate to a vanishing condensateamplitude in every other layer. Such a scenario contradicts experimental findings and,consequently, justifies our hypothesis. Incorporating these ingredients into Eq. (3.2), wefindF(φ) = F0 + 2B0ψ2 [− cos(2θ) cosφ+K cos(2φ)] , (3.4)where K = Cψ2/B0 and F0 contains all terms independent of φ. It is now a simple matterto minimize F(φ) for a given twist angle θ. An interesting scenario arises when K > 0,as indicated by Fig. 3.1(c). Although symmetry alone does not constrain the sign ofK, a simple physical argument suggests that K should, in general, be positive. One caninterpret the ψ21ψ∗22 term in Eq. (3.2) as representing the coherent tunneling of two Cooperpairs between the layers. Up to an overall scale factor, its coefficient should, therefore,be proportional to the square of the coefficient B, which itself corresponds to single pairtunneling. Coefficient C (and, hence, also K) will then be positive, regardless of the sign ofB. We demonstrate that microscopic models fully support this argument, yielding K > 0.At θ = 45o and for K > 0, the free energy minimum clearly occurs at φmin = ±π/2which indicates a d± id′ superconducting state. For θ ≠ 45o, the free energy is minimizedbyφmin = arccos(cos 2θ4K). (3.5)As illustrated in Fig. 3.1(f), a nontrivial phase difference occurs for a range of anglesθ−c < θ < θ+c with θ±c =12 arccos (∓4K). This indicates a T -broken phase with orderparameter dx2−y2+e±iφmindxy which is fully gapped and topologically non-trivial as long asφmin ̸= 0, π, see also Fig. 3.1(a,b). As K depends on temperature through its dependenceon ψ, the critical temperature T˜c of the topological phase will be a function of the twistangle θ. If we adopt the standard GL temperature dependence for the order parameterψ(T ) = ψ0√1− T/Tc, it is straightforward to deduceT˜c(θ) = Tc(1− | cos 2θ|4K0), θ−c < θ < θ+c , (3.6)where K0 = Cψ20/B0. This defines the phase diagram in Fig. 3.1(d). We observe that,remarkably, T˜c(45o) coincides with the Tc of a single monolayer which can be as high as 90Kin carefully prepared Bi2212 flakes. Away from the ideal 45o twist, the critical temperaturefalls approximately as T˜c(θ) ≃ Tc(1− |θ − π/4|/2K0). Microscopic models with realisticBi2212 parameters discussed in the following give typical values K0 ≃ 0.1− 0.2, implying asignificant extent for the topological phase, as indicated in Fig. 3.1(d).3.3 Microscopic modelsWe now turn to the microscopic theory. Although the pairing mechanism in high-Tccuprates remains a subject of debate, it is widely accepted that most physical propertiesof the superconducting state in the optimally doped and overdoped regime are accuratelydescribed within the framework of the standard BCS theory with a dx2−y2 order parameter.143.3. Microscopic models0.0 0.2 0.4 0.6 0.8 1.0 1.202468min (meV)= 53.13o= 43.60o0.0 0.2 0.4 0.6 0.8 1.0 1.2T / Tc0.00.10.20.30.40.5= 53.13o= 43.60o010203040max (meV)θ2,5 = 43.600a) b) d) c) Figure 3.2: Lattice model results. Panels (a,b) show the temperature dependence of theminimum gap ∆min, the maximum gap ∆max and phase φ, based on a self-consistentlattice calculation for coupled layers with commensurate twist angles θ1,2 ≃ 53.13o, andθ2,5 ≃ 43.60o corresponding to a unit cell with 10 and 58 sites, respectively. Panel (c)shows zero temperature phase diagrams of the system at θ2,5 twist as a function of chemicalpotential µ and interlayer coupling g0. The range of chemical potentials µ ∈ (−1.6t,−1.0t)correspond to fillings n ∈ (0.033, 0.04) near optimal doping. Each data point in the phasediagram corresponds to an independent self-consistent solution and the radius of markers isproportional to the size of the minimum gap ∆min. The color indicates the Chern numberas shown in the legend. Panel (d) illustrates the edge modes for C = 4 topological phase,in θ2,5 configuration for parameters µ = −1.3t and g0 = 20meV. The energy spectrum isshown for a bilayer system in the infinite strip geometry with a 90 unit cell width. Thecolor scale represents the normalized position expectation value of the eigenstate along thedirection of finite length.Therefore, we begin by modeling the twisted bilayer using a simple continuum modelof coupled d-wave superconductors, and complement this with a calculation based onan attractive Hubbard model on the square lattice. The continuum formulation has theadvantage of being applicable to an arbitrary twist angle θ, similar to the Bistritzer-MacDonald model for graphene [31]. In contrast with the continuum model above, thelattice calculation captures various microscopic details of CuO2 planes more effectively,such as the shape of the Fermi surface. It is, however, limited to commensurate twistangles that produce relatively small moire´ unit cells due to practical constraints.In keeping with the strategy outlined above, we consider a single CuO2 plane permonolayer, as for Bi2201. The continuum theory for such a bilayer is defined by theHamiltonianH =∑kσaξkac†kσackσa + g∑kσ(c†kσ1ckσ2 + h.c.)(3.7)+∑ka(∆kac†k↑ac†−k↓a + h.c.)−∑ka∆ka⟨c†k↑ac†−k↓a⟩.Here, c†kσa creates an electron with crystal momentum k and spin σ in layer a = 1, 2153.3. Microscopic modelswhile ξka and g represent the in-plane kinetic energy and inter-plane tunneling amplitude,respectively. The superconducting order parameter in layer a can then be expressed as∆ka =∑p′V(a)kp ⟨c−p↓acp↑a⟩. (3.8)The prime on the summation indicates a restriction to momentum states with energywithin ϵc of the Fermi level. The Hamiltonian (3.7) should be regarded as a mean-fieldapproximation to the BCS pairing Hamiltonian with an interaction term∑kp V(a)kp c†k↑ac†−k↓ac−p↓acp↑a, where V(a)kp denotes the interaction matrix element in layer a. We shall use asimple separable formV(a)kp = −2V cos (2αk) cos (2αp), (3.9)where αk represents the polar angle of the vector k. This is known to yield a robustsolution with dx2−y2 symmetry for a single CuO2 layer, namely ∆ka = ∆d cos (2αk).In order to incorporate the twist, we take the interaction in layer 1 as in Eq. (3.9),but we rotate the interaction in layer 2 by angle θ such that V(2)kp = −2V cos (2αk − 2θ)cos (2αp − 2θ). For the sake of simplicity, we consider a circular Fermi surface generatedby ξka = ℏ2k2/2m− µ that remains invariant under rotation, see Fig. 3.1(e). The problemposed by Hamiltonian (3.7) is then solved by defining a four-component Nambu spinorΨk = (ck↑1, c†−k↓1, ck↑2, c†−k↓2)T in terms of which H = ∑kΨ†khkΨk + E0. Here, theBogoliubov-de Gennes (BdG) Hamiltonian readshk =ξk ∆k1 g 0∆∗k1 −ξk 0 −gg 0 ξk ∆k20 −g ∆∗k2 −ξk (3.10)and E0 =∑k 2ξk −∑ka∆ka⟨c†k↑ac†−k↓a⟩. Diagonalizing hk gives two pairs of energyeigenvalues ±Ekα (α = 1, 2) for each momentum k. With the Hamiltonian (3.7) expressedin diagonal form, the free energy of the system can be calculated from the standardexpressionFBdG = E0 − 2kBT∑kαln [2 cosh (Ekα/2kBT )] . (3.11)By performing a systematic expansion of FBdG in powers of the order parameter amplitudes∆d (see Appendix B.1), it is possible to ascertain various coefficients entering the GL freeenergy (3.2). This confirms the form of the GL coefficients B and C deduced previouslyon the basis of general symmetry arguments, and that K in Eq. (3.4) is positive.We proceed with a fully self-consistent solution which follows from minimizing FBdGwith respect to ∆ka and can be performed for any given V , twist angle θ and temperatureT . To facilitate this calculation we assume that the order parameters in the two layers havethe same amplitude ∆d, but can differ in phase as in Eq. (3.3). Free energy computed fromEq. (3.11) then shows the same qualitative behaviour as the GL theory (3.4) with K > 0:for small twist angles, FBdG(φ) has a single minimum at φ = 0, while for twist anglesclose to 45o there are always two degenerate minima at φ = ±φmin, indicating formation163.3. Microscopic modelsof a T -broken topological phase. Fig. 3.1(f) shows the calculated φmin as a function ofthe twist angle θ for realistic Bi2212 parameters [94–96] with g = 30meV at T = 0. In aBCS superconductor, because phonons enable the formation of Cooper pairs, the energycut-off ϵc is set by the Debye frequency. In cuprates, however, electronic correlations arebelieved to drive the superconducting transition, so Fermi energy is the natural scale. Weuse ϵc = 60meV to ensure an experimentally comparable gap amplitude ∆d = 40meV.The self-consistent calculation is qualitatively insensitive to the choice of ϵc. An excellentagreement with the GL theory prediction is observed when K = 0.125. The BdG theoryalso allows us to predict the minimum excitation gap which, in this simple model, can be aslarge as 20meV when θ is close to 45o. In generic multi-band superconductors that breakT , inter-band pairing can give rise to Fermi surfaces of Boguliubov quasiparticles [97].In our model, that would correspond to inter -layer pairing. The literature on cupratesprovides a good reason to believe that realistic mechanisms cannot bind electrons intoCooper pairs across layers.The BdG theory can be formulated directly on the lattice by starting from the Hubbardmodel with nearest-neighbour attraction. This is known to produce a dx2−y2 superconductorwhen applied to a single CuO2 monolayer. The Hamiltonian isH = −∑ij,σatijc†iσacjσa − µ∑iσaniσa +∑ij,aVijnianja −∑ijσgijc†iσ1cjσ2, (3.12)where tij encodes the normal-state band structure of the single layer, gij describes theinterlayer tunneling and Vij denotes density-density interactions. The mean-field calcula-tions are performed in bilayer geometries characterized by a twist vector v = (m,n) andcommensurate twist angle θm,n = 2arctan (m/n) as explained in Appendix B.2. Fig. 3.2summarizes our main results. These lend further support to our conclusions drawn fromthe GL and continuum BdG approaches.The lattice model confirms the onset of the T -breaking phase below critical temperatureT˜c(θ) whose dependence on θ is consistent with the GL prediction, see Fig. 3.2(a,b). Theresults are in reasonable quantitative agreement with the continuum BdG theory but thelattice model gives a smaller spectral gap. This may be attributed to more complicatedFermi surface geometry resulting from Brillouin zone folding that accompanies the largereal-space Moire´ unit cell. The lattice model allows for a direct evaluation of the Chernnumber C. In addition to the C = 4 topological phase anticipated on the basis of thecontinuum BdG theory, the T = 0 phase diagram in Fig. 3.2(c) reveals the existence ofC = 2 and C = 0 gapped phases that can be reached by varying the chemical potential µand the interlayer coupling strength g0. In an experiment, the former can be tuned over awide range by oxygen annealing [72], while the latter may depend on the twist angle θ andcould be further controlled by applying hydrostatic pressure as demonstrated in twistedgraphene [74]. The C = 0 phase corresponds to a d+ is superconductor, while the C = 2phase can be effectively thought of as a single-layer d+ id′ superconductor. For a stripgeometry, the lattice model predicts protected chiral edge modes traversing the bulk gap,Fig. 3.2(d), confirming the topological nature of the phase.173.4. Discussion3.4 DiscussionOur results establish a new avenue for engineering high-temperature topological super-conductors by combining monolayers of van der Waals-bonded d-wave superconductingcuprates, such as Bi2Sr2CuO6+δ or Bi2Sr2CaCu2O8+δ [72]. The topological phase is fullygapped, breaks time-reversal symmetry and supports protected chiral Majorana modesat the sample boundary. The chiral d± id′ phase arises robustly, driven by Cooper pairtunneling between the layers when they are twisted by an angle close to 45o. Unlike ingraphene, where the twist angle must be tuned very accurately, we predict that a widerange of angles (45o ± 10o) will produce nontrivial physics in these systems.The resulting fully gapped topological phase can be probed by conventional spectroscopicor transport techniques. For example, on lowering the sample temperature, angle resolvedphotoemission or tunneling spectroscopies [95, 96] will detect the transition at T = T˜c(θ)from the gapless dx2−y2 superconductor with Dirac nodes to a fully gapped excitationspectrum. This change should also be observable via thermodynamic probes such asthose measuring the specific heat [98] and superfluid density [99, 100]. Signatures of thebroken T -symmetry characteristic of the d± id′ and d+ is phases may be identified bycircular birefringence and dichroism experiments [101]. In transport, the non-zero Chernnumber manifests as quantized electronic contribution to the thermal Hall conductanceκxy = C(π2k2BT/3h) [89]. In mesoscopic samples, longitudinal thermal conductance κxxwill also be quantized due to the ballistic contribution of the edge states. The d ± id′phase exhibits anomalous interlayer Josephson effect due to the two-minima structureof the free energy depicted in Fig. 3.1(c); the Josephson current IJ = (2e/ℏ)∂F/∂φ willdevelop pronounced deviations from the usual sinusoidal dependence on the phase bias∆φ = φ− φmin.It is worth noting that existing experimental results on c-axis twist Josephson tunnelingin cuprates already offer some support for the ideas advanced here. In the following chapter,we discuss the signatures of time-reversal symmetry breaking in transport experiments.18Chapter 4Josephson effects in twistedcuprate bilayers4.1 IntroductionDissipationless current driven by a phase difference between superconductors – the Joseph-son effect [102] – constitutes a quintessential manifestation of a coherent quantum state atmacroscopic scales. Measurements of the critical current Jc as a function of temperature T ,magnetic field B and applied electromagnetic radiation provide important insights into thephenomenology and microscopic mechanisms of a great variety of superconductors rangingfrom conventional metals to cuprate and iron-based superconductors. In this chapter, wediscuss various aspects of the Josephson effect between two samples of high-Tc cupratesstacked with a relative twist angle θ, as illustrated in Fig. 4.1(a). It is well known thatcuprates are d-wave superconductors [103], implying a strongly anisotropic order parameter.One thus expects significant dependence of the Josephson current on the twist angle, whichwas indeed noted in the classic works on this topic [90, 104–108]. For reasons that remainnot well understood, this expectation was not borne out in early [109, 110] and more recent[49] experimental works which showed twist-independent critical currents (see however Ref.[111], for a contrary result on cross-whisker junctions). Very recent experimental work ontwisted thin Bi2Sr2CaCu2O8+δ (Bi2212) flakes assembled at cryogenic temperatures [50]finally reported the critical current anisotropy indicative of d-wave order parameter.The experimental finding [72] that Bi2212 superconducts up to Tc ≃ 90K whenexfoliated down to a single monolayer inspired renewed interest in twisted cuprate bilayers.In addition to the nontrivial topology near the 45◦ twist we discussed previously, theoreticalstudies also predict the emergence of novel spontaneously T -broken phases in the vicinityof ‘magic angle’ θM [48] that depends on model parameters but is generally in the rangeof few degrees. The magic angle in cuprates is marked by a pair of Dirac cones merginginto a quadratic band crossing. In analogy to twisted bilayer graphene [31, 32, 73–76],interactions become important near θM and can also lead to spontaneous T breaking.On the other hand, T -breaking near 45◦ arises from the co-tunneling process (coherenttunneling of two Cooper pairs) when ordinary single-pair tunneling is suppressed by thed-wave symmetry of the pair wavefunction. A typical phase diagram, predicted on the basisof continuum Bogoliubov-de Gennes (BdG) theory, is shown in Fig. 4.1(b). It illustratesthe central theoretical finding that at a 45◦ twist the topological T -broken d+ id′ phase hasa potential to persist up to high temperatures, approaching the native Tc of the cupratematerial.The goal of this chapter is to elucidate how these T -broken phases are manifested in194.1. Introduction0 15 30 45 60 75 900.00.51.0T/Tc0θdSCd + isd + id′ θMThe complete phase diagram of a bilayer ec,v1,v2,v0,g=60,.12,.0,.043,10.5 (meV)θJBℏωa)b)0 15 30 45 60 75 900.00.51.0T/Tc0θdSCd + isd + id′ θMThe complete phase diagram of a bilayer ec,v1,v2,v0,g=60,.12,.0,.043,10.5 (meV)θJBℏωa)b)0 15 30 45 60 75 900.00.51.0T/Tc0θdSCd + isd + id′ θMThe complete phase diagram of a bilayer ec,v1,v2,v0,g=60,.12,.0,.043,10.5 (meV)θJBℏωa)b)Figure 4.1: (a) Schematic of the Josephson effect measurement on a twisted cuprate bilayer.Magnetic field B parallel to the interface and electromagnetic radiation with frequencyω can be used to probe the current-phase relationship. (b) Typical phase diagram ofthe system calculated from the microscopic model discussed in Sec. 4.2. It is assumedthat d+ is order is nucleated near the magic angle θM with the s component driven bya weak on-site attractive interaction. Model parameters are ϵc = 60meV, g = 10.5meV,NFV = 0.12 and NFVs = 0.043.the Josephson effect and, specifically, what type of measurement can deliver solid evidenceof their existence. We also address the fate of the T -broken phases in twisted structurescomposed of thicker flakes which might be easier to assemble in the laboratory. Our mainfinding in this regard is that the basic phenomenology of T -breaking remains in placealthough the topological gap is strongly suppressed as the number of monolayers comprisingthe flake increases.Much of the phenomenology underlying the Josephson effect can be understood fromthe Ginzburg-Landau theory of Sec. 3.2, which we reproduce here for completeness. Fortwo weakly coupled monolayers, the corresponding GL free energy takes the formF [ψ1, ψ2] = F0[ψ1] + F0[ψ2] +A|ψ1|2|ψ2|2 +B(ψ1ψ∗2 + c.c.) + C(ψ21ψ∗22 + c.c.), (4.1)where ψa are complex scalars representing the dx2−y2 order parameters in layers a = 1, 2and F0[ψ] = α|ψ|2 + 12β|ψ|4 is the free energy of a monolayer. Physically, the B and Cterms represent coherent tunneling of single and double Cooper pairs between the layers,respectively.For two identical layers, the order parameters can only differ by a phase which allowsus to write ψ1 = ψ and ψ2 = ψeiφ, where we take ψ to be real and positive. d-wavesymmetry additionally implies that GL parameter B must change sign when the twist isincreased by π/2. We henceforth assume the simplest angle dependence consistent withthis condition, B = −B0 cos(2θ), where B0 > 0 is taken to ensure that untwisted layersare at zero phase. It is also generally true that C > 0. With these ingredients, we maywrite the corresponding Josephson free energyF(φ) = E0 − ℏ2e[Jc1 cosφ− 12Jc2 cos(2φ)], (4.2)204.1. Introductionwhere Jc1 = (4eB0/ℏ)ψ2 cos(2θ), Jc2 = (2eC/ℏ)ψ4 and E0 collects terms that are indepen-dent of φ. The competition between the cosφ and cos 2φ terms in Eq. (4.2) underlies theemergence of the spontaneously T -broken phase near the 45◦ twist. When θ is close to zerothe conventional Josephson tunneling term Jc1 dominates and the free energy minimumoccurs at φ = 0. Increasing the twist, however, decreases Jc1 ∼ cos(2θ). Eventually,for twist angle approaching 45◦, the Jc2 term begins to dominate, and F(φ) developstwo distinct minima at ±φmin signalling the T -broken phase. The other T -broken phasethat occurs near the magic angle θM depends on the structure of low-energy quasiparticleexcitations and interaction physics. As such it cannot be understood based on the simpleGL theory formulated above.The equilibrium current between the layers follows from the Josephson relation J(φ) =(2e/ℏ)(dF/dφ), which yields a simple but all-important current-phase relationJ(φ) = Jc1 sinφ− Jc2 sin(2φ) (4.3)that will form the basis for much of our analysis. Generically, coefficients Jc1,2 depend onboth the twist angle θ and temperature T . Within the GL theory, temperature dependencefollows from expressions below Eq. (4.2) together with ψ(T ) = ψ0√1− T/Tc, whichholds close to Tc. We note that other studies have looked at the general consequencesof a dominant second harmonic in the current-phase relation (cf. Ref. [112]). Since ourparameters are informed by the microscopic model, the results will be more directly relevantto twisted cuprates.The critical current is a physical observable that is straightforward to measure experi-mentally Jc = maxφ[J(φ)]. We find that T and θ dependence of Jc show some interestingfeatures in the low-T regime but, perhaps surprisingly, do not contain any clear signaturesof the T -broken phases; Jc(T, θ) is a smooth function of its arguments across the transitionto the T -broken phase. The transition, however, is signalled by the vanishing phase stiffness,defined asρs =dJ(φ)dφ∣∣∣∣φ=0. (4.4)Mathematically, the condition ρs = 0 marks the point at which the free energy minimumat φ = 0 becomes a local maximum. Unfortunately, ρs is not easily measurable. For theconventional sinusoidal current-phase relationship (i.e. when Jc2 = 0), it is easy to see thatρs and Jc coincide. We will highlight departures from the Jc = ρs equality to quantifydeviations from the conventional Josephson behavior in twisted cuprate bilayers. TreatingJc2 as a constant independent of θ and writing Jc1(θ) = Jc1(0) cos(2θ), Eq. (4.4) can beused to determine the critical twist angle beyond which T breaking occursθc =12arccos(2Jc2Jc1(0)). (4.5)When a magnetic field B is applied parallel to the junction plane, the SC phasedifference between varies spatially and gives rise to the well-known Fraunhofer oscillationswith the critical current given byJc(Φ) = Jc(0)∣∣∣∣sin(NπΦ/Φ0)NπΦ/Φ0∣∣∣∣ , (4.6)214.2. Modelwhere Φ = BS is flux through the effective junction area S and Φ0 = hc/2e is thesuperconducting flux quantum. For an ordinary Josephson relation we have N = 1, butwhen the current-phase relation (4.2) is dominated by the sin(2φ) term N = 2. Whenboth terms are present in the current-phase relation a more detailed analysis, given in Sec.4.7, predicts a crossover from the conventional N = 1 Fraunhofer pattern at small twistangles to the N = 2 behavior as θ approaches 45◦. We conclude that the field dependenceof the critical current can be tested to probe for spontaneous T breaking in twisted cupratebilayers. Similar results are obtained for Shapiro steps in Sec. 4.8, where theory predictsthat fractional steps appear in the current-voltage characteristics when a T broken state issubjected to an external bias in the form of radio-frequency radiation.4.2 ModelTo capture essential physics of the twisted bilayer, we work with a continuum microscopicmodel of coupled d-wave superconducting monolayers [47] which we check against the moreaccurate lattice model. Within the continuum model, it is straightforward to obtain thecurrent-phase relation for any temperature T and twist angle θ. We extract coefficientsJc1,2 and then employ Eq. (4.3) to study Fraunhofer oscillations and Shapiro steps thatcan be experimentally probed to reveal spontaneous T breaking.Assuming periodic boundary conditions, the model is defined by the second-quantizedHamiltonian H =∑kΨ†khkΨk + E0, where Ψk = (ck↑1, c†−k↓1, ck↑2, c†−k↓2)T represents afour-component Nambu spinor and c†kσa creates an electron with spin σ in layer a of thetwisted bilayer. The BdG Hamiltonian is a 4× 4 matrixhk =ξk ∆k1 g 0∆∗k1 −ξk 0 −gg 0 ξk ∆k20 −g ∆∗k2 −ξk (4.7)with normal-state dispersion ξk = ℏ2k2/2m− µ, interlayer coupling g, which conserves thein-plane momentum k, and d-wave pairing amplitudes∆k1 = ∆eiφ/2 cos(2αk − θ),∆k2 = ∆e−iφ/2 cos(2αk + θ). (4.8)Here ∆ is taken as real positive, αk denotes the polar angle of vector k, and ±θ termsencode the twist between the layers. Finally, E0 =∑k 2ξk − V−1∑ka |∆ka|2, where thelast term results from the standard mean-field decoupling of the pairing interaction in thed-wave channel with strength V. The corresponding free energy is given byFBdG = E0 − 2kBT∑kαln [2 cosh (Ekα/2kBT )] , (4.9)where the sum extends over all positive energy eigenvalues Ekα of hk.We are interested in interlayer Josephson current driven by externally imposed phasebias. The order parameter amplitude ∆ will be determined self-consistently through the224.2. ModelJ(φ)J(φ)φφa) θ = 0ob) θ = 44o0 0.2 0.4 0.6 0.8 1T/Tc0.0000.0020.0040.006J cαJc1Jc2Jc3 x100c) θ = 44oJ cαT/0 /2 0.000.050.100.150.200 /2 -0.0020.0000.0020.0040.0060.008Tc0J(φ)J(φ)φφa) θ = 0ob) θ = 44o0 0.2 0.4 0.6 0.8 1T/Tc0.0000.0020.0040.006J cαJc1Jc2Jc3 x100c) θ = 44oJ cαT/0 /2 0.000.050.100.150.200 /2 -0.0020.0000.0020.0040.0060.008Tc0J(φ)J(φ)φφa) θ = 0ob) θ = 44o0 0.2 0.4 0.6 0.8 1T/Tc0.0000.0020.0040.006J cαJc1Jc2Jc3 x100c) θ = 44oJ cαT/0 /2 0.000.050.100.150.200 /2 -0.0020.0000.0020.0040.0060.008Tc0Figure 4.2: (a,b) Typical current-phase relations obtained from the continuum model. Datapoints represented as black dots are obtained by evaluating Eq. (4.13) at fixed twist angleθ for a range of temperatures between zero and just above Tc. The color-coded solid linesrepresent the best fit to the relation Eq. (4.3) amended by a third order term −Jc3 sin(3φ)to achieve a close fit. Model parameters are ϵc = 6 meV, g = 7meV, NFV = 0.12. Panel (c)shows the temperature dependence of the best-fit coefficients Jcα (α = 1, 2, 3) for θ = 44◦.minimization of FBdG for a given bias. To this end, it is convenient to perform a globalgauge rotation (ck1, ck2) → (eiφ/4ck1, e−iφ/4ck2) which moves the phase factor from theorder parameter terms in hk to the interlayer coupling and defines a transformed BdGHamiltonianh˜k =ξk ∆˜k1 ge−iφ/2 0∆˜k1 −ξk 0 −geiφ/2geiφ/2 0 ξk ∆˜k20 −ge−iφ/2 ∆˜k2 −ξk . (4.10)Here ∆˜ka are defined as in Eqs. (4.8) but with the phase factors e±iφ/2 omitted. Unlessotherwise noted, we use this representation of the BdG Hamiltonian, and, for simplicity,we drop the tilde sign.The gap equation follows from ∂FBdG/∂∆ = 0 and reads∆ = 2V∑kα∂Ekα∂∆tanh12βEkα, (4.11)where β = 1/kBT is the inverse temperature. By noting that Ekα = ⟨kα|hk|kα⟩, where|kα⟩ is an eigenstate of hk, Eq. (4.11) can be rewritten in a form that is more suitable fornumerical evaluation,∆ = 2V∑kα⟨kα|∂hk∂∆|kα⟩ tanh 12βEkα. (4.12)Here ∂hk/∂∆ is a fixed 4× 4 matrix that follows from Eq. (4.10) and the required matrixelement is easily evaluated from the knowledge of the eigenstates. For any chosen phase φand temperature T the gap equation is then solved by iteration starting from a suitableguess for ∆.Similarly, by differentiating the free energy with respect to φ, it is possible to derive aconvenient expression for the interlayer supercurrentJ(φ) = −2eℏ∑kα⟨kα|∂hk∂φ|kα⟩ tanh 12βEkα. (4.13)234.3. Temperature dependence of critical currentIt is worth noting that this expression is non-perturbative and valid to all orders in theinterlayer coupling g. This is in contrast to many previous works on twisted junctions[90, 104, 105] which typically include only the leading ∼ g2 term. Retaining higher-ordercontributions is crucial for understanding the physics near θ = 45◦. For example, thedouble-pair tunneling process that gives the leading non-vanishing contribution in thisregime comes in at order g4.Henceforth, we will denote all quantities in units where 2e/ℏ = 1. And for purposes ofnumerical evaluation, the momentum sums in Eqs. (4.11-4.13) are converted into integralsusing the standard BCS procedure, which assumes a constant density of states NF at theFermi energy and introduces an energy cutoff ϵc to account for the states surrounding theFermi level that form the condensate. These details are summarized in Appendix C.3.4.3 Temperature dependence of critical currentCurrent-phase relationIn cuprates, the layers are very weakly coupled. Experimentally, it has now been establishedthat the superconducting transition temperature of a monolayer is identical to that of abulk crystal [72]. Correspondingly, the interlayer coupling is a weak perturbation on theSC order parameter amplitude and indeed we find that self-consistently determined ∆ isessentially independent of g. In other words, at a given temperature, the relevant degree offreedom is the phase difference φ and, without loss of generality, one can fix the amplitude.Typical current-phase relations obtained for parameters relevant to Bi2212 are displayedin Fig. 4.2. At zero twist angle, not surprisingly, we obtain a sinusoidal looking J(φ) at alltemperatures. A fit to Eq. (4.3) reveals that a small Jc2 term is required to capture thedata at low T . In order to achieve a close fit we add a third harmonic, −Jc3 sin(3φ). Therequired Jc3 is generally two orders of magnitude smaller than Jc1,2 and has no appreciableeffect on the physics. At twist angles close to 45◦, see Fig. 4.2(b), a significant secondharmonic term is required to capture the form of J(φ) obtained from the microscopicmodel. The negative slope of J(φ) at the origin observed for the low T curves signals thespontaneously T -broken phase: φ = 0 corresponds to a local maximum of F(φ) while twominima occur at a non-zero phase ±φmin.Fig. 4.2(c) shows an example of the dependence of parameters Jcα with α = 1, 2, 3 ontemperature, typical for twist angles close to 45◦. We observe that at low T the secondharmonic Jc2 dominates over the regular Josephson tunneling Jc1. The third harmonicis always very small and can be safely neglected when making qualitative arguments.By differentiating Eq. (4.3) it is easy to see that the condition ρs = 0 is attained whenJc1 = 2Jc2, which implies that the system is in the T -broken phase when Jc1 < 2Jc2. Thisis the criterion we used to establish the phase diagram shown in Fig. 4.1.Critical current and superfluid densityHere we summarize our results for the T -dependent interlayer critical current Jc(T ) andsuperfluid stiffness ρs(T ) defined in Eq. (4.4). At intermediate twist angles, Jc in Fig.4.3(a) exhibits an increase with increasing temperature. Such behavior is anomalous –244.3. Temperature dependence of critical current0 0.5 1T/Tc0.000.050.100.150.20Jc(T)0 0.5 1T/Tc0.000.050.100.150.20ρs(T)ec=60, v1=0.12, v2=0.00, v0=0.00, g=7.Jc( )0 0.5 1T/Tc0.000.050.100.150.20ρs(T)ec=60, v1=0.12, v2=0.00, v0=0.043, g=7a) c)b)θ = 0o5o10o15o25o35o20o30oFigure 4.3: (a) Critical current Jc and (b) superfluid stiffness ρs as a function of temperaturefor various twist angles. The traces for twist angles above 40◦ are shown in 1◦ increments.Red lines at θ = 10◦ correspond to the magic angle θM. Blue lines show results forθ = 42− 45◦ where time reversal is spontaneously broken at low T . Panel (c) shows theresults for the superfluid stiffness calculated in the presence of the secondary s-wave orderparameter that arises at low T close to θM. In these plots the same model parameters areused as in Fig. 4.2.normally thermal excitations cause Cooper pair breaking which depletes the SC condensateand results in a monotonic decrease of Jc as a function of temperature – and it reflectsthe sign-changing nature of the d-wave order parameter and provides an experimentallyaccessible signature of unconventional superconductivity in cuprates. Note that ρs showssimilar behavior which is even more pronounced, especially near the magic angle θM ≃ 15◦.The anomalous increase in ρs(T ) and Jc(T ) can be understood as follows. We focus onρs(T ) but similar arguments apply to Jc(T ). As shown in Appendix C.1, the superfluidstiffness can be expressed asρs(T ) =∑k∆k1∆k2Ω(ξk, T ), (4.14)where Ω(ξk, T ) ≥ 0. The key observation is that the product of the two gap functions∆k1∆k2 = ∆2 cos(2αk − θ) cos(2αk + θ) becomes negative for angles αk close to ±π/4when the twist angle is non-zero. At T = 0 these nodal regions of k-space, therefore, makea negative contribution to ρs. This negative contribution grows with increasing twist whichaccounts for decreasing ρs(0) observed in Fig. 4.3(b). Thermal excitations cause Cooperpair breaking and at low T this happens predominantly in the nodal regions of k-spacewhere the low-lying excitations reside. Since Cooper pairs in the nodal regions make anegative contribution to ρs, removing these leads to an anomalous increase of ρs(T ) withincreasing temperature.The T -broken phase is signalled by ρs turning negative. For our parameters, thisoccurs when θ > 42◦ and corresponds to blue curves in Fig. 4.3. Unfortunately, whilestraightforward to calculate, ρs is not easily measurable. Critical current, on the other hand,is experimentally accessible but does not show any clear signatures of the spontaneously254.4. Estimate of interlayer couplingT -broken phase. Physically, Jc represents the largest current that the system can sustainbefore going normal. At such a large externally imposed current time-reversal symmetry isexplicitly and strongly broken, so Jc will not be sensitive to a small amount of T -breakingalready present in the unbiased system.4.4 Estimate of interlayer couplingThe interlayer coupling g is a key parameter whose strength determines the position of themagic angle θM as well as the width of the T -broken topological phase near the 45◦ twist.In the literature on high-Tc cuprates, one can find many estimates of g, but they vary widelyfor different materials. For Bi2212 one finds values (obtained from both experimentalfits and theoretical modelling) ranging between 5− 120 meV depending on the techniqueused [113]. This wide range can be attributed to the complicated crystal and electronicstructure of Bi2212, as well as the fact that g is likely momentum-dependent and differentprobes are sensitive to different parts of the Brillouin zone.To provide quantitative predictions for transport in twisted Bi2212 bilayers, we nowattempt to extract the likely value of g by comparing our theoretical models to the recentcritical current measurements in twisted Bi2212 flakes [50]. Given various uncertainties,we give two separate estimates: one based on the GL theory valid at high temperaturesand one based on microscopic modeling relevant in the T → 0 limit.The experimental measurement (Fig. 2D in Ref. [50]) gives temperature dependenceof JcRN for various twist angles. Here RN is the normal-state junction resistance whichserves to normalize Jc for different junctions, effectively providing a measure of the criticalcurrent density that is independent of the junction area and geometry. Our GL estimaterelies on extracting the critical angle θc from Eq. (4.5) based on the data and then backingout the value of g by matching to the microscopic phase diagram shown in Fig. 4.1. Wechose T = 50 K as our reference temperature and read off JcRN = (10.4, 8.2, 3.5, 0.2) mV atθ = (0◦, 29◦, 39◦, 44.9◦). We assume that for twist angles θ < 40◦ the current-phase relation(4.6) is dominated by the first term and we can approximate Jc(θ) ≃ Jc1(θ) = Jc1(0) cos(2θ).Close to 45◦ the second term dominates and we take Jc2 ≃ Jc(44.9◦). Eq. (4.5) then givesθc = (44.1± 0.2)◦. Finally, comparing to the microscopic phase diagram in Fig. 4.1 we seethat this range of critical angles at T = 50 K is well captured by taking g ≃ 10.5 meV.This corresponds to θc ≃ 39◦ and θM ≃ 15◦.A more direct but technically somewhat more challenging estimate of g can be givenstarting from the microscopic expression for J(φ) given in Eq. (C.2) of Appendix C.1.Because at θ = 0 the current is dominated by the ordinary Josephson tunneling whilenear 45◦ double Cooper pair tunneling dominates, one expects the leading term in theT → 0 limit to behave in the two cases as ∼ g2 and ∼ g4, respectively. Indeed as shown inAppendix C.2 one can estimate, at T = 0,Jc(θ = 0◦) ≃ 2πNF eg2ℏC2 (4.15)Jc(θ = 45◦) ≃ 2πNF eg4ℏ∆2C4. (4.16)264.5. Interlayer coupling dependence of critical current0.0 0.2 0.4 0.6 0.8 1.0 1.2T/Tc0.000.250.500.751.001.25Jc£10°312.7±16.3±22.6±28.1±36.9±46.4±0.0 0.2 0.4 0.6 0.8 1.0 1.2T/Tc02468Jc£10°312.7±16.3±22.6±28.1±36.9±46.4±a) b)Figure 4.4: Behavior of the critical current for various accessible commensurate twistangles in the lattice model (details in Appendix B.2). The parameters used are µ = −1.3t,g0 = 10 meV in (a) and µ = −1.4t, g0 = 30 meV in (b).Here C2 and C4 are dimensionless constants of order one which depend on model detailssuch as the Fermi surface shape and the pairing interaction cutoff scale ϵc. For simplicityand concreteness we assume a circular Fermi surface and work in the limit ∆≪ ϵc. Underthese conditions C2 = 1 and C4 ≃ 0.55. Taking a ratio one can determine g asg ≃ ∆√C2C4Jc(θ = 45◦)Jc(θ = 0◦). (4.17)From the experimental data (Fig. 2D in Ref. [50]) we estimate JcRN (θ = 45◦)/JcRN (θ =0◦) ≃ 0.010, which together with ∆ = 45 meV gives g ≃ 6.2 meV.Given various simplifying assumptions, the two estimates agree reasonably well. Inview of the values found in the literature (typically 10s of meV), we conclude that g around10 meV is an appropriate value to use when modelling Josephson effects in twisted Bi2212.4.5 Interlayer coupling dependence of critical currentResults reported in Fig. 4.3 were obtained using the simple continuum model for a d-waveSC. Because the interesting low-T behavior arises from the physics of nodal excitations,we expect this model to provide an accurate description of twisted bilayer cuprates. Asa further check, we compute critical currents from the more realistic lattice model thataccurately captures the hole-like Fermi surface of Bi2212 near optimal doping. Herethe interlayer tunneling, while still momentum conserving, has a soft cut-off and decaysexponentially in the distance between the sites. Fig. 4.4 shows the temperature dependenceof Jc for a range of commensurate angles θm,n = 2arctan(m/n) with integers m,n. Forthese angles, the twisted bilayer forms a moire´ pattern with 2(m2 + n2) sites per unit cell.As in the continuum model, an anomalous increase in Jc at low temperatures forintermediate twists is observed in Fig. 4.4(a). Such behavior is observed only if the interlayertunneling is small. Specifically, as shown in Fig. 4.4(b), when g is increased the curvesresemble the familiar Ambegaokar-Baratoff form [114] where the current monotonicallydecreases with T . Note that this dependence on tunneling strength is a feature also ofthe continuum model discussed in the previous sections. We attribute this change in274.6. Magic angles and secondary order parametersbehavior to non-perturbative effects: when g becomes comparable to ∆ simple perturbativearguments underlying our reasoning in Sec. 4.3 no longer apply. Nevertheless, the observedagreement in the characteristic behavior of the two models lends support to our predictionsfor Jc(T ) in twisted cuprates.4.6 Magic angles and secondary order parametersAn interesting situation occurs at small twist angles when Dirac points originating from thetwo layers merge pairwise and give rise to a quadratic band crossing (QBC) in each quadrantof the Brillouin zone. In analogy to the flat-band formation in twisted bilayer graphene,the twist at which that QBC occurs has been called the magic angle. As discussed in Ref.[48], such a merger of Dirac cones is unavoidable as long as time reversal and two-foldrotation (along the in-plane diagonal) are respected and it occurs when θ = θM ∼ 2g/∆.Similar to bilayer graphene, the system becomes susceptible to interactions for twistangles close to θM. In a Dirac dispersion, the density of states (DOS) vanishes at zero energyand interactions are perturbatively irrelevant. At a QBC, however, the DOS becomesconstant at low energies and the system behaves like a metal of Bogoliubov quasiparticles.Such a metal will be unstable with respect to residual interactions that have not beenincluded in the original mean-field treatment of d-wave superconductivity. For instance, asdiscussed in Ref. [48], in the presence of attractive interactions in s or dxy channels, thesystem will develop a secondary SC order parameter for twist angles close to θM.To illustrate this behavior we focus on a secondary instability in the s-wave channel andinvestigate its effect on the behavior of Jc(T ) and ρs(T ). Working in the gauge describedin Eq. (4.10) the pair amplitudes can be expressed as∆k1 = ∆cos(2αk − θ) + i∆s, (4.18)∆k2 = ∆cos(2αk + θ) + i∆s, (4.19)where ∆ and ∆s are real. The imaginary unity accompanying ∆s is required to break T ;the QBC is protected by T and a gap can only open when time reversal is broken. Thegap equation for the secondary order parameter follows from minimizing the free energyEq. (4.9) with respect to ∆s and reads∆s = 2Vs∑kα⟨kα| ∂hk∂∆s|kα⟩ tanh 12βEkα, (4.20)where Vs denotes interaction strength in the s channel.Fig. 4.1 shows the typical phase diagram obtained by numerically solving the coupledgap equations (4.12) and (4.20). We observe a broad dome of s order parameter (peakedat θ = θM), in addition to the previously discussed d+ id′ phase. It is to be noted thatparameters in Fig. 4.1 are chosen to maximize the secondary instability; if Vs were chosenany larger we would find nonzero ∆s even at zero twist angle, contrary to experimentalobservations which indicate pure dx2−y2 order parameter in untwisted crystals and films.These considerations suggest that secondary order parameter physics would be visible onlyat low temperatures, as compared to the native Tc of cuprates. Similar results are obtainedfor a secondary instability in the dxy channel.284.7. Fraunhofer patternJc(T ) and ρs(T ) can be calculated as before even in the presence of the secondaryorder parameter. We find that ρs(T ) shows a clear signature of the secondary order, asillustrated in Fig. 4.3(c). The enhancement at low T can be attributed to the contributionof s-wave Cooper pairs to the tunneling between the layers. Jc(T ), on the other hand, isunaffected by the secondary order parameter. This is because critical current is achievedwhen the phase difference φ is close to π/2; we find that in this regime Eqs. (4.12) and(4.20) produce a self-consistent solution with ∆s strongly suppressed compared to its valueat zero phase. Thus, one more time we find that the critical current alone is not a suitableprobe for establishing phase transitions.4.7 Fraunhofer patternFraunhofer interference occurs when a Josephson junction is subjected to an externalmagnetic field parallel to the junction plane [115, 116]. To describe the effect of the fieldon inter-plane Josephson phase and current distribution we must allow for spatial variationof the order parameters ψa in the plane. The starting point of the calculation is Eq. (4.1),except that the interlayer Cooper pair hopping terms now contain magnetic phases pickedup when tunneling between the layers. We thus consider a free energy densityf [ψ1, ψ2] = f0[ψ1] + f0[ψ2] +A|ψ1|2|ψ2|2 + (Bψ1ψ∗2e−iqx + Cψ21ψ∗22 e−2iqx + c.c.) (4.21)where q = Bd/Φ0 encodes the magnetic flux per unit length that is threaded between thelayers separated by distance d, andf0[ψ] = α|ψ|2 + 12β|ψ|4 + γ|∇ψ|2 (4.22)describes each layer. We choose a gauge ψ2 = ψei(φx−qx) in Eq. (4.21), where φx is thespatially dependent Josephson phase difference and ψ is assumed constant and real. Thisleads to a free energy densityf(φx) = f0 + γ˜(∂xφx − q)2 − ℏ2eV[Jc1 cosφx − Jc22cos(2φx)], (4.23)where V denotes the 2D volume of the system, γ˜ = γψ2, and we omitted higher harmoniccontributions (Jc3) that are found to be very small, see Fig. 4.2(c).Notice that in addition to the parameters Jcα, one also requires the in-plane superfluidstiffness, related to the parameter γ in the free energy Eq. (4.22). The calculation of thisquantity for a single-layer dSC is standard in the high-Tc literature [117–120], but for thesake of completeness we include it for our specific model in Appendix C.3. Using this input,one can obtain results for the phase evolution and winding of the SC phase difference φxas a function of both the applied in-plane magnetic flux and temperature.For a set of zero-field model parameters γ and Jc1,2, it is easy to discretize andnumerically determine solutions φx that minimize the above free energy for a given totalflux Φ/Φ0 = qW in a twist junction of width W (measured in the direction perpendicularto the in-plane field). For small to moderate flux densities, we find that the local Josephson294.7. Fraunhofer patternFigure 4.5: Fraunhofer patterns. Critical current Ic(Φ, T ) vs applied magnetic flux Φin the twist junction for various temperatures T , with twist angle θ = 42◦ (left) andθ = 44◦ (right). The curves are computed from a Lawrence-Doniach type model of thebilayer junction, taking into account the first and second harmonic terms of the current inEq. (4.24) with coefficients Jc1,c2 extracted from the continuum model.phase increases approximately linearly, φx ≈ φ0 + qx, as expected in short Josephsonjunctions [112, 121]. This trivially minimizes the kinetic term and leaves one with the(local) Josephson current-phase relation,jx(φx) = jc1 sin(φ0 + qx)− jc2 sin(2φ0 + 2qx), (4.24)where jca = Jca/V denote the respective densities. We then calculate the critical currentby integrating the above Josephson current density and maximizing with respect to thephase offset φ0, as described by Goldobin et al [112]. The final expression for the totalcritical current (normalized by V ) then takes the formJc(θ,Φ, T ) = Ic1∣∣∣∣sin(πΦ/Φ0)πΦ/Φ0∣∣∣∣+ Ic2 ∣∣∣∣sin(2πΦ/Φ0)2πΦ/Φ0∣∣∣∣ , (4.25)where Icα = Icα(θ, T ) ∼ Jcα(θ, T ). As already alluded to, the Fraunhofer pattern thusinterpolates between the regular and period-halved cases, as a function of both twist angleand temperature.We show data for twist angles θ = 42◦ and 44◦ and with varying temperature T inFig. 4.5. Since the d + id′ phase is prominent only for a range of twist angles aroundθ = 45◦, the period-halved contribution in the Fraunhofer pattern is rapidly lost whenmoving to θ ≲ 40◦. In this case the pattern comprises slightly deformed regular periodiclobes which might be difficult to identify experimentally. However, as θ → 45◦ the period-halved contribution becomes dominant. With the presently accessible twist angle control,assuming that samples have a fairly uniform twist angle distribution, it should thus befeasible to observe period-halved Fraunhofer patterns.To estimate the magnetic field that is required to thread a single flux quantumΦ0 ≈ 2.067 · 10−3 Tµm2 through the junction, we take the c-axis lattice constant ofBSCCO d ≈ 3.1 nm as the effective junction thickness. For a junction of widthW ≈ 10 µm,cf. Refs. [50, 72], the magnetic field then comes out to be B = Φ0/(dW ) ≈ 67 mT. Suchweak parallel magnetic fields are perfectly compatible with maintaining superconductivity304.8. Shapiro steps°2 °1 0 1 2Jdc°2°1012V(h¯! 2e)µ = 45o, Jrf = 0.70! = 0.6! = 0.9! = 1.20 1 2 3Jrf°2°1012V(h¯! 2e)µ = 45o0120 1 2 3Jrf°2°1012V(h¯! 2e)µ = 15o012a) b) c)Figure 4.6: Shapiro steps. (a) Current-Voltage curves in the 45◦ twisted configuration inthe presence of an external rf drive with amplitude Jrf = 0.7 at three different frequencies.The voltage is scaled in units of ℏω/2e to highlight the Shapiro steps and the lowest fourfractional steps are indicated with the horizontal dashed lines. (b-c) Voltage as a functionof drive current for two twists that are in the topological and trivial phases respectively,with ω = 0.6. The color bar indicates step amplitudes. Since the numerical values of Jc1,2at different twists vary widely, to aid comparison, J(φ) is normalized such that the criticalcurrent in the absence of rf drive (Jrf = 0) is unity (or equivalently the zeroth order stephas width 2). Clearly, fractional steps appear only for twist angles close to 45◦, as in panel(b). At high rf currents, one enters a regime where the step amplitudes oscillate. Thetemperature is fixed at T = 0.5Tc and R = 0.7 in all simulations.in BSCCO flakes. We also note that experimental devices usually consist of multipleBSCCO layers per flake, and the effective junction thickness d in twisted junctions is foundto be in the range of tens of nm. This further eases magnetic field requirements.Finally, we note that while the precise magnitudes and dependencies of parametersIc1,2(θ, T ) in Eq. (4.25) depend on the assumptions made in this section (e.g. linear growthof φx), the functional form with first and second harmonics in the Fraunhofer pattern isgeneric.4.8 Shapiro stepsTopological superconductivity with a dominant second harmonic in the current-phaserelation can also be made manifest by subjecting the system to an external drive. If aJosephson junction is irradiated with a radio-frequency source, an AC voltage is inducedthat frequency modulates the AC Josephson current. When the dynamics of the junctionis phase-locked to the rf drive, the supercurrent shows constant voltage Shapiro steps atvoltages Vn = nℏω/2e, where ω is radiation frequency and n ∈ Z is the step index [122]. Ina π periodic junction, additional steps appear at fractional values Vn/2, which are a directexpression of a dominant second harmonic in the current-phase relation.In a usual experiment, the time-averaged voltage is measured in response to a currentbias. To study the dynamics of such a current-driven system, it is convenient to work withina semiclassical framework where current is carried in three parallel channels: supercurrent314.8. Shapiro stepscarried by the Cooper pairs, a resistive path for dissipative current and a capacitive channelthat accounts for charge build up on the superconducting leads. If the junction has anegligible geometric capacitance, one is in an over-damped regime where the resistiveshunt provides the only impedance. In this so-called resistively shunted junction (RSJ)model [116], the total current through the system can be written as J = JR + J(φ), whereJR = V/R is the current through the resistor R and J(φ) is the current-phase relation(4.3). Finally, utilizing the universally valid superconducting phase evolution relationdφ/dt = 2eV/ℏ and imposing the current drive J = Jdc + Jrf sin(ωt), one obtains the firstorder differential equationℏ2eR∂φ∂t+ J(φ) = Jdc + Jrf sin(ωt). (4.26)For a given twist, temperature, rf-drive parameters, R and Jdc, the time evolution ofthe phase and, via a time derivative, the voltage can be obtained by solving the aboveequation numerically through a routine Runge-Kutta algorithm. The representative time-averaged voltage behavior as a function of the direct current is depicted in Fig. 4.6(a),where fractional steps corresponding to the halved period in the T -broken phase are seen.Another way to visualize Shapiro physics is to study the dependence of the step widthsas a function of rf current, where the steps are revealed as maxima at quantized voltages,see Fig. 4.6(b-c). Since Prf = J2rfRrf , Jrf also serves as a proxy for rf power. Steps appearprogressively, starting with low values of n, as the drive amplitude is increased and athigher powers one observes an oscillatory pattern [115, 123].The constant voltage steps may be intuitively understood as the tunneling of Cooperpairs across the barrier that is accompanied by an exchange of photons with the radiation.Whenever the potential energy across the junction is equal to the photon energy ℏω, or amultiple thereof, a Cooper pair can absorb (emit) photons from (to) the radiation field.The nth Shapiro step corresponds to n photons being exchanged. The fractional steps, onthe other hand, coincide with two Cooper pairs tunnelling across the junction.With all else fixed, step amplitudes are proportional to the critical currents Jc1,2 thatare informed by the microscopic model. To observe all steps and avoid interference, onehas to ensure that they are sufficiently apart on the Jdc axis; the separation is controlled bythe amplitude and frequency of the drive. While Jc2 is non-zero at any point in the phasediagram, we find empirically that having Jc2/Jc1 > 2 is required to discern the fractionalsteps. This is a stronger condition than Jc2/Jc1 > 1/2 which, according to our discussionin Sec. 4.3, defines the T -broken phase. Therefore, appearance of fractional Shapiro stepsin experimental data can be regarded as a strong sign of the T -broken phase in the system.Temperature implicitly enters our model through microscopic parameters Jc1,2. To fullyaccount for the thermal effects, one needs to add a noise term to the bias current [124].The main consequence of such a treatment is the ‘rounding’ of steps in the current-voltagecharacteristic: In the presence of fluctuations, the switch from a constant voltage plateauto a dissipative state would no longer be sharp, resulting in smaller step amplitudes.Nevertheless, noise does not alter the location of the steps and the qualitative featuresdiscussed above continue to hold.It is to be noted that fractional Shapiro steps can arise in physical systems that areremoved from the context of topological superconductivity. For instance, two-dimensional324.9. Thick flakes0 1 2 3ϕ-0.0100.01J(ϕ)45o44o0 0.5 1T/Tc0.000.050.100.150.200.25Jc(T)M=1M=2M=40 0.5 1T/Tc0.000.050.100.150.200.25ρs(T)ec=60, v1=0.12, v2=0.0, v0=0.0, g=10.5a) M=12 b) c) 0o10o20o30o40o45oFigure 4.7: (a) Current-phase relation computed for two flakes, each with M = 12monolayers, for twist angles near 45◦. Various curves are for temperatures T betweenzero and Tc in increments Tc/20. (b) Critical current Jc(T ) and (c) phase stiffness ρs(T )extracted from the calculated J(φ). We obtained results for M up to 24 but we find thatthe curves become essentially independent of M when M ≥ 4.Josephson junction arrays may show half-integer steps due to a skewed current-phase rela-tion [125, 126]. A sin(2φ) dependence of the current also appears in magnetic Josephsonjunctions [127]. In the present setting, given our theoretical understanding, topological su-perconductivity should be the primary candidate for the physics underlying unconventionalShapiro physics.4.9 Thick flakesWe briefly consider twisted structures composed of thicker flakes that might be easierto assemble and probe in the lab [49, 50]. Specifically, we study the Josephson currentbetween two flakes each composed of M cuprate monolayers. The microscopic Hamiltoniandescribing this situation can be constructed as a straightforward extension of Eq. (4.7).For example the M = 2 system is represented by an 8× 8 matrix BdG HamiltonianHk =hk1 Γ 0 0Γ hk1 Γ 00 Γ hk2 Γ0 0 Γ hk2 , (4.27)wherehka =(ξk ∆ka∆∗ka −ξk), Γ =(g 00 −g), (4.28)describe the individual monolayers and their coupling, respectively. The twist betweenthe flakes is encoded in the structure of the order parameter given in Eqs. (4.8). For ageneric M the Hamiltonian becomes a 4M × 4M matrix with a structure that follows asan obvious generalization of Eq. (4.27).The supercurrent flowing between the two flakes can be calculated from Eq. (4.13). Wefind that many of the features observed previously for two coupled monolayers (the M = 1334.10. Discussioncase) persist in thicker flakes. Fig. 4.7 shows some representative results. In panel (a)we observe that the characteristic crossover to π-periodic oscillations in J(φ) as the twistapproaches 45◦ persists for M > 1. The critical current Jc and phase stiffness ρs likewiseshow similar characteristic temperature dependence forM > 1 as monolayer devices, exceptfor a slight suppression of the amplitude at low T . Interestingly, the largest change occursbetween M = 1 and M = 2 cases; for M ≥ 2 the system appears to have reached the bulklimit and we see no perceptible change in Jc(T ) and ρs(T ) as M is increased further.A quantity that strongly depends on thickness M is the spectral gap Egap which wedefine as the energy of the lowest quasiparticle excitation above the ground state, obtainedas the smallest positive eigenvalue of Hk. The spectral gap is measurable through variousspectroscopic probes, such as electron tunneling or angle-resolved photoemission. Fig. 4.8depicts the spectral gap as a function of twist angle for several values of M . The gap islargest at θ = 45◦, where the T -breaking induced by the interlayer tunneling is expectedto be maximal. The maximum gap Egap(45◦) is seen to rapidly decay with increasing M ;the inset suggests an exponential dependence on M . Indeed this is to be expected: in thelimit of thick flakes the tunneling becomes a surface perturbation to a 3D system thatcannot open a bulk gap. Nevertheless, the T broken phase persists in this limit as can beseen from the behavior of the critical twist angle θc also shown in the inset. θc is definedas the twist angle at which the free energy FBdG first develops minima away from φ = 0,signalling the onset of T breaking. The inset shows that with the increasing flake thicknessθc quickly approaches saturation, leading to a stable range of twist angles |θ − 45◦| < θcin which the system ground state exhibits spontaneously broken time-reversal symmetry.In this regime, technically, the system is in the topological phase with non-zero Chernnumber and protected chiral edge modes. However, for larger M this topological phase isprotected by a gap that becomes exponentially small, thus limiting its potential usefulness.We may conclude that even though the T -breaking gap quickly becomes small as thethickness of the flakes increases, the current-phase relation continues to display signaturesof the spontaneous T breaking in this limit. Intuitively, this behavior can be understood bynoting that the supercurrent is controlled by interface properties and critical currents areonly weakly affected by conditions at the surface. We thus expect Fraunhofer interferencepatterns and fractional Shapiro steps discussed in the context of monolayer-thin flakes toremain a useful probe of spontaneous T breaking in the limit of thicker flakes. Spectro-scopies, on the other hand, will be most useful when applied in the limit of monolayer-thinflakes.4.10 DiscussionTheory suggests that twisted cuprate bilayers can form topological phases with supercon-ducting critical temperatures that are close to the bulk critical temperature of the parentcompounds. These exciting predictions beg the question: what would be an unequivocalsignature of time reversal symmetry breaking and of a high-temperature topological super-conducting phase? The full answer likely requires complementary evidence from differentexperimental quarters. Yet, quantum transport through the bilayer might be a naturalstarting point.344.10. Discussion0 15 30 45twist angle [deg]0.00.51.01.52.02.5E gap [meV]M=1M=2M=40 4 8 12M0.010.11E gap(45o )30354045θ c[deg]Figure 4.8: Evolution of the spectral Egap with the flake thickness M . Note that Egapvalues below ∼ 0.1 meV are comparable to the finite-size gap in our system (500 × 500grid of k points) and are therefore consistent with gapless behavior. Inset shows scaling ofthe maximum gap Egap(45◦) with M , suggesting exponential decay, and the critical twistangle θc beyond which the T -broken phase sets in.Because each layer is a standalone superconductor, the twisted structure forms aJosephson junction, provided that the interlayer coupling is weak. At the fundamentallevel, T breaking becomes apparent in the current phase relation J(φ) which reflects thedouble-minimum structure of the Josephson free energy F(φ). But when the experimentallymeasurable critical current is analyzed as a function of temperature and twist, somewhatsurprisingly, we find the behavior to be smooth. In other words, a transition to thetopological regime cannot be earmarked by simply looking at the critical current of anotherwise unperturbed junction.When the bilayer is subjected to certain external perturbations, however, the higherharmonic term corresponding to double Cooper pair tunneling can become manifest. Underthe influence of an in-plane magnetic field, the junction shows a Fraunhofer oscillationpattern consistent with the new π-periodic term. Similarly, when the junction is drivenwith electromagnetic radiation, fractional Shapiro steps appear in the current-voltagecharacteristic.It is to be noted that, strictly speaking, detection of the π-periodic term is by itself notsufficient to identify the T -broken topological phase in our twisted system. An additionalrequirement is that the second harmonic term Jc2 comes with the opposite sign relativeto the fundamental term Jc1, as indicated in Eq. (4.3) and that 2Jc2 > Jc1, which is thecondition for the onset of the T -broken phase. Our theoretical considerations show thatnatural models of twisted bilayers produce the correct signs of Jc1,2. In addition, we findthat clearly observable deviations from conventional Fraunhofer and Shapiro responsesoccur only when the condition 2Jc2 > Jc1 is well satisfied. This suggests that such effectscan be taken as strong signatures of the topological phase.For the sake of simplicity and to keep the discussion self-contained we neglected theeffect of random disorder throughout this work. It is however known that disorder can354.10. Discussionplay an important role in the physics of cuprates, especially for properties that dependsensitively on interlayer tunneling. Experimental evidence from the c-axis transport in(untwisted) cuprate crystals suggests that ‘impurity assisted’ hopping could be an importantchannel for interlayer transport [128]. The underlying physics is captured by a familyof incoherent tunneling models [129–132] which generally yield better agreement withthe data than models where momentum is fully conserved. One might thus expect thatdisorder effects will play a role in twisted bilayers and may in some cases even dominate theinterlayer transport. Our investigation in Ref. [133] shows that the basic phenomenologyobtained on the basis of simple clean-limit models discussed here, including the spontaneousT -breaking near 45◦ twist and period doubled current-phase relationship, remains in placefor incoherent tunneling models, at least when disorder varies slowly on the lattice scale.Understanding fully the effect of disorder on twisted cuprate bilayers will be importantwhen making detailed comparisons to experimental data.The current-phase relation describes supercurrent through the 2D bulk of the constituenttwisted superconductors. In this sense, Fraunhofer patterns and Shapiro steps are amaterialization of bulk physics. Chiral edge modes, mandated by the bulk boundarycorrespondence, and finite orbital magnetism of the quasiparticles are others features ofthe topological phase. Spontaneous magnetic fields associated with the edge currents[68] could potentially be measured through local probes such as scanning SQUID. Directexamination of the edges through transport or spectroscopic measurements would provideuseful complementary insights into the problem. From the viewpoint of experiments,however, this is arguably a more difficult task.36Chapter 5Twisted multilayer nodalsuperconductors5.1 IntroductionWithin the framework of Bardeen–Cooper–Schrieffer (BCS) mean-field theory, we haveshown that twisted bilayers of 2D nodal d-wave superconductors are a candidate platform torealize topological superconductivity. However, it stands to reason that the phenomenologyof twist-induced spontaneous T -breaking is not specific to cuprates, but applies to nodalsuperconductors in general. Advances in exfoliation and material growth techniquescould potentially reveal new nodal superconductors in organics [134], quasi-2D limits ofiron-based [135] or heavy-fermion [136] materials, or in transition metal dichalcogenidesfollowing recent theoretical proposals [57, 137, 138]. There is thus hope that 2D nodalsuperconductors could form the basis for various heterostructures with novel properties [92],possibly opening new routes towards Majorana zero-modes [139].Here we explore the physics of twisted structures comprising a generic number L ofnodal superconducting layers. Going beyond bilayers affords two key advantages: (i) agreater variety of possible phases with additional tunability of physical properties, and (ii)an effective renormalization of parameters that stabilizes topological superconductivityfor weaker interlayer tunneling as well as smaller twist angles far from 45◦, potentiallyfacilitating experimental realizations.As the twist degree of freedom of each layer can be tuned independently, multilayersadmit several inequivalent arrangements. In this chapter, we consider two high-symmetryconfigurations that arise by fixing a unique twist angle θ between neighboring layers: (i)alternating twists (AT), where the relative twists between successive layers differ in signand (ii) chiral twists (CT), where neighboring layers are rotated by the same angle suchthat a net chirality may be associated with the system.To shed light on phenomenology, we forgo materials-specific considerations and the-oretical complications stemming from the treatment of strong interactions and adopt aBCS mean-field approach as a first exposition of this problem. We work mainly within aself-consistent continuum formulation, where each layer is treated as a 2D nodal supercon-ductor with dx2−y2 symmetry coupled to neighboring layers via electron tunneling. Wethen connect our results to a twisted multilayer lattice model, which gives us access totopological indices and edge modes. Contrasting the phase diagrams of the two stackings,we obtain the following key results.The physics of AT multilayers can be understood via a mapping to multiple copies ofthe bilayer problem [38]. T -breaking superconductivity occurs for twist angles θ close to375.2. Model45◦, but the spectrum is fully gapped only when the number of layers L is even. This isbecause for odd L the system comprises an effectively decoupled monolayer sector with anodal d-wave order parameter. In analogy to twisted bilayers, quadratic band touchingsoccur at a set of ‘magic angles’. We find that such band touchings are susceptible toT -breaking secondary instabilities through two distinct mechanisms. First, the interlayertunneling between twisted layers tends to nucleate a d + id′ order [47] – an effect thatis enhanced in multilayers due to a renormalization of the tunneling strength. Second,residual interactions left out of the mean-field BCS treatment become marginal at theQBTs and can nucleate different symmetry-breaking orders [48].CT multilayers spontaneously break both time-reversal T , but preserve the product C2T ,where C2 is a π rotation about the diagonals of the middle layer. Their superconductingstate is fully gapped and topological, with chiral Majorana modes propagating on sampleboundaries, for generic large angles different from 45◦. For small twist angles, the Diraccones in the BdG quasiparticle dispersion merge to form either quadratic or cubic bandtouchings, depending on L. The latter has an enhanced susceptibility to gap opening dueto a power-law divergent density of states at zero energy, as opposed to a constant densityof states in the quadratic case.The rest of this chapter is organized as follows. After describing the model in Sec. 5.2,we analyze the physics of twisted trilayers in Sec. 5.3 (AT stacking) and Sec. 5.4 (CTstacking). Sec. 5.5 presents a generalization to multilayers with L > 3. Finally, in Sec. 5.6,we comment on experimental signatures of the phases identified in this work and speculateon their relevance to ongoing investigations in twisted cuprates.5.2 ModelWe consider a continuum model consisting of nodal d-wave superconducting layers thatare twisted relative to each other and coupled through single electron tunneling. In a BCSmean-field description this can be written asH =∑kσlξkc†kσlckσl +∑kσ⟨l,m⟩glmk c†kσlckσm +∑kl∆kl(c†k↑lc†−k↓l + h.c.)−∑kl∆kl⟨c†k↑lc†−k↓l⟩.(5.1)The operator c†kσl creates an electron with momentum k and spin σ in layer l. We assumea rotationally invariant Fermi surface, with the in-plane kinetic energy ξk = k2/2me − µ(taking ℏ = 1 throughout), and glmk denotes the electron tunneling amplitude betweenneighboring layers l and m. The momentum sums span an energy interval ξk ∈ [−ϵc, ϵc]around the Fermi level, with ϵc a high-energy cutoff fixed at 60meV. We also considera lattice model with a Fermi surface more closely resembling that of high-Tc cupratesuperconductors and arrive at qualitatively similar results.The last two terms in Eq. (5.1) originate from a mean-field BCS decoupling of an intra-layer, zero pair-momentum attractive interaction term HI = 1N∑kpl Vlkpc†k↑lc†−k↓lc−p↓lcp↑lthat respects C4 rotation symmetry. Here N denotes the number of unit cells. We usea separable form V lkp = −2V cos(2αk − 2θl) cos(2αp − 2θl) where V > 0 is the interaction385.3. Trilayer: Alternating twiststrength, θl the rigid rotation angle of layer l and αk is the polar angle of the wavevector k.This choice greatly simplifies our analysis; we however expect our results to be representativeof nodal superconductors more generally. The gap function is then defined as∆kl =1N∑pV lkp⟨c−p↓lcp↑l⟩ = ∆l cos(2αk − 2θl), (5.2)where ∆l ≡ −2VN∑p cos(2αp − 2θl)⟨c−p↓lcp↑l⟩ is the complex order parameter of layerl. The form of the tunneling matrix elements glmk depends on the symmetries of thematerial under consideration, and in particular of the orbitals involved in various interlayertunneling processes [94, 140]. For simplicity we focus on a momentum-independent formglmk = g [47, 48] and comment on its relevance for cuprate multilayers towards the end.With these assumptions the Hamiltonian (5.1) can be written in the BdG formalism asH =∑kΨ†khkΨk + E0, where E0 =∑kl ξk + N2V ∑l |∆l|2. Trilayers are described by theNambu spinor ΨTk =(ck↑1, ck↑2, ck↑3, c†−k↓1, c†−k↓2, c†−k↓3)andhk =ξk g 0 ∆k1 0 0g ξk g 0 ∆k2 00 g ξk 0 0 ∆k3∆∗k1 0 0 −ξk −g 00 ∆∗k2 0 −g −ξk −g0 0 ∆∗k3 0 −g −ξk . (5.3)Analogous forms follow for multilayers with L > 3. The order parameters ∆l defined byEq. 5.2 are obtained self-consistently by minimizing the free energyF = E0 − 2β−1∑kαln[2 cosh(βEkα2)], (5.4)where β = 1/kBT is the inverse temperature and Ekα is the positive energy band αobtained by diagonalizing hk. Numerically, it is convenient to express the minimizationconditions ∂F/∂∆∗l = 0 as∆l =2VN∑kαtanh(βEkα2)⟨kα| ∂hk∂∆∗l|kα⟩, (5.5)where Ekα = ⟨kα|hk|kα⟩ has been used to recast the minimization conditions in terms ofthe matrices ∂hk/∂∆∗l with l = 1 . . . L.This continuum mean-field formulation is particularly attractive because it allows thetwist angle to be tuned continuously, in contrast to lattice formulations in twisted geometrieswhich in practice can only be studied at large commensurate angles (see Appendix D.5 fordetails).5.3 Trilayer: Alternating twistWe first consider the case of trilayers in the alternating twist configuration, where the topand bottom layers are twisted in the same direction, as shown in Fig. 5.1(a). This stacking395.3. Trilayer: Alternating twistFigure 5.1: Spontaneous T breaking in twisted trilayer nodal superconductors. The toppanels (a-d) depict data for the alternating twist (AT) stacking and the bottom panels(e-h) for the chiral twist (CT) stacking. (a,e): Lattice geometries. (a): The AT stackingis obtained by twisting successive layers by ±θ, and is invariant under an out-of-planemirror reflection Mz about the middle layer. (e): The CT stacking is obtained by twistingsuccessive layers by the same angle θ with respect to a common origin in the plane and isinvariant under the C2 symmetry that implements a π rotation about the diagonals of themiddle layers. (b,f): Contour plot of the BdG free energy (in arbitrary units) for θ = 45◦in the space of the phases φ1, φ3 of the top and bottom layer order parameters (settingφ2 = 0). (c,g): Low-energy spectrum in one quadrant of the BZ for large twist angles. TheCT stacking is fully gapped, whereas the AT stacking is gapless due to the Dirac conedispersion of the decoupled monolayer. (d,h): Phase differences φ1 and φ3 that minimizethe BdG free energy at zero temperature, and the spectral gap (grey curve) as a functionof the twist angle. In (d), the physics is similar to that of a standalone bilayer, except thatthe BdG spectrum is gapless owing to the decoupled Dirac cone. In (h), the non-trivialphases break T but preserve the product C2T . The corresponding dispersion is gappedexcept at 45o. The parameters are chosen as ϵc = 60meV, |∆l| ≈ 40meV and g = 20meV.405.3. Trilayer: Alternating twistis invariant under a mirror reflection symmetry Mz that interchanges the top and bottomlayers.5.3.1 Spontaneous T -breakingIn the AT geometry, the self-consistency conditions Eq. 5.5 yield a solution with |∆1| = |∆3|as expected from the mirror symmetryMz, while |∆2| is reduced by a factor |∆1| − |∆2| ∼g2/|∆1|. Setting the order parameter of the middle layer to be real (φ2 = 0), we find thatthe preferred solution is always of the form φ1 = φ3. Much like the bilayer, non-trivialphase differences φ1 ̸= 0, π that spontaneously break T emerge for twist angles aroundθ = π/4. Interestingly, however, the BdG spectrum always remains gapless.These observations, summarized in Fig. 5.1(a-d), can be understood by adapting resultsfrom alternating-twist multilayer graphene [38]. Through a unitary transformation detailedin Appendix D.1, the Hamiltonian describing the AT trilayer can be reduced to twodecoupled blocks describing (i) a twisted bilayer with a renormalized interlayer couplingg˜ =√2g and (ii) a monolayer d-wave SC with the usual four Dirac cones in its quasiparticlespectrum. The monolayer and bilayer blocks are respectively odd and even under themirror reflection Mz [38], and thus cannot hybridize as long as Mz is preserved. It thenfollows that the bilayer block admits a T -broken phase near θ = 45◦ via the phenomenologyoutlined in the previous chapters. Since the monolayer contributes spectator Dirac cones,the superconducting state however remains gapless throughout the phase diagram.In the T -broken phase, the bilayer sector is characterized by a non-zero Chern number|C| = 2 or 4, depending on the value of parameters [47]. Evaluating C using the latticemodel for twisted trilayers confirms this expectation, while an infinite strip geometry showschiral edge modes as shown in Fig. 5.2(a). One subtlety in interpreting the edge modespectra is that the monolayer sector contributes zero-energy edge modes that connectprojections of the bulk Dirac cones, in analogy to zigzag graphene nanoribbons [141].Chiral propagating states and zero-energy modes carry opposite eigenvalues under Mzand thus cannot scatter onto each other by symmetry-preserving impurities.5.3.2 Quadratic band touchingThe decoupled twisted bilayer sector also implies the presence of magic angles θM whereDirac cones merge pairwise to form quadratic band touchings (QBTs), as described inRef. [48]. To see this we linearize the Hamiltonian (5.3) around the nodes by assumingsmall angles and real order parameters. In a monolayer d-wave superconductor, the Diracnodes occur at the Fermi momenta kF where the order parameter changes sign. Expandingin small momenta q around the node in the first quadrant of the BZ, the low-energy theoryreadshq = (v∆q⊥)τ1 + (vF q∥)τ3, (5.6)where τ are Pauli matrices acting in Nambu space and q∥ (q⊥) denotes momentumcomponent parallel (perpendicular) to kF at the node. The Fermi velocity is defined asvF = kF /me and the velocity in the perpendicular direction v∆ = ∂αk∆k/kF is evaluated at415.3. Trilayer: Alternating twistFigure 5.2: Spectrum of the AT trilayer (left) and the CT trilayer (right) on an infinitecylinder geometry for a commensurate twist angle θ1,2 = 36.9◦ (see Appendix D.5 formodelling details). The color scale shows the expectation value ⟨yˆ⟩ of the eigenstates, withthe length y along the open direction normalized to unity. The AT trilayer shows 4 chiraledge modes (contributed by the bilayer sector) that coexist with zero-energy edge states(contributed by the d-wave monolayer sector) connecting the projection of the bulk Diraccones. In contrast, the CT trilayer shows a fully gapped bulk with 6 chiral edge modes.The unit cells in the two configurations have 15 and 75 sites respectively. Both calculationsuse 90 unit cells along the open direction with parameters t = 1, µ = −1.3, ∆ = 0.6 andg = 0.5.the node. The low-energy theory of a twisted multilayer is obtained by similarly linearizingall the nodes and coupling them via an interlayer tunneling matrix.As the twist angle is tuned, two of the Dirac points in each quadrant merge into aquadratic band touching point (QBT), illustrated in Fig. 5.3(a). The magic angle is givenby θM = 2g˜/v∆kF , where g˜ =√2g is the effective tunneling parameter [48]. As outlined inAppendix D.2, the low-energy theory describing this situation can be obtained throughdegenerate perturbation theory and the effective Hamiltonian is given byheffq = d1(q)η1 + d3(q)η3, (5.7)where η are Pauli matrices acting in the degenerate subspace of zero-energy solutions atthe QBT, andd1(q) =q2∥v2F − q2⊥v2∆2g˜, d3(q) =q∥q⊥vF v∆g˜. (5.8)The low-energy spectrum is quadratic and anisotropic,Eq =12g˜(q2∥v2F + q2⊥v2∆). (5.9)Such two-dimensional QBTs have a non-vanishing density of states that renders themsusceptible to gap opening through spontaneous generation of a mass term [142, 143]. From425.3. Trilayer: Alternating twistFigure 5.3: Dispersion of BdG quasiparticles in one quadrant of the Brillouin zone, absentany gap-opening instabilities, as the system is tuned across the magic angle θM for trilayersin the AT (top row) and CT (bottom row) configurations. (a): The decoupling of ATtrilayers in two sectors is manifest by the merging of two Dirac cones in a quadraticband touching (QBT) point, contributing a finite density of states at E = 0 as shown in(b), along with a spectator Dirac cone. (c): The merging of three Dirac cones along thediagonals at q = 0 is protected by the C2 symmetry of the CT stacking and leads to a cubicband touching (CBT) with a divergent density of states D(E) ∼ E−1/3 at low energies, asshown in (d).the perspective of the effective Hamiltonian (5.7) such a mass Mη2 breaks time-reversalsymmetry and can be generated in two ways.First is via residual interactions – that is, interactions not taken into account by the BCSmean-field treatment of nodal superconductivity. While such sub-dominant interactions areirrelevant in the renormalization group sense for a Dirac dispersion characteristic of d-wavesuperconductors, they become marginal at a QBT. In the simplest case of an attractives-wave channel with strength U , a mean-field analysis in Appendix D.3 shows that thefree energy is minimized by a d+ is phase, which is topologically trivial with zero Chernnumber although it clearly breaks T . The induced gap scales as e−1/UD(0) for small U ,where D(0) = 2g˜/πvF v∆ ∼ g˜/∆ is the density of states at zero energy at the QBT (seeAppendix D.4). The gap is thus exponentially small in both the interaction U and theeffective tunneling strength g˜/∆. It underlies the numerical difficulty in stabilizing thisphase, unless using an unphysically large U that is on verge of destabilizing the primaryd-wave order [144].A T -breaking gap can also open without assuming additional pairing channels beyondthe leading dx2−y2 channel. Indeed, Cooper pair tunneling between twisted layers generatesan effective dxy component that can combine with the native dx2−y2 order with a relativephase, via the phenomenology that underlies T -breaking at large twist angles [47]. For435.3. Trilayer: Alternating twistadb cfeFigure 5.4: Self-consistent mean-field results for alternating-twist (a-c) and chiral-twist(d-f) trilayers. (a,d): Zero-temperature phase diagrams for g = 12meV, showing thelayer-resolved order parameter amplitudes |∆l|, phase factor φ1 and the spectral gap Egapas a function of θ. The T -breaking phase is nucleated near the magic angles θM (verticaldotted gray lines), where the density of states is maximal, and persists to large twist anglesup to π/2 − θM. For AT trilayers, the spectrum remains gapless throughout due to thedecoupled monolayer sector. (b,e): Phase diagrams showing φ1/π in the g − θ plane. Theanalytical magic angle conditions are indicated by dashed gray lines and track the onsetof the topological regime for small twist angles. The orange and purple lines show theg = 12meV cut used in (c,f) and the g = 20meV cut used in Fig. 5.1, respectively. (c,f):Phase diagrams showing φ1/π in the T − θ plane, with the temperature normalized by thecritical temperature Tc for zero twist angle. The T -broken domes where topological SC isstabilized around the magic angles are delineated in (f).generic small angles, this mechanism is not operative because the dxy perturbation isirrelevant at the Dirac points. However, near the magic angles the finite density of statesat the QBTs leads to the nucleation of a d+ id′ phase. For large enough g ≳ 10meV thisdome merges with that around 45◦ to give a topological region that extends for all twistangles θM < θ < π/2− θM, as seen in Fig. 5.4(b). This effect is enhanced in AT trilayers ascompared to twisted bilayers because of the√2 increase in the effective tunneling strength,which leads to a larger magic angle and an enhanced density of states at the QBT.Ultimately, the nature of the superconducting order parameter near the QBT willdepend on the competition between the twist-angle-induced dxy channel and other possiblesub-dominant pairing channels specific to the material under consideration.445.4. Trilayer: Chiral twist5.4 Trilayer: Chiral twist5.4.1 Spontaneous T -breakingWe now consider trilayer systems with a constant twist θ between successive layers. Thechiral stacking has a C2 symmetry corresponding to a π rotation about diagonals of themiddle layer; see Fig. 5.1(e). In such a geometry, the problem can no longer be simplifiedthrough a unitary transformation, as in the alternating twist case, and is thus qualitativelydifferent from twisted bilayers.Once again choosing the order parameter of the middle layer to be real (φ2 = 0), wefind that as θ is increased the free energy develops two equivalent, time-reversed minimaat φ1 = −φ3 ̸= 0, as shown in Fig. 5.1(f). This non-trivial phase structure describes aT -broken state which nevertheless preserves the product C2T . The system is gapped exceptat the special point θ = 45◦, which is gapless because the CT and AT stackings are thenidentical up to a rotation of the top layer by π/2, which simply contributes an additionalphase difference of π given the C4 symmetry of the d-wave order parameter. We note thatfor identical parameters g, V and ϵc, the CT stacking exhibits T -breaking in a larger rangeof twist angles than the AT stacking, a trend which also extends to multilayers as discussedin Sec. 5.5.In the T broken phase, the lattice model shows a non-zero Chern number that takesvalues as high as |C| = 6, depending on the choice of parameters. The correspondingchiral edge modes that traverse the bulk gap are seen in Fig. 5.2. Similar to the physics oftwisted bilayers (whose superconducting state also respects C2T symmetry), the Chernnumber assignment can be understood by noting that each layer contributes C = ±2when in a d± id′ phase. Contrast this to AT stacking where, due to the presence of Mzmirror symmetry and aforementioned unitary decoupling, only two out of the three layerscontribute to topological superconductivity and give rise to a maximal Chern number|C| = 4.5.4.2 Cubic band touchingWe now turn to CT trilayers with small twist angles. The C2 symmetry dictates that thethree Dirac cones in a BZ quadrant must be arranged symmetrically about the diagonalskx = ±ky. As illustrated in Fig. 5.3(c), at small twist angles all three Dirac cones are onthe diagonals, while at large twist angles one Dirac cone remains on the diagonal and theother two are located at an equal distance on either side. The transition between thesetwo cases occurs at the magic angle θM =√2g/v∆kF where the three Dirac cones mergeinto a cubic band touching (CBT).Starting from Eq. (5.6) we can project down to the zero-energy subspace at the CBTusing degenerate perturbation theory (outlined in Appendix D.2). The resulting low-energyeffective Hamiltonian reads heffq = d1(q)η1 + d3(q)η3 withd1(q) =q⊥v∆(q2⊥v2∆ − 3q2∥v2F )4g2, d3(q) =q∥vF (q2∥v2F − 3q2⊥v2∆)4g2. (5.10)455.4. Trilayer: Chiral twistFigure 5.5: Self-consistent mean-field results for quadrilayers in the AT (left) and CT (right)configuration, using g = 12meV, |∆l| ∼ 40meV and ϵc = 60meV and zero temperature. Fortwist angles around θ = π/4 the SC phase breaks T for both stackings. (Left:) The magicangles are determined by the effective tunneling in the bilayer blocks g˜1,2 = (√5∓ 1)g/2and are denoted by dashed gray lines. For these parameters T -breaking sets in only atthe second magic angle. (Right:) For CT quadrilayers, chiral SC which preserves the C2Tsymmetry is stabilized for a large range of twist angles.Note that unlike in the AT twist, here the bare interlayer tunneling g appears and thespectrum takes the anisotropic cubic formEq =14g2(q2∥v2F + q2⊥v2∆)3/2. (5.11)Such a CBT in two dimensions admits a divergent density of states at low energyD(E) ∼ νE−1/3 with ν = g4/3/vF v∆, as illustrated in Fig. 5.3(d), which makes the systemsensitive to a gap opening. Indeed, as derived in Appendix D.4, in the presence of asub-dominant s-wave pairing channel of strength U a T -breaking gap scaling as a power law(Uν)3 ∼ U3g4/(v∆vF )3 is induced. This is a stronger dependence than the exponentiallysmall gap expected for the QBTs. However, the gap is small in practice because of thesuppression by a large power of both U and the interlayer tunneling g.In the absence of additional pairing channels, the self-consistent numerical solutionshows that the system develops a d+ id′ phase induced by interlayer tunneling for mosttwist angles, as shown in Fig. 5.4(e). For sufficiently large g the magic angle θM marksthe boundary of the d+ id′ phase, which survives all the way to π/2 − θM. For smallerg there are two separate pockets of d + id′ superconductivity (one centered around θMand the other centered around θ = 45◦) that are separated by a topologically trivial(T -preserving) region. The topological superconductivity near θ = 45◦ is robust, with acritical temperature close to the native Tc of a monolayer, while the domes surroundingthe magic angles in Fig. 5.4(f) are more fragile and persist only up to ∼ 0.05Tc for optimalg and θ.Finally, we stress that because interactions are expected to be relevant at a CBT, atreatment of interactions beyond mean-field might be necessary to determine the fate ofthe system. In particular, the scale of the interaction-induced gap could be larger thanthat anticipated through our simple BCS treatment.465.5. Generalization to multilayersa b c dFigure 5.6: Phase diagram comparison of chiral twist multilayers. The color scale indicatesthe T -breaking gap in (a) bilayers, (b) trilayers, (c) quadrilayers and (d) pentalayers. Theextent of the topological region of the phase diagram increases with L, and the maximalgaps are obtained away from θ = π/4. We choose parameters ϵc = 60meV, |∆l| ≈ 40meVand T = 0. The grey areas in c and d denote regions where numerical convergence isdifficult due to the small energy scales involved.5.5 Generalization to multilayersWith an understanding of spontaneous T breaking and higher-order band touchings intrilayers, we generalize our analysis to multilayer stacks.5.5.1 Alternating twistIn an AT stacking, layer l is rotated by θl = (−1)lθ/2 such that neighboring layers haverelative twists ±θ. For a trilayer, the problem boiled down to a twisted bilayer and aspectator monolayer. More generally, as outlined in Appendix D.1, the decoupling unitarytransformation holds for any number of layers and the system reduces to ⌊L/2⌋ bilayerscharacterized by the renormalized interlayer tunneling amplitudesg˜m = 2g cos(πmL+ 1). (5.12)Here m = 1 . . . ⌊L/2⌋ and ⌊.⌋ is the floor function. Additionally, in the case of odd L onehas a decoupled monolayer sector. When the twist angle is close to 45◦, each of the bilayerblocks is expected to spontaneously break T , leaving any decoupled layer as is – the layersalways gap out in pairs. In the same vein, the physics at small twist angles carries overand each bilayer is characterized by a QBT at a different magic angle θmM = 2g˜m/v∆kF .The decoupling unitary transformation relies on the assumption that the order pa-rameters in all layers of a given parity (even/odd) carry the same phase and amplitude.While such an assumption is not enforced by symmetry, it is approximately respected inthe self-consistent solutions of AT multilayers, with small differences between |∆l| of orderg2/∆. The lack of a symmetry constraint also means that the bilayer blocks can weaklyhybridize with each other, which introduces quantitative differences in the numerical valuesof magic angles, but does not qualitatively alter the physics.An analysis of quadrilayers, summarized in Fig. 5.5(a), bears out these expectations.We find that the condition for unitary decoupling is approximately respected: ∆1 ≈ ∆3475.5. Generalization to multilayersL NM Low-energy spectrum at θM1 0 n/a2 1 QBT3 1 CBT4 2 DQBT, DQBT5 2 DQBT + Dirac cone, DQBT + Dirac cone6 3 DQBT, QBT, DQBT7 3 DQBT, CBT, DQBT8 4 DQBT, DQBT, DQBT, DQBTTable 5.1: In chirally twisted multilayers with weak interlayer tunneling, the low energyspectra at the sequence of magic angles follow a mod 4 pattern as a function of the numberof layers L. Here NM denotes the number of magic angles, while QBT and CBT standfor quadratic and cubic band touchings at q = 0, respectively. DQBT stands for doubleQBT, and denotes two QBTs that occur at the same twist angle, i.e., four Dirac conesmerge pairwise as shown in Fig. 5.7. When more than one magic angle occurs we list theircorresponding low-energy spectra in order of increasing θM.and ∆2 ≈ ∆4. Close to θ = 45◦ a non-trivial phase difference develops between even andodd numbered layers, accompanied by a T -breaking gap. With our choice of parameters,T -breaking appears only at the second magic angle, presumably because the density ofstates at the first magic angle, which scales as D(0) ∼ g˜m, is too small to nucleate a d+ id′phase. In a twisted bilayer the phase difference determines the sign of the Chern number.Given the small phase difference within the odd (even) layer subspace, we deduce that bothbilayer blocks break T in the same way such that the total Chern number C = ±4⌊L/2⌋.This expectation is indeed borne out by our lattice calculations, which find a Chern number|C| = 8 in the T -broken phase of AT quadrilayers.5.5.2 Chiral twistDue to the lack of a simplifying unitary transformation, the generalization to multilayersis less straightforward in the chiral twist case, which corresponds toθl =(l − L+ 12)θ. (5.13)We thus explore this configuration numerically, and find that for a wide range of twistangles the system spontaneously breaks T while respecting C2T , as depicted in Fig. 5.5(b)for quadrilayers. Interestingly, the extent of the T -breaking region of the phase diagramincreases with L (see Fig. 5.6). In contrast to a bilayer, the optimal twist angle fortopological superconductivity occurs away from θ = 45◦. The topological gap at optimaltwist angle decreases with the number of layers L but the dependence is weak. Guidedby the intuition that in a chirally-twisted stack all layers (except those at the very topand bottom) experience effectively the same local environment, we expect each layer tocontribute the same value Cl = ±2 to the total Chern number. Supported by lattice485.6. OutlookFigure 5.7: Evolution of the BdG quasiparticle dispersion of chirally-twisted quadrilayers atsmall angles. There are two magic angles θ(1)M and θ(2)M where two quadratic band touchingsoccur simultaneously at momenta ±q(1)M and ±q(2)M respectively.calculations on trilayers, we thus conjecture that generically C = ±2L in the T -brokenphase.At small twist angles, the continuum model exhibits an intriguing pattern of bandtouchings linked to L mod 4, summarized in Table 5.1. We find that CBTs occur onlywhen L mod 4 = 3, and also note the occurrence of double QBTs: two quadratic bandtouchings in a BZ quadrant occurring at the same twist angle. In a quadrilayer, for instance,there are two magic angles θ(1)M and θ(2)M where double QBTs arise (see Fig. 5.7).5.6 OutlookTwisted nodal superconductors are a promising platform to realize topological supercon-ducting phases. Here we extended these ideas to multilayers in two different geometriesand mapped their phase diagrams.At small twist angles, the quasiparticle spectra of multilayers exhibit a sequence ofmagic angles with either quadratic, cubic or double quadratic dispersions, depending on thestacking and the number of layers L. The constant (power-law divergent) density of statesof the quadratic (cubic) band touchings renders these systems susceptible to secondarysymmetry-breaking transitions in the superconducting state, which can be triggered eitherby interlayer tunneling or residual interactions. While interlayer tunneling genericallyfavors T breaking chiral superconductivity with an order parameter of d+ id′ symmetry,the ordering tendencies due to the residual interactions will depend on microscopics. Afull accounting of instabilities in pairing and various particle-hole channels would require adetailed understanding of interactions in a given material and might necessitate a treatmentbeyond mean-field.The cubic and double quadratic band touchings constitute new features that go beyondtwisted bilayers. Cubic band touchings have a divergent density of states at zero energy,D(E) ∼ E−1/3, and residual interactions are expected to be relevant from a renormalizationgroup perspective. This is reminiscent of higher-order van Hove singularities occurring, forinstance, in the normal state dispersion of cuprates [145], Moire´ surface states of topologicalinsulators [146] and transition metal dichalcogenides [147]. At a double quadratic bandtouching, four Dirac cones in the same BZ quadrant merge pairwise to give rise to two495.6. Outlookquadratic crossings. The proximity of these nodes in momentum space could potentiallylead to finite-momentum pairing instabilities.Closer to θ = 45◦, both configurations spontaneously break T but in qualitativelydifferent ways. AT multilayers can be reduced to a set of bilayers with renormalizedinterlayer tunneling strengths (alongside a decoupled monolayer for odd L). Their physicsis therefore similar to that of a twisted bilayer [47, 140, 148]. Nevertheless, in view of theestimates for the tunneling strength g ∼ 5−10meV inferred from recent transport studies [50,144, 149], such a renormalization might help stabilize T -broken superconductivity. Incontrast, the problem posed by CT multilayers cannot be simplified in such a manner.While T and a π rotation C2 along the diagonals are individually broken, their productC2T is preserved. The system exhibits chiral topological superconductivity in an extendedregion of the twist-angle/interlayer tunneling phase diagram with a maximal Chern numberC = ±2L.The connection of our results to cuprate heterostructures, explored in recent transportexperiments [49–51], relies on two key simplifications: the mean-field treatment of super-conductivity and use of a momentum-independent interlayer tunneling form factor. Thelatter point was the subject of recent discussion in the context of twisted Bi2212 [140].Close to 45◦ twist, the fate of the SC state will depend sensitively on the symmetries ofthe orbitals involved in interlayer tunneling processes. In clean and aligned bilayers thedirect (momentum-independent) tunneling between d orbitals vanishes by symmetry at 45◦twist angle. There the dominant contribution is expected to come from s-orbital mediatedtunneling, which is described by a form factor [94] glmk ∼ cos(2αk − θl) cos(2αk − θm) thatunfortunately vanishes along the all-important nodal directions where the T -breaking gapopens [47]. In the context of chirally-twisted multilayers, the vanishing of direct tunnelingat 45◦ might be less of a problem because, as illustrated in Fig. 5.6, the region of the phasediagram with robust topological superconductivity is not centered around 45◦.In bulk cuprates, incoherent tunneling mediated by disorder in layers interspacingCuO2 planes is expected to play a crucial role in interlayer transport [129–132]. Therefore,incoherent tunneling in twisted layers is of potential interest. Because disorder breakstranslation invariance, tunneling processes are not subject to symmetry constraints invokedin Ref. [140] and may allow for substantial twist-induced gaps even close to 45◦. A simpleestimate from an impurity-assisted tunneling model shows that a gap of order 1− 10meVcould be generated in twisted bilayers [133].Also, note that the analysis presented here assumes one CuO plane per layer. Currentexperiments, however, involve Bi2212 which has two CuO planes per monolayer. Based onthe previous analysis, we expect the phenomenology of T -breaking to hold in this situation,albeit with two quantitative modifications: (i) the scale of the gap, which is determinedby a competition between inter-monolayer tunneling and intra-monolayer physics, will godown, and (ii) Chern numbers in the T -broken phases will be doubled.In terms of experimental signatures of the topological phases, perhaps the most promis-ing for cuprate heterostructures is optical Hall conductivity probed by the polar Kerr effect,which was predicted to occur with a large magnitude owing to the gap scales involved [150].Importantly, the Kerr effect is sensitive not only to T -breaking but also to the breaking ofmirror symmetries, which allows it to distinguish between a chiral state such as d+ id′505.6. Outlookand a trivial state such as d + is. The predicted topological gap could additionally bedetected through spectroscopies such as angle-resolved photoemission, electron tunneling,or optically using Raman response. Transport experiments such as the superconductingdiode effect, whereby the critical supercurrent becomes direction-dependent, can also inprinciple discriminate between chiral and non-chiral T -broken states [151].Taking a broad view, the phenomenology discussed here could be relevant to other2D nodal superconductors that can be exfoliated and stacked. The material-specificmodelling of such platforms is left for future work. An interesting question concerns thecreation of new heterostructures that realize an odd number of Majorana modes, eitherpropagating along boundaries or localized in vortex cores (see for instance Ref. [139]).The limit of a large number of layers L is also of potential interest. On the one hand,similar to alternating-twist graphene [38], the ‘cascade’ of magic angles in alternating-twistmultilayers as L is increased flattens the low-energy BdG bands, which could make it easierto stabilize symmetry-broken phases induced by interactions. On the other hand, theinfinite-layer limit of the chiral twist stacking exhibits a non-symmorphic symmetry givenby a combination of twist and translation in the stacking direction. A generalization of ‘3Dtwistronics,’ developed in the context of chirally-twisted graphene multilayers [152–154],to nodal superconductors thus represents another direction for future studies.51Chapter 6ConclusionsA superconductor with px + ipy order has long fascinated the physics community becausevortex defects in such a system host Majorana zero modes. In Chapter 2 we propose asimple construction of a chiral superconductor using proximitized quantum wires and twistangle engineering as basic ingredients. A weakly coupled parallel array of these wiresforms a gapless p-wave superconductor. Two such arrays, stacked on top of one anotherwith a twist angle close to 90◦, spontaneously break time-reversal symmetry and producea px + ipy superconductor. We have mapped the topological phases of this setup anddemonstrated the existence of Majorana zero modes in vortices.As another application of twist angle engineering, in Chapter 3 we consider structurescomposed of two monolayer-thin d-wave superconductors with a twist angle θ that canbe realized by mechanically exfoliating van der Waals-bonded high-Tc copper oxide mate-rials, such as Bi2Sr2CaCu2O8+δ (BSCCO). Based on symmetry arguments and detailedmicroscopic modelling, we show that for a range of twist angles in the vicinity of 45o,such bilayers form a robust, fully gapped topological phase with spontaneously brokentime-reversal symmetry and protected chiral Majorana edge modes. Importantly, at 45othe topological phase can, in principle, set in close to the bulk Tc that could be as high as90K in BSCCO.Given these predictions and ongoing experiments, a key challenge is to identify un-ambiguous signatures of these topological phases. With this goal in mind, Chapter 4explores a suite of Josephson phenomena between twisted cuprate layers. At intermediatetwist angles, we find an unusual non-monotonic temperature dependence of the criticalcurrent which can be attributed to the unconventional sign structure of the d-wave orderparameter. The onset of the T -broken phase near the 45◦ twist is marked by a crossoverfrom the conventional 2π-periodic Josephson relation to a π-periodic function because thesingle-pair tunneling becomes dominated by a second order co-tunneling process. Despitethis fundamental change, the critical current remains a smooth function of the twist angleand temperature implying that a measurement of Jc alone will not be a litmus test forthe T -broken phase. We show that clear signatures of the T -broken phase appear whenJc is measured in the presence of an applied magnetic field or radio-frequency drive: theresulting Fraunhofer oscillations acquire additional nodes and fractional Shapiro stepsbecome visible. These results are broadly consistent with recent experiments on twistedbilayers of BSCCO.In Chapter 5, we extend the twist analysis to multilayers of nodal superconductors,focusing on two high-symmetry stackings with alternating (±θ) and constant (θ) twistangles. In analogy to alternating-twist multilayer graphene, the former can be mappedto twisted bilayers with renormalized interlayer couplings, along with a remnant gaplessmonolayer when the number of layers L is odd. In contrast, the latter exhibits physics52Chapter 6. Conclusionsbeyond twisted bilayers, including the occurrence of magic angles characterized by cubicband crossings when L mod 4 = 3. Owing to their power-law divergent density of states,such multilayers are highly susceptible to secondary instabilities. Within a BCS mean-fieldtheory, defined in the continuum and on a lattice, we find that both stackings host chiraltopological superconductivity in extended regions of their phase diagrams.Before closing, we comment on a few open problems and possible avenues for futureresearch: Theoretically, it is a well-known fact that cuprates are not quantitatively describedby BCS theory. A more sophisticated treatment of correlations in twisted bilayers using at-J model, and the inclusion of an interlayer tunneling form factor appropriate to Bi2212,suggests a topologically trivial phase near 45◦ that is either gapless or has a very smallT -breaking gap [140]. On the other hand, numerical treatment of the Hubbard modelon a twisted square lattice at the commensurate angle θ = 53.1◦ via the variationalcluster approximation finds a T -broken phase that is gapped yet topologically trivial [148].Discrepancies between these theoretical models highlight the importance of complex orbitalstructures and strong electronic correlations inherent to cuprates.Moreover, interlayer transport in cuprates indicates the importance of disorder mediatedincoherent tunneling which evades symmetry constraints present in an idealized crystal[120, 128]. This arises even in the cleanest single-crystal samples through oxygen vacanciesin layers separating the CuO2 planes. These are introduced to achieve the hole dopingnecessary for superconductivity. In a separate work [133], we assessed the influence ofincoherent tunneling on the phase diagram of a twisted bilayer and show that the modelcontinues to support a fully gapped topological phase with broken T . Compared to themodel with a momentum conserving interlayer coupling, the extent of the topological phasearound the 45◦ twist decreases with increasing incoherence but remains robustly presentfor parameters likely relevant to BSCCO.So far, we have focused on 2D superconductors whose 1D boundaries host gapless statesdue to nontrivial topology of the bulk. As per the nomenclature in Ref. [155], these aretermed first order topological superconductors. Second order topological superconductorsin 2D would have gapless modes at corners. It would be worthwhile to consider whethertwisted structures could exhibit higher order topology. Higher-order topological protectionrelies on spatial symmetries and constitutes an extension of the notion of topologicalcrystalline phases of matter [156]. One potential candidate is the cuprate bilayer exactlyat 45◦ twist, where the lattice forms an eight-fold rotation symmetric quasicrystal.While Majorana edge states in the chiral d-wave superconductor are not directlyapplicable in quantum computation, it has been shown that when a twisted bilayer in thed+ id′ phase is proximitized with a strong topological insulator, vortex defects again hostMajorana zero modes [139]. This motivates the consideration of other hetero-structureswhich take advantage of the topological gap that can persist to high temperatures. Forinstance, Majoranas were predicted in proximitized topological insulator nanowires [157].Proximity coupling to the chiral state could similarly engender interesting physics.On the experimental front, transport measurements in few-layer thick samples ofBSCCO find that close to 45◦ twist angle the interlayer Josephson current is dominatedby its second harmonic, which signals co-tunneling of Cooper pairs [50]. While indicativeof a non-trivial phase difference between the two superconducting layers, and hence T53Chapter 6. Conclusionsbreaking, this does not directly reveal the topological nature of the state [144, 149]. Alongsimilar lines, Ref. [51] notes a vanishing critical current at 45◦. Another experimentalwork, however, interprets transport data as evidence for a subdominant s-wave pairingchannel [49]. Clearly, the findings and interpretations are at odds with each other. Onemay attribute this to the differences in sample thickness, fabrication techniques andinterface quality. Spontaneous time-reversal symmetry breaking is essentially an interfacephenomenon, meaning that it is most visible at the layers closest to the twist interface.Probing the interface is harder in thick samples. Further, unlike graphene, cupratesare sensitive to the atmosphere, so special techniques, such as low temperature samplefabrication, are necessary to preserve the properties of crystals [50]. While the observationof second harmonic in Josephson phenomena is a promising start, additional probes such aspolar Kerr angle measurements [150, 158, 159] and, possibly, edge mode detection throughthermal Hall measurements will be necessary to ascertain the topological character.54Bibliography[1] K. v. Klitzing, G. Dorda, and M. Pepper, Phys. Rev. Lett. 45, 494 (1980).[2] Q. Niu, D. J. Thouless, and Y.-S. Wu, Phys. Rev. B 31, 3372 (1985).[3] F. D. M. Haldane, Phys. Rev. Lett. 61, 2015 (1988).[4] C. L. Kane and E. J. Mele, Phys. Rev. Lett. 95 (2005), 10.1103/physrevlett.95.226801.[5] M. Konig, S. Wiedmann, C. Brune, A. Roth, H. Buhmann, L. W. Molenkamp, X.-L.Qi, and S.-C. Zhang, Science 318, 766 (2007).[6] I. Knez, R. R. Du, and G. Sullivan, Phys. Rev. B 81 (2010), 10.1103/phys-revb.81.201301.[7] C.-Z. Chang, J. Zhang, X. Feng, J. Shen, Z. Zhang, M. Guo, K. Li, Y. Ou, P. Wei,L.-L. Wang, Z.-Q. Ji, Y. Feng, S. Ji, X. Chen, J. Jia, X. Dai, Z. Fang, S.-C. Zhang,K. He, Y. Wang, L. Lu, X.-C. Ma, and Q.-K. Xue, Science 340, 167 (2013).[8] M. Serlin, C. L. Tschirhart, H. Polshyn, Y. Zhang, J. Zhu, K. Watanabe, T. Taniguchi,L. Balents, and A. F. Young, Science 367, 900 (2020).[9] G. Chen, A. L. Sharpe, E. J. Fox, Y.-H. Zhang, S. Wang, L. Jiang, B. Lyu, H. Li,K. Watanabe, T. Taniguchi, Z. Shi, T. Senthil, D. Goldhaber-Gordon, Y. Zhang,and F. Wang, Nature 579, 56 (2020).[10] A. Altland and M. R. Zirnbauer, Phys. Rev. B 55, 1142 (1997).[11] S. Ryu, A. P. Schnyder, A. Furusaki, and A. W. W. Ludwig, New Journal of Physics12, 065010 (2010).[12] D. A. Ivanov, Phys. Rev. Phys. Rev. 86, 268 (2001).[13] M. Sato and Y. Ando, Reports on Progress in Physics 80, 076501 (2017).[14] A. Y. Kitaev, Physics-Uspekhi 44, 131 (2001).[15] C. Nayak, S. H. Simon, A. Stern, M. Freedman, and S. D. Sarma, Reviews of ModernPhysics 80, 1083 (2008).[16] H. Tou, Y. Kitaoka, K. Asayama, N. Kimura, Y. Onuki, E. Yamamoto, andK. Maezawa, Phys. Rev. Phys. Rev. 77, 1374 (1996).[17] L. Fu and E. Berg, Phys. Rev. Lett. 105, 097001 (2010).55Bibliography[18] Y. S. Hor, A. J. Williams, J. G. Checkelsky, P. Roushan, J. Seo, Q. Xu, H. W.Zandbergen, A. Yazdani, N. P. Ong, and R. J. Cava, Phys. Rev. Phys. Rev. 104(2010).[19] P. Zhang, K. Yaji, T. Hashimoto, Y. Ota, T. Kondo, K. Okazaki, Z. Wang, J. Wen,G. D. Gu, H. Ding, and S. Shin, Science 360, 182 (2018).[20] Z. Wang, J. O. Rodriguez, L. Jiao, S. Howard, M. Graham, G. D. Gu, T. L. Hughes,D. K. Morr, and V. Madhavan, Science 367, 104 (2020).[21] X.-L. Qi, T. L. Hughes, and S.-C. Zhang, Phys. Rev. B 82, 184516 (2010).[22] J. Wang, Q. Zhou, B. Lian, and S.-C. Zhang, Phys. Rev. B 92, 064520 (2015).[23] K. Kuboki and M. Sigrist, Journal of the Physical Society of Japan 65, 361 (1996).[24] R. Nandkishore, L. S. Levitov, and A. V. Chubukov, Nature Physics 8, 158 (2012).[25] Y. Oreg, G. Refael, and F. von Oppen, Phys. Rev. Phys. Rev. 105 (2010).[26] R. M. Lutchyn, J. D. Sau, and S. D. Sarma, Phys. Rev. Phys. Rev. 105 (2010).[27] V. Mourik, K. Zuo, S. M. Frolov, S. R. Plissard, E. P. A. M. Bakkers, and L. P.Kouwenhoven, Science 336, 1003 (2012).[28] L. Fu and C. L. Kane, Phys. Rev. Lett. 100, 096407 (2008).[29] P. Yu, J. Chen, M. Gomanko, G. Badawy, E. P. A. M. Bakkers, K. Zuo, V. Mourik,and S. M. Frolov, Nature Physics 17, 482 (2021).[30] Q. L. He, L. Pan, A. L. Stern, E. C. Burks, X. Che, G. Yin, J. Wang, B. Lian,Q. Zhou, E. S. Choi, K. Murata, X. Kou, Z. Chen, T. Nie, Q. Shao, Y. Fan, S.-C.Zhang, K. Liu, J. Xia, and K. L. Wang, Science 357, 294 (2017).[31] R. Bistritzer and A. H. MacDonald, Proceedings of the National Academy of Sciences108, 12233 (2011).[32] Y. Cao, V. Fatemi, A. Demir, S. Fang, S. L. Tomarken, J. Y. Luo, J. D. Sanchez-Yamagishi, K. Watanabe, T. Taniguchi, E. Kaxiras, R. C. Ashoori, and P. Jarillo-Herrero, Nature 556, 80 (2018).[33] K. S. Novoselov, A. Mishchenko, A. Carvalho, and A. H. C. Neto, Science 353(2016).[34] M. Koshino, Phys. Rev. B 99, 235406 (2019).[35] N. R. Chebrolu, B. L. Chittari, and J. Jung, Phys. Rev. B 99, 235417 (2019).[36] J. Liu, Z. Ma, J. Gao, and X. Dai, Phys. Rev. X 9, 031021 (2019).[37] G. W. Burg, J. Zhu, T. Taniguchi, K. Watanabe, A. H. MacDonald, and E. Tutuc,Phys. Rev. Lett. 123, 197702 (2019).56Bibliography[38] E. Khalaf, A. J. Kruchkov, G. Tarnopolsky, and A. Vishwanath, Phys. Rev. B 100,085109 (2019).[39] J. M. Park, Y. Cao, K. Watanabe, T. Taniguchi, and P. Jarillo-Herrero, Nature 590,249 (2021).[40] Z. Hao, A. M. Zimmerman, P. Ledwith, E. Khalaf, D. H. Najafabadi, K. Watanabe,T. Taniguchi, A. Vishwanath, and P. Kim, Science 371, 1133 (2021).[41] F. Wu, T. Lovorn, E. Tutuc, and A. H. MacDonald, Phys. Rev. Lett. 121, 026402(2018).[42] E. C. Regan, D. Wang, C. Jin, M. I. B. Utama, B. Gao, X. Wei, S. Zhao, W. Zhao,Z. Zhang, K. Yumigeta, M. Blei, J. D. Carlstro¨m, K. Watanabe, T. Taniguchi,S. Tongay, M. Crommie, A. Zettl, and F. Wang, Nature 579, 359 (2020).[43] L. Wang, E.-M. Shih, A. Ghiotto, L. Xian, D. A. Rhodes, C. Tan, M. Claassen, D. M.Kennes, Y. Bai, B. Kim, K. Watanabe, T. Taniguchi, X. Zhu, J. Hone, A. Rubio,A. N. Pasupathy, and C. R. Dean, Nature Materials 19, 861 (2020).[44] Z. Zhang, Y. Wang, K. Watanabe, T. Taniguchi, K. Ueno, E. Tutuc, and B. J.LeRoy, Nat. Phys. 16, 1093 (2020).[45] Q. Tong, F. Liu, J. Xiao, and W. Yao, Nano Lett. 18, 7194 (2018).[46] K. Hejazi, Z.-X. Luo, and L. Balents, PNAS 117, 10721 (2020).[47] O. Can, T. Tummuru, R. P. Day, I. Elfimov, A. Damascelli, and M. Franz, Nat.Phys. 17, 519 (2021).[48] P. A. Volkov, J. H. Wilson, and J. H. Pixley, “Magic angles and current-inducedtopology in twisted nodal superconductors,” (2020), arXiv:2012.07860 .[49] Y. Zhu, M. Liao, Q. Zhang, H.-Y. Xie, F. Meng, Y. Liu, Z. Bai, S. Ji, J. Zhang,K. Jiang, R. Zhong, J. Schneeloch, G. Gu, L. Gu, X. Ma, D. Zhang, and Q.-K. Xue,Phys. Rev. X 11, 031011 (2021).[50] S. Y. F. Zhao, N. Poccia, X. Cui, P. A. Volkov, H. Yoo1, R. Engelke, Y. Ronen,R. Zhong, G. Gu, S. Plugge, T. Tummuru, M. Franz, J. H. Pixley, and P. Kim,arXiv preprint arXiv:2108.13455 (2021).[51] J. Lee, W. Lee, G.-Y. Kim, Y.-B. Choi, J. Park, S. Jang, G. Gu, S.-Y. Choi, G. Y.Cho, G.-H. Lee, and H.-J. Lee, Nano Lett. 21, 10469 (2021).[52] A. Kitaev, Annals of Physics 303, 2 (2003).[53] A. Das, Y. Ronen, Y. Most, Y. Oreg, M. Heiblum, and H. Shtrikman, Nature Physics8, 887 (2012).[54] L. P. Rokhinson, X. Liu, and J. K. Furdyna, Nature Physics 8, 795 (2012).57Bibliography[55] A. D. K. Finck, D. J. Van Harlingen, P. K. Mohseni, K. Jung, and X. Li, Phys. Rev.Lett. 110, 126406 (2013).[56] M. T. Deng, S. Vaitiekenas, E. B. Hansen, J. Danon, M. Leijnse, K. Flensberg,J. Nyg˚ard, P. Krogstrup, and C. M. Marcus, Science 354, 1557 (2016).[57] L. Wang, T. O. Rosdahl, and D. Sticlet, Phys. Rev. B 98, 205411 (2018).[58] L. Kong, S. Zhu, M. Papaj, H. Chen, L. Cao, H. Isobe, Y. Xing, W. Liu, D. Wang,P. Fan, Y. Sun, S. Du, J. Schneeloch, R. Zhong, G. Gu, L. Fu, H.-J. Gao, andH. Ding, Nature Physics 15, 1181 (2019).[59] T. Machida, Y. Sun, S. Pyon, S. Takeda, Y. Kohsaka, T. Hanaguri, T. Sasagawa,and T. Tamegai, Nature Materials 18, 811 (2019).[60] S. Zhu, L. Kong, L. Cao, H. Chen, M. Papaj, S. Du, Y. Xing, W. Liu, D. Wang,C. Shen, F. Yang, J. Schneeloch, R. Zhong, G. Gu, L. Fu, Y.-Y. Zhang, H. Ding,and H.-J. Gao, Science 367, 189 (2020).[61] J. Alicea, Reports on Progress in Physics 75, 076501 (2012).[62] C. Beenakker, Annual Review of Condensed Matter Physics 4, 113 (2013).[63] M. Leijnse and K. Flensberg, Semiconductor Science and Technology 27, 124003(2012).[64] T. D. Stanescu and S. Tewari, Journal of Physics: Condensed Matter 25, 233201(2013).[65] S. R. Elliott and M. Franz, Reviews of Modern Physics 87, 137 (2015).[66] R. M. Lutchyn, E. P. A. M. Bakkers, L. P. Kouwenhoven, P. Krogstrup, C. M.Marcus, and Y. Oreg, Nature Reviews Materials 3, 52 (2018).[67] N. Read and D. Green, Phys. Rev. B 61, 10267 (2000).[68] C. Kallin and J. Berlinsky, Reports on Progress in Physics 79, 054502 (2016).[69] A. Pustogow, Y. Luo, A. Chronister, Y.-S. Su, D. A. Sokolov, F. Jerzembeck, A. P.Mackenzie, C. W. Hicks, N. Kikugawa, S. Raghu, E. D. Bauer, and S. E. Brown,Nature 574, 72 (2019).[70] B. A. Bernevig and T. L. Hughes, Topological insulators and topological superconduc-tors (Princeton university press, 2013).[71] J.-H. Kang, A. Grivnin, E. Bor, J. Reiner, N. Avraham, Y. Ronen, Y. Cohen,P. Kacman, H. Shtrikman, and H. Beidenkopf, Nano Letters 17, 7520 (2017).[72] Y. Yu, L. Ma, P. Cai, R. Zhong, C. Ye, J. Shen, G. D. Gu, X. H. Chen, and Y. Zhang,Nature 575, 156 (2019).58Bibliography[73] Y. Cao, V. Fatemi, S. Fang, K. Watanabe, T. Taniguchi, E. Kaxiras, and P. Jarillo-Herrero, Nature 556, 43 (2018).[74] M. Yankowitz, S. Chen, H. Polshyn, Y. Zhang, K. Watanabe, T. Taniguchi, D. Graf,A. F. Young, and C. R. Dean, Science 363, 1059 (2019).[75] A. L. Sharpe, E. J. Fox, A. W. Barnard, J. Finney, K. Watanabe, T. Taniguchi,M. A. Kastner, and D. Goldhaber-Gordon, Science 365, 605 (2019).[76] X. Lu, P. Stepanov, W. Yang, M. Xie, M. A. Aamir, I. Das, C. Urgell, K. Watanabe,T. Taniguchi, G. Zhang, A. Bachtold, A. H. MacDonald, and D. K. Efetov, Nature574, 653 (2019).[77] N. N. T. Nam and M. Koshino, Phys. Rev. B 96, 075311 (2017).[78] J. Kang and O. Vafek, Phys. Rev. X 8, 031088 (2018).[79] L. Zou, H. C. Po, A. Vishwanath, and T. Senthil, Phys. Rev. B 98, 085435 (2018).[80] H. C. Po, L. Zou, A. Vishwanath, and T. Senthil, Phys. Rev. X 8, 031089 (2018).[81] M. Xie and A. H. MacDonald, Phys. Rev. Lett. 124, 097601 (2020).[82] M. Koshino, N. F. Q. Yuan, T. Koretsune, M. Ochi, K. Kuroki, and L. Fu, Phys.Rev. X 8, 031087 (2018).[83] F. Guinea and N. R. Walet, Proceedings of the National Academy of Science 115,13174 (2018).[84] J. Kang and O. Vafek, Phys. Rev. Lett. 122, 246401 (2019).[85] K. Hejazi, C. Liu, H. Shapourian, X. Chen, and L. Balents, Phys. Rev. B 99, 035111(2019).[86] E. Y. Andrei and A. H. MacDonald, Nat. Mater. 19, 1265 (2020).[87] R. B. Laughlin, Phys. Rev. Lett. 80, 5188 (1998).[88] M. Franz and Z. Tesˇanovic´, Phys. Rev. Lett. 80, 4763 (1998).[89] A. Vishwanath, Phys. Rev. Lett. 87, 217004 (2001).[90] A. Bille, R. A. Klemm, and K. Scharnberg, Phys. Rev. B 64, 174507 (2001).[91] T. Yokoyama, S. Kawabata, T. Kato, and Y. Tanaka, Phys. Rev. B 76, 134501(2007).[92] Z. Yang, S. Qin, Q. Zhang, C. Fang, and J. Hu, Phys. Rev. B 98, 104515 (2018).[93] S. Pathak, V. B. Shenoy, and G. Baskaran, Phys. Rev. B 81, 085431 (2010).59Bibliography[94] O. Andersen, A. Liechtenstein, O. Jepsen, and F. Paulsen, Journal of Physics andChemistry of Solids 56, 1573 (1995).[95] A. Damascelli, Z. Hussain, and Z.-X. Shen, Rev. Mod. Phys. 75, 473 (2003).[96] O. Fischer, M. Kugler, I. Maggio-Aprile, C. Berthod, and C. Renner, Rev. Mod.Phys. 79, 353 (2007).[97] D. Agterberg, P. Brydon, and C. Timm, Physical Review Letters 118 (2017),10.1103/physrevlett.118.127001.[98] S. C. Riggs, O. Vafek, J. B. Kemper, J. B. Betts, A. Migliori, F. F. Balakirev, W. N.Hardy, R. Liang, D. A. Bonn, and G. S. Boebinger, Nature Physics 7, 332 (2011).[99] J. L. Tallon, J. W. Loram, J. R. Cooper, C. Panagopoulos, and C. Bernhard, Phys.Rev. B 68, 180501 (2003).[100] A. Hosseini, D. M. Broun, D. E. Sheehy, T. P. Davis, M. Franz, W. N. Hardy,R. Liang, and D. A. Bonn, Phys. Rev. Lett. 93, 107003 (2004).[101] S. K. Yip and J. A. Sauls, Journal of Low Temperature Physics 86, 257 (1992).[102] B. D. Josephson, Physics Letters 1, 251 (1962).[103] C. C. Tsuei and J. R. Kirtley, Rev. Mod. Phys. 72, 969 (2000).[104] R. A. Klemm, G. Arnold, C. T. Rieck, and K. Scharnberg, Phys. Rev. B 58, 14203(1998).[105] G. B. Arnold and R. A. Klemm, Phys. Rev. B 62, 661 (2000).[106] S. Kawabata, S. Kashiwaya, Y. Asano, and Y. Tanaka, Phys. Rev. B 70 (2004).[107] E. Il'ichev, M. Grajcar, R. Hlubina, R. P. J. IJsselsteijn, H. E. Hoenig, H.-G. Meyer,A. Golubov, M. H. S. Amin, A. M. Zagoskin, A. N. Omelyanchouk, and M. Y.Kupriyanov, Phys. Rev. Lett. 86, 5369 (2001).[108] A. A. Golubov, M. Y. Kupriyanov, and E. Il’ichev, Reviews of Modern Physics 76,411 (2004).[109] Q. Li, Y. Tsay, M. Suenaga, G. Gu, and N. Koshizuka, Physica C: Superconductivity282-287, 1495 (1997).[110] Q. Li, Y. N. Tsay, M. Suenaga, R. A. Klemm, G. D. Gu, and N. Koshizuka, Phys.Rev. Lett. 83, 4160 (1999).[111] Y. Takano, T. Hatano, A. Fukuyo, A. Ishii, M. Ohmori, S. Arisawa, K. Togano, andM. Tachiki, Phys. Rev. B 65, 140513 (2002).[112] E. Goldobin, D. Koelle, R. Kleiner, and A. Buzdin, Phys. Rev. B 76 (2007).60Bibliography[113] R. S. Markiewicz, S. Sahrakorpi, M. Lindroos, H. Lin, and A. Bansil, Phys. Rev. B72, 054519 (2005).[114] V. Ambegaokar and A. Baratoff, Phys. Rev. Lett. 11, 104 (1963).[115] M. Tinkham, Introduction to superconductivity (Dover Publications, 2004).[116] A. Barone and G. Paterno`, Physics and Applications of the Josephson Effect (Wiley,1982).[117] D. J. Scalapino, S. R. White, and S. Zhang, Phys. Rev. B 47, 7995 (1993).[118] P. J. Hirschfeld and N. Goldenfeld, Phys. Rev. B 48, 4219 (1993).[119] Q.-H. Wang, J. H. Han, and D.-H. Lee, Phys. Rev. Lett. 87, 077004 (2001).[120] D. E. Sheehy, T. P. Davis, and M. Franz, Phys. Rev. B 70, 054510 (2004).[121] L. N. Bulaevskii, J. R. Clem, and L. I. Glazman, Phys. Rev. B 46, 350 (1992).[122] S. Shapiro, Phys. Rev. Lett. 11, 80 (1963).[123] P. Russer, Journal of Applied Physics 43, 2008 (1972).[124] V. Ambegaokar and B. I. Halperin, Phys. Rev. Lett. 23, 274 (1969).[125] S. P. Benz, M. S. Rzchowski, M. Tinkham, and C. J. Lobb, Phys. Rev. Lett. 64,693 (1990).[126] R. Panghotra, B. Raes, C. C. de Souza Silva, I. Cools, W. Keijers, J. E. Scheerder,V. V. Moshchalkov, and J. V. de Vondel, Communications Physics 3 (2020).[127] H. Sellier, C. Baraduc, F. Lefloch, and R. Calemczuk, Phys. Rev. Lett. 92 (2004).[128] S. L. Cooper and K. E. Gray, Physical Properties of High Temperature SuperconductorsIV , edited by D. M. Ginsberg (World Scientific, 1994).[129] M. J. Graf, D. Rainer, and J. A. Sauls, Phys. Rev. B 47, 12089 (1993).[130] R. Radtke and K. Levin, Physica C: Superconductivity 250, 282 (1995).[131] R. J. Radtke, V. N. Kostur, and K. Levin, Phys. Rev. B 53, R522 (1996).[132] M. Turlakov and A. J. Leggett, Phys. Rev. B 63, 064518 (2001).[133] R. Haenel, T. Tummuru, and M. Franz, Phys. Rev. B 106, 104505 (2022).[134] K. Clark, A. Hassanien, S. Khan, K.-F. Braun, H. Tanaka, and S.-W. Hla, Nat.Nanotechnol. 5, 261 (2010).[135] P. J. Hirschfeld, M. M. Korshunov, and I. I. Mazin, Rep. Prog. Phys. 74, 124508(2011).61Bibliography[136] B. White, J. Thompson, and M. Maple, Physica C 514, 246 (2015).[137] W.-Y. He, B. T. Zhou, J. J. He, N. F. Q. Yuan, T. Zhang, and K. T. Law, Commun.Phys. 1 (2018).[138] D. Shaffer, J. Kang, F. J. Burnell, and R. M. Fernandes, Phys. Rev. B 101, 224503(2020).[139] A. Mercado, S. Sahoo, and M. Franz, Phys. Rev. Lett. 128, 137002 (2022).[140] X.-Y. Song, Y.-H. Zhang, and A. Vishwanath, Phys. Rev. B 105, L201102 (2022).[141] S. Ryu and Y. Hatsugai, Phys. Rev. Lett. 89 (2002).[142] S. Uebelacker and C. Honerkamp, Phys. Rev. B 84 (2011).[143] K. Sun, H. Yao, E. Fradkin, and S. A. Kivelson, Phys. Rev. Lett. 103 (2009).[144] T. Tummuru, S. Plugge, and M. Franz, Phys. Rev. B 105, 064501 (2022).[145] R. S. Markiewicz, B. Singh, C. Lane, and A. Bansil, “High-order van hove singularitiesin cuprates and related high-tc superconductors,” (2021), arXiv:2105.04546 .[146] T. Wang, N. F. Q. Yuan, and L. Fu, Phys. Rev. X 11, 021024 (2021).[147] Y.-T. Hsu, F. Wu, and S. Das Sarma, Phys. Rev. B 104, 195134 (2021).[148] X. Lu and D. Se´ne´chal, “Doping phase diagram of a hubbard model for twistedbilayer cuprates,” (2021), arXiv:2112.00487 .[149] P. A. Volkov, S. Y. F. Zhao, N. Poccia, X. Cui, P. Kim, and J. H. Pixley, “Josephsoneffects in twisted nodal superconductors,” (2021), arXiv:2108.13456 .[150] O. Can, X.-X. Zhang, C. Kallin, and M. Franz, Phys. Rev. Lett. 127, 157001 (2021).[151] B. Zinkl, K. Hamamoto, and M. Sigrist, “Symmetry conditions for the supercon-ducting diode effect in chiral superconductors,” (2021), arXiv:2111.05340 .[152] F. Wu, R.-X. Zhang, and S. Das Sarma, Phys. Rev. Research 2, 022010 (2020).[153] T. Cea, N. R. Walet, and F. Guinea, Nano Letters 19, 8683 (2019).[154] L. Xian, A. Fischer, M. Claassen, J. Zhang, A. Rubio, and D. M. Kennes, “Engi-neering three dimensional moire´ flat bands,” (2020), arXiv:2012.09649 .[155] F. Schindler, A. M. Cook, M. G. Vergniory, Z. Wang, S. S. P. Parkin, B. A. Bernevig,and T. Neupert, Science Advances 4 (2018).[156] L. Fu, Phys. Rev. Lett. 106 (2011).[157] A. Cook and M. Franz, Phys. Rev. B 84 (2011).62[158] J. Xia, Y. Maeno, P. T. Beyersdorf, M. M. Fejer, and A. Kapitulnik, Phys. Rev.Lett. 97, 167002 (2006).[159] A. Kapitulnik, J. Xia, E. Schemm, and A. Palevski, New J. Phys. 11, 055060 (2009).[160] W. Harrison, Electronic Structure and the Properties of Solids: The Physics of theChemical Bond, Dover Books on Physics (Dover Publications, 2012).[161] A. Bansil, M. Lindroos, S. Sahrakorpi, and R. S. Markiewicz, Phys. Rev. B 71,012503 (2005).[162] J. A. Slezak, J. Lee, M. Wang, K. McElroy, K. Fujita, B. M. Andersen, P. J. Hirschfeld,H. Eisaki, S. Uchida, and J. C. Davis, Proceedings of the National Academy ofSciences 105, 3203 (2008).[163] P. Coleman, Introduction to many-body physics (Cambridge University Press, 2015).63Appendix AAppendices for Chapter 2A.1 Chern numberTo reliably determine the topology of a band structure, we frequently turn to the Chernnumber, which is defined as an integral of the Berry curvature over the full Brillouin zone:C = 12π∫BZΩ(k) d2k. Computation of the Berry curvature directly from the eigenstates iscomplicated by the arbitrary phases that accompany numerical diagonalization. To counterthis, we rely on a gauge independent formulation where [70]Ω(k) =∑m ̸=n⟨n|∇khk|m⟩ × ⟨m|∇khk|n⟩(Em − En)2 . (A.1)Here |n⟩ and |m⟩ are eigenstates of the Bloch Hamiltonian hk with energies En and Emrespectively.A.2 Construction from nanowiresIn this section, we construct a more realistic model of the twisted bilayer structure thattakes full account of the electron spin degree of freedom. We use the now standarddescription of a proximitized semiconductor nanowire as a 1D electron gas with strong spinorbit coupling (SOC) and induced superconducting order [25–27, 53–56]. When placed inan external magnetic field, such systems are known to realize the physics of Kitaev chains.We find that a coupled array of such wires forms a 2D p-wave superconductor and, in atwisted bilayer geometry, there is a range of parameters that supports the chiral px ± ipyphase.In the normal state, each standalone layer (labelled by a = 1, 2) that is composed ofweakly coupled parallel nanowires can be described by the Hamiltonianh0a(k) = ξka +B · σ + α(Ea × k) · σ, (A.2)where ξka is the electronic dispersion as in Eq. (2.5), B = (Bx, By, 0) is the in-planeexternal magnetic field, σ = (σx, σy, σz) is the vector of Pauli matrices in spin space and αcharacterizes the strength of SOC that originates from the internal electric field Ea. Withthe structure of SOC given above, one can account for both Rashba and Dresselhaus effects.In the most general form, SOC for electrons in the x-y plane (kz = 0) can be written ashSOC = α [Ez(σykx − σxky) + σz(Exky − Eykx)] . (A.3)For concreteness, we only retain the σz term in the above expression. If the magneticfield B is chosen to lie in the x-y plane then the SOC term will anti-commute with64A.2. Construction from nanowires0.30.20.10.00.10.20.3E/tlayer 1layer 20kx0.30.20.10.00.10.20.3E/tlayers 1&2Figure A.1: Energy spectra of the 8-band model projected onto kx momentum component.Top panel shows the two layers decoupled (g = 0), while in the bottom panel g = 0.3. Themodel parameters are t′ = 0.3, ∆1 = −i∆2 = 0.4, µ = −1.9, α = 0.5 andB = (0.35, 0.35, 0).Effective electric fields for SOC are chosen to be E1 = (0, 1, 0) and E2 = (1, 0, 0), suchthat SOC has a dispersion only along the wires. The phase difference between ∆1 and ∆2results in a gapped spectrum (≈ 0.035) with a Chern number C = 1.hB = B · σ, as in the 1D nanowire model, resulting in a gap at k = 0. More generally, asmall Ez can be tolerated and will not change the physics as long as it does not cause gapclosing.To obtain a 2D p-wave superconductor, we introduce s-wave superconductivity ineach layer and work with the Nambu spinor Ψa = (ck↑a, ck↓a, c†−k↓a,−c†−k↑a)T , where ckσaannihilates a fermion with momentum k and spin σ in layer a. In this basis, the BdGHamiltonian of a single layer isha =(h0a(k) ∆a∆∗a −σyh0a(−k)∗σy), (A.4)where ∆a = (∆′a + i∆′′a) denotes the singlet pairing SC order parameter in layer a. Itshould be viewed as a 2× 2 matrix in spin space.In the case where SOC is only along the wires, we have E1 = (0, Ey, 0) and E2 =(Ex, 0, 0) and the Hamiltonians of the two layers areh1 = [ξk1 − αEy sin(kx)σz]τz + hB +∆′1τx +∆′′1τyh2 = [ξk2 + αEx sin(ky)σz]τz + hB +∆′2τx +∆′′2τy,where (τx, τy, τz) are Pauli matrices in Nambu space and the linear dependence of SOChas been regularized with the substitution kx(y) → sin kx(y). This is justified because weintend to capture the physics near k = 0 in each wire.65A.2. Construction from nanowires0.30.20.10.00.10.20.3E/tlayer 1layer 20kx0.30.20.10.00.10.20.3E/tlayers 1&2Figure A.2: The same as Fig. A.1 except that the effective electric fields for SOC areE1 = (0.3, 1, 0) and E2 = (1, 0.3, 0). As we adiabatically turn on SOC perpendicular tothe nanowires, we find that the gap does not close and it is in fact larger (≈ 0.039), ascompared to that in Fig. A.1. The system remains in the topological phase C = 1.The two layers may now be coupled via nearest-neighbor interlayer tunnelings of theformh12 = g∑σ[c†kσ1ckσ2 + h.c.], (A.5)which results in a 8× 8 BdG Hamiltonian. We find that this model supports qualitativelysimilar phenomenology as the spinless fermion model discussed in the main text. Theparameter space for the present model is much larger, so in the following we only highlightthe most interesting phase characterized by Chern number C = 1.Fig. A.1 shows the spectrum for the case where SOC is dispersing only along the wires.We find that it is possible to obtain a stable px + ipy phase with C = 1 for a range ofparameters. As in the spinless fermion model, this occurs when the chemical potential µ istuned such that the underlying normal-state dispersion defines a single electron-like Fermisurface.We also consider a more general case (see Fig. A.2) where SOC has a componentthat disperses in a direction perpendicular to the nanowires, as determined by Eq. (A.3).Interestingly, as we turn on this effect, the system remains in the C = 1 topological phase.An increase in µ brings the system to phases with C = 2 and C = 0, as before. Thus, weconclude that the structure of the phase diagram deduced from the simple spinless fermionmodel in the main text is likely quite robust and remains valid for systems composed ofquantum wires.66A.3. Details of the vortex calculationA.3 Details of the vortex calculationIn order to simulate a vortex we generalize Hamiltonian (2.3) to include spatial dependenceof the order parameter by replacing ∆→ ∆rδ. We take ∆rδ as indicated in Eq. (2.6) andwrite the Hamiltonian in a matrix representationH =∑r,r′Ψ†rhrr′Ψr′ , (A.6)where Ψr = (ψr1, ψ†r1, ψr2, ψ†r2)T and ψra represents a column vector composed of operatorscr,a. For an L× L system hrr′ is a Hermitian matrix of dimension 4L2. From it, one canobtain the energy eigenvalues and eigenfunctions of the BdG quasiparticles, which we plotin Fig. 2.3.We make a plausible assumption that the two states with energies near zero are linearcombinations of Majorana zero modes localized in the vortex core and at the edge anddenote the corresponding Majorana operators as γc and γe, respectively. The energysplitting due to the small overlap in their respective wavefuctions can then be modeled byan effective low-energy Hamiltonianheff = iϵ0γcγe =12ϵ0Γ†σyΓ, (A.7)where Γ = (γc, γe)T and σy is the Pauli matrix. We note that the form of heff is uniquelydetermined by the requirement of Hermiticity and the self-adjoint property γ† = γ ofMajorana operators. Eq. (A.7) implies that the wavefunctions ϕ± belonging to eigenvalues±ϵ0 are eigenstates of σy, that is, they are related to the two Majorana wavefunctionsϕc and ϕe through the relation ϕ± = (ϕc ± iϕe)/√2. A numerical diagonalization givesus ϕ± and it is a simple matter to invert the above relation and obtain the Majoranawavefunctions ϕc and ϕe. These are plotted as insets in Fig. 2.3. The fact that theyare clearly localized near the vortex core and the edge, respectively, validates our initialassumption.A more complete approach to the vortex problem would include the external magneticfield to stabilize the vortex solution and a self-consistently determined order parameterprofile ∆rδ. As we are mainly interested in qualitative features of the system we do notpursue this line of inquiry here. We note, however, that these additional ingredients cannotaffect the existence of the Majorana zero mode (although they could change details such asthe spatial profile of its wavefunction, for instance). Once an isolated Majorana zero modeis established to exist in the vortex core it cannot be removed by any small perturbation,this being the essence of topological protection. The only way to remove it would be tobring the vortex to the edge of the system (where another zero mode resides) or merge itwith another vortex.67Appendix BAppendices for Chapter 3B.1 Microscopic derivation of GL coefficientsFrom the continuum BdG model, defined by the matrix Hamiltonian (3.10), it is possibleto derive the corresponding GL theory by directly expanding the free energy FBdG givenin Eq. (3.11) in the powers of the superconducting order parameters. This procedure ismade feasible by the fact that explicit expressions for the two positive-energy eigenvaluesof hk can be obtained and are given byEkα =√(∆2k1 +∆2k2)/2 + ξ2k + g2 + (−)αDk (B.1)where D2k = (∆2k1 −∆2k2)2/4 + g2(∆2k1 +∆2k2 + 4ξ2k)− 2g2∆k1∆k2 cosφ and α = 1, 2. Wenow express the two d-wave gap functions as∆k1 = ψ cos(2αk), ∆k2 = ψ cos(2αk − 2θ), (B.2)substitute into (B.1) and expand the free energy (3.11) in the powers of the order parameteramplitude ψ to fourth order. The expansion results in many terms and is best carried outwith the help of Mathematica or similar software capable of symbolic manipulations.We are primarily interested in terms governing the phase difference between the twolayers, that is, coefficients B and C in the GL free energy (3.2). These can be isolated byfocusing on terms containing powers of cosφ. Coefficient B is thus found as a prefactor ofψ2 cosφ in the above expansion and readsB = cos (2θ)NF∑s=±∫ ϵc0dξ−πsgξ|ξ − sg|tanh|ξ − sg|2kBT. (B.3)In deriving this expression we replaced the crystal momentum sum using the usual prescrip-tion∑k → NF∫ ϵc−ϵc dξ∫ 2π0 dϕ where NF = L2m/ℏ2 is the normal-state density of states atthe Fermi level, L is the system size and the angular integral has been performed. Theremaining integral over the energy variable cannot be completed in the closed form butone can show that the integrand fB(ξ) is negative for all values of ξ and thus the integralitself is also negative, as expected on physical grounds. This is illustrated in Fig. B.1(a).Coefficient C can likewise be identified as a prefactor of the ψ4 cos 2φ term in theexpansion and has the form C = (2 + cos 4θ)IC . Here IC = NF∫ ϵc0 fC(ξ)dξ is thecorresponding energy integral similar to Eq. (B.3). Its integrand fC(ξ) has a morecomplicated structure and for the sake of brevity we do not show it explicitly. ImportantlyfC(ξ) is positive for all ξ and, as anticipated in the main text, is approximately equalto the integrand of Eq. (B.3) squared, as shown in Fig. B.1(b). We thus conclude that68B.2. Self-consistent lattice mean-field theoryOut[!]= fB(ξ)0 1 2 3 4 5-5-4-3-2-10Out[!]=fC(ξ)fB2(ξ)0 1 2 3 4 50.00.20.40.60.81.0Out[!]= fB(ξ)0 1 2 3 4 5-5-4-3-2-10Out[!]=fC(ξ)fB2(ξ)0 1 2 3 4 50.00.20.40.60.81.0a) b) ξ/g ξ/gFigure B.1: Analysis of energy integrands. (a) Integrand fB(ξ) entering the expressionEq. (B.3) for the GL theory coefficient B. (b) Integrand fC(ξ) entering the expression forthe GL coefficient C compared to fB(ξ)2. Note that the amplitude of the latter is scaledby a constant factor of 1/30 for easier comparison. The temperature is chosen such thatkBT = 0.2g but the qualitative properties of the integrands remain the same in a widerange of temperatures.the continuum BdG model indeed supports the GL phenomenology deduced from generalphysical arguments. The calculation explicitly shows that prefactor K of the cos 2φ termin Eq. (3.4) is positive which we argued was the key condition underlying the emergence ofthe T -broken topological phase in the twisted dSC bilayers.B.2 Self-consistent lattice mean-field theoryTo describe the two-layer system microscopically, we employ the Hubbard Hamiltonian(3.12) and treat the interaction term at the mean-field level to account for d-wave super-conductivity in each layer. The MF Hamiltonian becomesH = −t∑⟨ij⟩σac†iσacjσa − t′∑⟨⟨ij⟩⟩σac†iσacjσa − µ∑iσaniσa+∑⟨ij⟩a(∆ij,ac†i↑ac†j↓a + h.c.)−∑ijσgijc†iσ1cjσ2, (B.4)where t and t′ are the hopping amplitudes between first and second neighbour sites onthe square lattice, µ is the chemical potential that controls on-site particle density niσaand ∆ij,a = V ⟨ci↑acj↓a⟩ is the complex order parameter. It is believed that the interlayertunneling processes in cuprates are mediated, to leading order, by the spherically symmetricCu 4s orbitals [94], although other orbitals such as those belonging to Bi and Sr atoms,probably also play an important role. With this in mind, we adopt a phenomenologicalform for the interlayer tunneling amplitude that decays exponentially with distance:gij = g0e−(rij−c)/ρ, (B.5)69B.2. Self-consistent lattice mean-field theorya) b) Figure B.2: Geometry of commensurate unit cells (shaded area) for bilayers with relativetwists of (a) θ1,2 ≃ 53.13o and (b) θ2,5 ≃ 43.60o.where r2ij = c2+ d2ij , with dij denoting the in-plane separation between two sites i and j ondifferent layers and c being the interlayer distance. In practice, we limit these terms to afew unit cells, beyond which the amplitudes are negligible. The characteristic radial extentof the Cu 4s orbital has been denoted by ρ. We have also considered interlayer couplings ofthe form gij = g0c2/r2ij [160]. Self-consistent calculations still result in topological phases,but this form of coupling, when compared to (B.5) for same the g0, effectively correspondsto stronger interlayer coupling as the decay of 1/r2 is slower compared to the exponential.While the above real space Hamiltonian is valid for any twist angle, the bilayer forms acrystal only for a discrete set of commensurate rotation angles. To carry out the latticemean-field calculations we focus on commensurate twist geometries that are amenable tonumerical study using the standard Bloch representation. We note that any commensuratetwist can be described by an integer-valued ‘twist’ vector v = (m,n). The correspondingmoire unit cell comprises of 2q = 2(m2 + n2) sites and the twist angle is determined by therelationθm,n = 2arctan (m/n). (B.6)This can be understood as a rotation of two perfectly aligned square lattices in oppositedirections by θm,n/2, which leads to the sites originally at the locations v and v¯ = (−m,n)to lie atop each other. Exemplary unit cells for commensurate twist angles are illustratedin Fig. B.2.To model the Fermi surface of Bi2212 we follow Ref. [90] and choose t = 153 meV,t′ = −0.45t with µ = −1.35t near optimal doping. Other parameters that determine theinterlayer tunneling amplitudes are as follows: the average spatial extent of electrons inthe 4s orbitals is ρ ≈ 4a.u. = 2.11A˚, the typical inter-layer distance is c ≈ 12A˚ [161] andthe in-plane separation between the Cu atoms is a ≈ 5.4A˚ [162].In the remainder of this section we outline the procedure for obtaining gap equationsfor the order parameters ∆ij,a using the imaginary time path integral formalism. Afterdecoupling the attractive pairing potential in the Cooper channel with the Hubbard-Stratonovich transformation, we obtain the action corresponding to the Hamiltonian (B.4)(cf. Ref. [163] for details)S =∑knΨ†kn[−iωn + hk]Ψkn +βNV∑i,j∈u.c.a∆ij,a∆∗ij,a. (B.7)70B.2. Self-consistent lattice mean-field theoryHere hk is the BdG Hamiltonian that follows from Fourier transforming Eq. (B.4), β =1/kBT is the inverse temperature, ωn = (2n+ 1)π/β are the Matsubara frequencies, V isthe strength of attractive nearest neighbour interactions, N is the number of unit cells andΨk = (Φk1,Φk2)T is the Nambu spinor such that Φka = (c(1)k↑a, . . . , c(q)k↑a, c(1)†−k↓a, . . . , c(q)†−k↓a)T ,with the superscripts 1, . . . , q labelling the sublattice degrees of freedom in layer a.The effective action can be determined by integrating out fermionic degrees of freedom:Seff = −∑knTr log[−G(k, iωn)−1] + βNV∑i,j∈u.c.a∆ij,a∆∗ij,a, (B.8)where we defined the Matsubara Green’s function G(k, iωn) = −(−iωn+hk)−1. To evaluatethe saddle point condition for this effective action, ∂Seff/∂∆∗ij,a = 0, we employ the identity∂∆(Tr logA) = Tr[∂∆AA−1]. After performing the Matsubara sum, this simplifies to∆ij,a = −VN∑kTr[∂hk∂∆∗ij,aUknF (Ek)U†k]. (B.9)In the above, Uk is the unitary matrix that diagonalizes the Bloch Hamiltonian as U†khkUk =Ek and nF (Ek) is a diagonal matrix where Fermi function is applied to the eigenvaluesEk. For decoupled layers (g0 = 0), with the tight binding parameters given above, aself-consistent solution converges to a d-wave form with maximal gap ∆max = 40 meV,when V = 146 meV. We use this value for calculations where two layers are coupled(g0 ≠ 0). Fig. 3.2 in the main text summarizes the results obtained from these calculations.71Appendix CAppendices for Chapter 4C.1 Interlayer supercurrent analysisAdditional insights into the temperature dependence of the interlayer supercurrent discussedin Sec. 4.2 can be gained by analyzing the Josephson relation in more detail. To this end itis useful to write down an explicit expression for the two positive eigenvalues of the BdGHamiltonian (4.7):Ek± =√(∆2k1 +∆2k2)/2 + ξ2k + g2 ±Dk(φ) (C.1)where D2k(φ) = (∆2k1 −∆2k2)2/4 + g2(∆2k1 +∆2k2 + 4ξ2k − 2∆k1∆k2 cosφ). Noting that thephase only enters through the cosine term in D2k(φ) it is possible, with use of Eq. (4.9), toexpress the supercurrent asJ(φ) = − sinφeg22ℏ∑k∆k1∆k2Dk(φ)∑a=±aEkatanh12βEka. (C.2)The superfluid stiffness follows from Eq. (4.4),ρs =eg22ℏ∑k∆k1∆k2Dk(0)∑a=±[ −aEkatanh12βEka]φ=0. (C.3)Noting that Ek+ > Ek− it is easy to see that the last term∑a[. . . ] in the above equationis positive for all β which justifies the expression given in Eq. (4.14).By linearizing the quasiparticle dispersion near the Dirac points – performing the ‘nodalapproximation’ – it is possible to extract from Eq. (C.3) the functional form of the low-Tbehavior,ρs(T ) ≃ ρs(0) + aθT − bθT 3 (C.4)with non-negative θ-dependent coefficients aθ and bθ. This expression is valid for twist anglesaway from θM and indicates that for aθ > 0 the ρs(T ) initially grows with temperature,as already observed in Figs. 4.3 and 4.4. The same calculation indicates a logarithmiccontribution ∼ ln (kBT/g) when θ = θM. This is also evident in Fig. 4.3(b).The analysis of the critical current is more complicated because J(φ) generally attainsits maximum at some generic phase angle φc which must be first determined. This ingeneral can only be done numerically. Some analytical progress can be made by notingthat at least for small twist angles φc tends to be close to π/2. One may thus approximateJc ≃ J(π/2) and then estimate J(π/2) using the nodal approximation as before. Onefinds exponentially activated behavior at the lowest temperatures reflecting the spectral72C.2. Expansion in powers of ggap that develops in the system at the large external phase bias. At somewhat highertemperatures J(π/2) behaves as in Eq. (C.4) but the log divergence near θM is now absent.This behavior is indeed observed in Fig. 4.3(a).C.2 Expansion in powers of gWe supply here some details leading to estimates quoted in Eqs. (4.15, 4.16). To estimateJc near zero twist angle we evaluate the current given in Eq. (C.2) at θ = 0 to leading orderin interlayer coupling g. In this regime, the current is dominated by ordinary Josephsontunneling, which is a g2 process, and we can thus set g = 0 inside the momentum sum. AtT = 0 this givesJ(φ) = sinφeg2ℏ∑k∆2 cos2 2αk2(ξ2k +∆2 cos2 2αk)3/2. (C.5)With this approximation, the maximum current occurs at φ = π/2 and thus Jc = J(π/2).We replace the momentum sum with an integral (as described in App. C.3) to obtainJc ≃ eg2ℏNF∫ 2π0dα∫ xc0dxcos2 2α(x2 + cos2 2α)3/2(C.6)with xc = ϵc/∆. In the physically relevant limit xc →∞, the integrals can be evaluatedanalytically which leads to Eq. (4.15) with C2 = 1.The same procedure applied at twist angle θ = 45◦ gives a vanishing current, consistentwith the notion that single Cooper pair tunneling is disallowed in this limit. To obtaina meaningful estimate for Jc one thus needs to carry the expansion to order g4 whichcaptures double-pair tunneling processes. Working once again at T = 0 we findJ(φ) = sinφeg2ℏ∑k∆2 sin 4αk2(Ek+ + Ek−)Ek+Ek−. (C.7)Expanding the denominator to order g2, collecting all non-vanishing terms, and notingthat in this case the current maximum occurs at φ = π/4 one arrives at Eq. (4.16) withthe constant given byC4 =∫ 2π0dα2πsin2 4α∫ xc0dx3(ϵ+ + ϵ−)2 − 2− 4x2(ϵ+ + ϵ−)3(ϵ+ϵ−)3, (C.8)where ϵ± =√2x2 + 1± cos 4α. Numerical integration in the limit xc →∞ gives C4 ≃ 0.55.C.3 In-plane superfluid stiffnessTo obtain the γ coefficient in the single-layer free energy Eq. (4.22) we begin from therelevant real-space BdG Hamiltonian,h =(h0 ∆(r)∆(r)∗ −h0)(C.9)73C.3. In-plane superfluid stiffness0.0 0.2 0.4 0.6 0.8 1.0 1.20.00.20.40.60.81.0T/TcΔ(T )Δ(0)γ(T )γ(0)Figure C.1: Temperature dependence of the superconducting gap amplitude (blue) calcu-lated by numerically iterating the gap equation (4.11), and the in-plane superfluid stiffnessγ(T ) (red) evaluated as the coefficient of the q2 term in Eq. (C.15).where h0 = −ℏ2∇2/2m − µ and ∆(r) = e2iqx∆d. The exponential factor in the orderparameter sets up a uniform superflow along the x direction with the amplitude proportionalto q. The coefficient γ measures the free energy cost of this supercurrent per unitvolume. To proceed it is useful to pass to a new gauge, taking c↑r → e−iqxc↑r andc†↓r → eiqxc†↓r. This removes the phase factor from ∆(r) and transforms the kinetic termas ±h0 → ∓(ℏ2/2m)(∇∓ iq)2 where q = (q, 0, 0). The advantage of this gauge is that theHamiltonian becomes translation invariant and can hence be expressed in the momentumspace ashk =(ξk−q ∆k∆∗k −ξk+q), (C.10)where ξk = ℏ2k2/2m− µ and ∆k = ∆cos(2αk). The standard expression for free energy(4.9) can now be used with the energy eigenvaluesEk± =ξk+q − ξk−q2±√(ξk+q + ξk−q2)2+∆2k. (C.11)The quantity of interest is the in-plane phase stiffness which follows from the γ termin Eq. (4.22) upon taking ψ(r) = ψ0eiφ(r). For the uniform supercurrent defined byφ(r) = 2qx this takes the form γ|∇ψ|2 = 4γψ20q2. Therefore, for a weak supercurrent, weseek the coefficient of q2 in the expansion of free energy in powers of q. Before we proceedwith the expansion it is useful to switch perspective slightly and view Eqs. (C.10-C.11)as describing a lattice model; specifically we may regard ξk and ∆k as lattice dispersionand gap function, respectively, defined in the Brillouin zone appropriate for the squareCuO2 lattice. This point of view avoids difficulties down the road where the continuummodel would show various ultraviolet-divergent k-space integrals. By contrast no suchdivergences appear in the lattice model because the BZ provides a natural UV cutoff.74C.3. In-plane superfluid stiffnessExpanding the free energy to second order in q we findF ≃ F0 − 12q2∑k[(∂2ξk∂k2x)(1− ξkϵktanhβϵk2)− β2(∂ξk∂kx)2sech2βϵk2], (C.12)where ϵk =√ξ2k +∆2k is the quasiparticle excitation energy. For a d-wave SC it is usefulto follow Ref. [120] and integrate the first term by parts. After some algebra we obtainF ≃ F0 + 12q2∑k(∂ξk∂kx)[∆2kϵ2k(∂ξk∂kx)− ∆kξkϵ2k(∂∆k∂kx)][1ϵktanhβϵk2− β2sech2βϵk2].(C.13)This expression, while seemingly more complicated than Eq. (C.12), has several desirablefeatures that simplify numerical evaluation. First, the coefficient of q2 explicitly vanisheswhen ∆k = 0; indeed one expects superfluid stiffness to be zero in the normal metal limit.Second, because of the powers of the excitation energy ϵk present in various denominators,it is clear that the largest contributions to the momentum sum come from the nodal regionswhere ϵk → 0. These properties facilitate a straightforward evaluation of the k-space sumand we can now safely revert to the continuum approximation. To this end, we considerfree energy density f = F/V , where V is the sample volume, and convert the k-space sumto an integral using the standard prescription1V∑k→ 1(2π)2∫d2k → 1(2π)2∫ 2π0dα∫kdk → NF∫ 2π0dα∫ ϵc−ϵcdξ. (C.14)Here NF = m/(2πℏ)2 is the density of states at the Fermi level and in the last stepwe restricted the integration to within a cutoff ϵc around the Fermi level which is nowa legitimate approximation because the integrand is strongly peaked here. With thesesimplifications the free energy density expansion can be written asf ≃ f0 + q2π2∫ π0dα(cosα cos 2α)2∫ ϵc0dξ∆2µϵ3[tanhβϵ2− βϵ2sech2βϵ2], (C.15)where ϵ =√ξ2 +∆2 cos2 2α. The integrals indicated in Eq. (C.15) lend themselves toa straightforward numerical evaluation and lead to the well-known curves for in-planesuperfluid stiffness of a d-wave SC with the characteristic T -linear dependence at lowtemperatures. An example of this behavior is given in Fig. C.1 along with the temperaturedependence of the gap function ∆(T ), which is required as input for this calculation.75Appendix DAppendices for Chapter 5D.1 Layer decoupling unitary transformationHere we derive a unitary transformation which, following the ideas in Ref. [38], reducesalternating twist multilayers to a combination of decoupled bilayers, with an additionalmonolayer for odd L. To begin, we reorganize the layer structure in the multilayerHamiltonian (5.1) according to an even/odd schemehk =ξkIo W ∆ok 0W † ξkIe 0 ∆ek(∆ok)∗0 −ξkIo −W0(∆ek)∗ −W † −ξkIe . (D.1)where Io (Ie) is an identity matrix of size no (ne) with no and ne the number of layerswith an odd and even index. While in the main text we take an isotropic form of ξk forconvenience, this is not necessary for the following derivation which holds as long as ξkis layer-independent. W is a matrix of size no × ne that accounts for arbitrary tunnelingamplitudes between neighboring layers:W =g12 0 0 0 . . .g23 g34 0 0 . . .0 g45 g56 0 . . .0 0 . . . . . . . . . . (D.2)The superconducting order parameters have been organized into diagonal matrices in theeven and odd layer subspaces as∆ok = diag(∆k,1,∆k,3, . . . ,∆k,2no−1,),∆ek = diag(∆k,2,∆k,4, . . . ,∆k,2ne .)(D.3)Given the 2D nature of superconductivity, weak interlayer tunneling is not expectedto significantly alter the order parameter amplitudes across the layers. Further, in analternating twist configuration, all layers in the odd (even) layer subspace are aligned. Theseobservations inform the assumption that order parameters in all the odd (even) numberedlayers can be treated as identical, that is, ∆ok = ∆okIo and ∆ek = ∆ekIe. Combinedwith the fact that the single-particle dispersion ξk is identical across the layers, thediscussion in Ref. [38] inspires the following trick: transform hk to a new basis usingV = diag (A,B,A,B), where A (B) is a unitary acting in the odd (even) layer subspace.76D.2. Low energy theory at magic anglesWe thus getV †hkV =ξkIo A†WB ∆okIo 0B†W †A ξkIe 0 ∆ekIe(∆ok)∗Io 0 −ξkIo −A†WB0 (∆ek)∗Ie −B†W †A −ξkIe . (D.4)The matrices A and B only act on the matrix of interlayer couplings W because all otherelements are ∼ I within the even and odd layer subspaces. We now choose A and B as inthe singular value decomposition of W , with W = AΛB† and Λ a no × ne diagonal matrixcontaining the ne real eigenvalues λm of√W †W . That is,V †hkV =ξkIo Λ ∆okIo 0Λ ξkIe 0 ∆ekIe(∆ok)∗Io 0 −ξkIo −Λ0 (∆ek)∗Ie −Λ −ξkIe . (D.5)This transformation has revealed a block structure with ne decoupled blocks describingbilayers with interlayer couplings λm. Note that the two ‘layers’ in these bilayer blocksdo not correspond to physical layers, but rather to a superposition of odd (even) layersdetermined by the A (B) matrices in the singular value decomposition. Furthermore,when L is odd (|no − ne| = 1), the transformation additionally gives rise to a decoupledmonolayer block characterized by a linear combination of the odd numbered layers.While the above derivation holds for a generic form of the coupling matrix W , itbecomes particularly simple when neighboring layers are coupled with identical strengthsW = g (δij + δi+1,j). The eigenvalues of√W †W are then given by [38]λm = g˜m = 2g cos(πmL+ 1)(D.6)with m = 1 . . . ne. For the cases L = 3 and L = 4 discussed in the main text, we haveg˜ = λ1 =√2g and (g˜1, g˜2) = (λ1, λ2) =(−1+√52 ,1+√52)g, respectively.D.2 Low energy theory at magic anglesTo arrive at a description of the multilayer systems near the magic angles we work in thelow energy limit where the Dirac cones in the spectrum are linearized. To elucidate thisfirst consider a monolayer d-wave superconductor. The Dirac cones are situated at theFermi momentum kF on the BZ diagonals, where the gap function vanishes. In the firstquadrant of the BZ, treating the node as the new origin and expanding in terms of smallmomenta q, one obtains the Dirac Hamiltonianhq = (v∆q⊥)τ1 + (vF q∥)τ3 (D.7)where τj denote Pauli matrices acting in Nambu space and vF (v∆) is the velocity in thedirection parallel (perpendicular) to the BZ diagonal.77D.2. Low energy theory at magic anglesQuadratic band touchingTo study magic angles in AT trilayers, the unitary transformation in Appendix D.1suggests that it suffices to work with a standalone bilayer with a renormalized interlayerlayer coupling g˜.For a moment consider an untwisted bilayer with no interlayer tunneling, i.e., g = 0. Insuch a scenario, Dirac cones from the two layers overlap with each other and lie on the BZdiagonal. Focusing again on the first quadrant, we treat this nodal point as our point ofreference and new origin q = 0. Now if the two layers are rotated slightly by ±θ/2, the newnode locations are given by ±QN = (0,±θkF /2) with respect to the node of an untwistedlayer. Turning on the interlayer tunneling, the full low energy Hamiltonian is given byh =(hq+QN tt hq−QN)(D.8)where the tunneling matrix t = g˜τ3. With the definition of a dimensionless parameterα = v∆kF θ/g˜, evaluating the analytical form of the dispersion is straightforward. Whilethe full expression is uninformative, we highlight its key features: When α > 2, the twoDirac nodes are positioned symmetrically off the diagonal. For α < 2 they both lie on thediagonal. Precisely at α = 2, however, the nodes merge and give rise to a quadratic bandtouching (QBT). Note that our definition of α differs from that in Ref. [48] by a factor of 2.An effective theory for small momenta near the QBT can be obtained using degenerateperturbation theory. Suppose that |α⟩ and |β⟩ are the degenerate zero-energy states of hprecisely at the node and |µ⟩ and |ν⟩ label the rest of the states. Up to third order, wethen have⟨β|heffq |α⟩ = ⟨β|h|α⟩+∑µ⟨β|h|µ⟩⟨µ|h|α⟩−Eµ +∑µ,ν⟨β|h|µ⟩⟨µ|h|ν⟩⟨ν|h|α⟩EµEν. (D.9)Keeping only the leading terms, which are quadratic in momentum, one arrives at theexpressions (5.7) and (5.8).Cubic band touchingFor a trilayer with a chiral twist, the Dirac cones of the middle layer always lie on the BZdiagonal. Treating the cone in the first quadrant as the reference point q = 0 as before,the other two nodal points are located at ±QN = (0,±θkF ). Accounting for the interlayertunneling, the full Hamiltonian readsh =hq+QN t 0t hq t0 t hq−QN (D.10)Unlike the bilayer, the dispersion now does not have a simple analytic form. Nevertheless,we find numerically that the Dirac cones move around as a function of α = v∆kF θ/g.Specifically, all three nodes lie on the normal state Fermi surface when α >√2, merge intoa cubic crossing at α =√2 and all three fall on the diagonal when α is decreased further.78D.3. Secondary instabilitiesUsing third-order perturbation theory following Eq. (D.9), the spectrum can be down-folded to retain just the two low-energy bands near the cubic crossing. The low-energyeffective Hamiltonian is defined by d1(q) and d3(q) in Eq. (5.10).D.3 Secondary instabilitiesIn addition to the primary d-wave pairing channel with strength V > 0, let us suppose thatan additional pairing channel with strength U > 0 exists in each layer. We take it to beisotropic (s-wave) for simplicity:Hs = −2UN∑kplc†k↑lc†−k↓lc−p↓lcp↑l. (D.11)Assuming U ≪ V, such that the d-wave component remains largely unaffected, we candefine the secondary s-wave order parameter as ∆sl = −2UN∑p⟨c−p↓lcp↑l⟩. The mean-fieldHamiltonian then reads H = E0+∑kΨ†khkΨk, where E0 =∑kl ξk+N∑l( |∆l|22V +|∆sl |22U).Correspondingly, the full order parameter of each layer is given by∆kl = ∆l cos(2αk + 2θl) + ∆sl , (D.12)with complex ∆sl . Recall that the free energy of the system at inverse temperature β, inthe BCS mean-field approach, readsF = E0 − 2β−1∑kαln [2 cosh(βEkα/2)] , (D.13)where Ekα are the positive-energy eigenvalues of the BdG Hamiltonian hk. The gapequation for the secondary order parameter is obtained by minimizing F with respect to∆s∗l :∆sl =2UN∑kαtanh(βEkα2)⟨kα| ∂hk∂∆s∗l|kα⟩. (D.14)In Fig. D.1 we show phase diagrams for trilayers in both the AT and CT stackings, inthe presence of the secondary instability. In our numerics we simultaneously solve for thes-wave and d-wave components, ∆sl and ∆l. We find that ∆s acquires a purely imaginaryvalue that results in a gapped spectrum, for both QBTs (as obtained in AT trilayers) andCBTs (as obtained in CT trilayers) at their respective magic angles. Crucially, we findthat the T breaking secondary instability is nucleated only for large U/V ∼ 0.3 in ournumerics, which can be understood in terms of an exponential (power-law) suppression ofthe induced gap for quadratic (cubic) band touchings, as we clarify in the following section.D.4 Density of states and gap scalingWorking in the low-energy limit, it is possible to derive simple analytical scaling formsfor the gap opening induced by residual interactions at the quadratic and cubic band79D.4. Density of states and gap scalingFigure D.1: Self-consistent mean-field results for the s-wave order parameter amplitude ∆s2in the middle layer (black lines) and the d-wave phase factor φ1 (red lines) as a functionof twist angle θ for trilayers in the AT (left) and CT (right) stacking. We consider asecondary s-wave pairing channel with relative strength U/V = 0.28, 0.32 and 0.36, anduse parameters g = 12meV and ϵc = 60meV, while the primary d-wave order parameter∆ ∼ 40meV. For U/V ≲ 0.30 we find that the d + is solution is not nucleated, and thed+ id′ solution first develops around the magic angles. For large values of U we observe acompetition between the two gapped T -broken phases, with the d+ id′ solution preferredat large angles and the d+ is solution preferred near the magic angles.touchings. Since self-consistent numerics from Appendix D.3 indicate that the T -breakings-wave order parameter is purely imaginary, the effective low energy Hamiltonian may bewritten asHeff = N2UM2 +∑q[d1(q)η1 +Mη2 + d3(q)η3], (D.15)where the q dependence of d1(q) and d3(q) is determined by the type of crossing asdescribed in the main text and Appendix D.2. Assuming that the solution in the d-wavechannel remains unchanged, one can minimize the effective free energy with respect to themass term as∂Feff∂M= NMU −∑q[tanh(βEq2)∂Eq∂M]= 0, (D.16)and thusM = U∫d2q(2π)2tanh(βEq2)∂Eq∂M, (D.17)where we transformed 1N∑q →∫d2q/(2π)2. We then express the integral in energydomain using∫d2q/(2π)2 → ∫ dE D(E) with the density of states D(E), and considerthe zero-temperature limit (β →∞) whereM = U∫dE D(E)∂E∂M. (D.18)Let us now focus on the two cases of interest.80D.4. Density of states and gap scalingQuadratic band touchingThe QBT, described at low energies by the effective Hamiltonian Eq. (5.8), consists of ananisotropic dispersion Eq =12g˜ (q2∥v2F +q2⊥v2∆). A generic cross section of this dispersion is anellipse, with major and minor axes given by a =√2g˜E/vF and b =√2g˜E/v∆ respectively.The density of states is easily obtained by a standard trick: let us consider the quasiparticledensity at energy E defined as ρ(E) =∫ E0 dE′D(E′), which is simply given by the areaenclosed by the ellipse at energy E normalized by the momentum spacing,ρ(E) = πab(12π)2=g˜E2πvF v∆. (D.19)Differentiating and taking into account the 4 QBTs in the Brillouin zone, we find the totalDOSD(E) =2g˜πvF v∆. (D.20)Using the derivative of the gapped spectrum Eq =√d1(q)2 + d3(q)2 +M2 with respectto the induced mass term ∂E/∂M =M/E, and noting that the integration over energyranges from the band bottom to a high-energy cutoff Λ≫M , this leads toπvF v∆2U g˜ = lnΛ− lnM, (D.21)which is just the BCS scaling formM = Λexp(−πvF v∆2U g˜). (D.22)This derivation shows that the induced gap M is exponentially small in the product of theattractive strength U and the effective interlayer tunneling g˜.Cubic band touchingThe CBT is described at low energies by the effective Hamiltonian Eq. (5.10) withanisotropic dispersion Eq =14g2(q2∥v2F + q2⊥v2∆)3/2. The cross section is again an el-lipse, but the major and minor axes scale differently with energy, a = (4g2E)1/3/vFand b = (4g2E)1/3/v∆ respectively. The number of states enclosed by the ellipse isρ(E) =ab4π=(4g2E)2/34πvF v∆, (D.23)which leads to the power-law divergent density of states (again counting the four QBTs inthe Brillouin zone)D(E) =2(2g)4/33πvF v∆E−1/3. (D.24)And since ∂E/∂M =M/E, we have3πvF v∆2(2g)4/3U =∫ ΛMdE1E4/3. (D.25)81D.5. Lattice modelThe above integral is UV convergent and evaluates to 3M−1/3 in the limit Λ/M →∞. Wethus findM =(2g)4(2U)3(πvF v∆)3 , (D.26)a power-law dependence of the T -breaking order parameter on the s-wave interactionstrength U . This power-law dependence is asymptotically stronger than the exponentialscaling for the QBT but is nevertheless strongly suppressed at small g due to the scalingM ∼ g4.D.5 Lattice modelFor the case of trilayers, we further use a twisted lattice model at commensurate twistangles to corroborate the findings of our continuum model analysis. This approach allowsus to compute the Chern number C of the T -broken phases, as well as probe the edgestates on an infinite cylinder geometry. Following the approach outlined in Ref. [47] weconsider a twisted square lattice Hubbard model with nearest neighbor density-densityelectron interactions. After a mean-field decoupling, one obtains a d-wave superconductorin each layer,H =− t∑⟨ij⟩σlc†iσlcjσl − t′∑⟨⟨ij⟩⟩σlc†iσlcjσl − µ∑iσlniσl+∑⟨ij⟩l(∆ij,lc†i↑lc†j↓l + h.c.)−∑ijσ⟨lm⟩gijc†iσlcjσm, (D.27)where l is a layer index, t and t′ are the nearest-neighbor and next-nearest-neighbor hoppingamplitudes on the square lattice, µ is the chemical potential that controls on-site particledensity niσl and ∆ij,l the complex order parameter on the bond connecting sites i and j.For the numerical results presented in the main text, we set t′ = 0 for simplicity, whichleads to a single electron-like Fermi surface centered around k = 0 for each monolayer.We also checked that our results are not qualitatively altered by modeling instead theFermi surface of near-optimally-doped Bi2212, which comprises four hole pockets centeredaround wavevectors (±π,±π), by setting the parameters t′ = −0.45t and µ = −1.35tfollowing Ref. [90]. For interlayer tunneling, we assume an exponentially decaying formgij = e−(rij−c)/ρ that is circularly symmetric and connects sites separated by rij ; the scaleis set by parameters c = 2.2 and ρ = 0.4 in units of the lattice constant [47].To make use of the Bloch representation of wavefunctions, one needs to work withcommensurate twist angles such that the multilayer system forms a crystallographic unitcell. For alternating twists, the commensurability condition, which is same as that forthe bilayer, can be specified in terms of a twist vector (m,n) with integers m,n. Thecorresponding twist angle is θm,n = 2arctan(m/n) and the number of sites in the unit cellis given by L(m2 + n2).In the chiral twist case, a periodic lattice forms between neighboring layers for the sameθm,n derived above; however, there is now the additional constraint that the non-adjacentlayers must also form a commensurate supercell. Let us consider a trilayer and denote82D.5. Lattice modelthe twist between layers 1-2 and 1-3 by vectors (m,n) and (p, q), respectively. We areinterested in the case where θp,q = 2θm,n, that is,arctan(p/q) = 2 arctan(m/n), (D.28)which, using the trigonometric expansion for tan(2θ), simplifies topq=2(m/n)1− (m/n)2 . (D.29)This implies that the supercell corresponding to the (m,n) chirally twist trilayer comprises3(p2+q2) lattice sites. When (m,n) = (1, 2), which corresponds to a twist angle θ1,2 = 36.9as in Fig. 5.2, (p, q) = (4, 3) and the smallest unit cell has 75 sites.This exercise can be extended to L = 4, for instance, by solving for the relationsθp,q = 2θm,n and θr,s = 3θm,n, where (r, s) now denotes the twist vector for layers 1-4.83
Thesis/Dissertation
eng
Physics
Vancouver : University of British Columbia Library
Attribution-NonCommercial-NoDerivatives 4.0 International
Graduate
Franz, Marcel
Topological superconductivity in twisted 2D structures
Text