# High-temperature topological superconductivity in twisted double layer
copper oxides

Oguzhan Can  Tarun Tummuru  Ryan P. Day  Ilya Elfimov  Andrea Damascelli
Marcel Franz  Department of Physics and Astronomy, University of British
Columbia, Vancouver BC, Canada V6T 1Z4  Quantum Matter Institute, University
of British Columbia, Vancouver BC, Canada V6T 1Z4

###### Abstract

A great variety of novel phenomena occur when two-dimensional materials, such
as graphene or transition metal dichalcogenides, are assembled into bilayers
with a twist between individual layers. As a new application of this paradigm,
we consider structures composed of two monolayer-thin $d$-wave superconductors
with a twist angle $\theta$ that can be realized by mechanically exfoliating
van der Waals-bonded high-$T_{c}$ copper oxide materials, such as
Bi2Sr2CaCu2O8+δ. On the basis of symmetry arguments and detailed microscopic
modelling, we predict that for a range of twist angles in the vicinity of
$45^{\rm o}$, such bilayers form a robust, fully gapped topological phase with
spontaneously broken time-reversal symmetry and protected chiral Majorana edge
modes. When $\theta\approx 45^{\rm o}$, the topological phase sets in at
temperatures close to the bulk $T_{c}\simeq 90$ K, thus furnishing a long
sought realization of a true high-temperature topological superconductor.

In a remarkable recent development Yu and co-workers [1] succeeded in
isolating and probing single monolayer samples of high-temperature
superconductor Bi2Sr2CaCu2O8+δ (Bi2212) with critical temperatures $T_{c}$
essentially indistinguishable from the bulk $T_{c}\simeq 90$ K of optimally
doped Bi2212 crystals. This result delivers the long-awaited proof that
high-$T_{c}$ superconductivity in Bi2212 (and possibly all the cuprates) is an
intrinsic property of a single two-dimensional monolayer. It also opens
exciting new possibilities for engineering structures with novel behaviour and
functionality, from individual high-$T_{c}$ cuprate monolayers. Motivated by
the results of Ref. [1] and recent developments in van der Waals-bonded
materials [2, 3, 4, 5, 6, 7], we consider here structures composed of two
cuprate monolayers with a twist angle $\theta$ between them.

Twisted double layer graphene famously exhibits a sequence of gate-voltage
controlled correlated insulating and superconducting phases [2, 3, 4, 5, 6]
when $\theta$ is close to the magic angle $1.1^{\rm o}$. Similar physics is
reported [7] in twisted bilayer WSe2 in the range of angles between $4.0^{\rm
o}$ and $5.1^{\rm o}$. In both cases, the observed behaviour is thought to
reflect a subtle interplay between the band structure energetics and topology
brought about by the Moiré patterns, electron-electron interactions and
disorder [8, 9, 10, 11, 12, 13, 14, 15, 16, 17]. As such, many aspects of the
reported experimental observations in twisted graphene and WSe2 still await
full theoretical understanding [18]. As we argue below, the problem of twisted
Bi2212 bilayers is potentially simpler to analyze theoretically, at least at
the level where each Bi2212 layer is treated as a BCS superconductor with a
$d_{x^{2}-y^{2}}$ order parameter. Results of such a theoretical analysis are,
nevertheless, non-trivial: over a continuous range of twist angles $\theta$ in
the vicinity of $45^{\rm o}$, we find that below a critical temperature
$\tilde{T}_{c}(\theta)$ the bilayer system spontaneously breaks time-reversal
symmetry ${\cal T}$ and forms a chiral topological phase that can be thought
of as an emergent $d_{x^{2}-y^{2}}\pm id_{xy}$ superconductor ($d\pm
id^{\prime}$ for short). Such a superconductor has a fully gapped bulk and is
characterized by Chern number $C=\pm 2$ per layer with a gap of
$d_{x^{2}-y^{2}}\pm id_{xy}$ symmetry, respectively. It exhibits two gapless,
topologically protected chiral edge modes per monolayer. These modes are
coherent superpositions of electron and hole states and can be therefore
classified as Majorana particles [19]. When $\theta\approx 45^{\rm o}$, our
results indicate that $\tilde{T}_{c}\simeq T_{c}$, suggesting that,
remarkably, in carefully assembled Bi2212 samples, the topological phase could
set in well above the liquid nitrogen temperature.

Figure 1: Schematic view of two copper-oxygen square lattices with twist angle
close to $45^{\rm o}$. Superconducting order parameter with a
$d_{x^{2}-y^{2}}$ symmetry in the twisted layer (red) resembles a $d_{xy}$
order parameter in the original untwisted coordinate frame (blue). Although
each Bi2212 monolayer contains two CuO2 planes, for simplicity and
concreteness we primarily focus on models with a single CuO2 plane per
monolayer, as would be the case in Bi2Sr2CuO6+δ and other high-$T_{c}$
compounds.

The physics underlying the formation of the $d\pm id^{\prime}$ topological
phase is illustrated in Fig. 1. In momentum space, the $d_{x^{2}-y^{2}}$ order
parameter that characterizes each monolayer can be represented as
$\Delta_{\bm{k}}=\Delta_{0}(\hat{k}_{x}^{2}-\hat{k}_{y}^{2})$, where
$\Delta_{0}$ denotes the amplitude and
$(\hat{k}_{x},\hat{k}_{y})=(k_{x},k_{y})/k$ are momentum components defined
relative to the principal crystal axes of the underlying copper-oxygen
lattice. When one of the layers is rotated by angle $\theta$, its order
parameter is transformed in the original coordinate frame to

$\Delta_{\bm{k}}^{(\theta)}=\Delta_{0}\left[\cos(2\theta)(\hat{k}_{x}^{2}-\hat{k}_{y}^{2})+\sin(2\theta)2\hat{k}_{x}\hat{k}_{y}\right].$
(1)

Figure 2: Twisted double layer $d$-wave superconductor. Panel (a) illustrates
the lattice geometry at small twist angle $\theta<\theta_{c}^{-}$ where the
free energy is minimized for interlayer phase difference $\varphi=0$. This
results in a gapless spectrum with weakly split Dirac points. For twist angles
$\theta>\theta_{c}^{-}$ in panel (b) the preferred state has $\varphi\neq 0$.
It breaks time reversal and is fully gapped as indicated by the massive Dirac
spectrum. The spectra are obtained by diagonalizing the BdG Hamiltonian (10);
$k_{\perp}$ and $k_{\parallel}$ denote momentum components perpendicular and
parallel to the Fermi surface near a nodal point. Panel (c) shows the GL free
energy ${\cal F}(\varphi)$ for twist angles $\theta$ ranging from $0^{\rm o}$
to $45^{\rm o}$ in $3^{\rm o}$ increments for ${\cal K}=0.125$. When
$\theta>\theta_{c}^{-}$ minima occur away from $\varphi=0$. Panel (d):
predicted phase diagram based on the the GL theory (I) and
$\tilde{T}_{c}(\theta)$ given by Eq. (6). TSC denotes the gapped topological
phase, $d$SC stands for gapless $d$-wave superconductor. Panel (e): Fermi
surface and $d$-wave gap in the continuum formulation of the microscopic
model. Panel (f) displays the relative phase between the layers $\varphi_{\rm
min}$ that minimizes the GL and the BdG free energy as a function of twist
angle $\theta$. Orange symbols represent the spectral gap.

From the vantage point of the unrotated monolayer, the second term in Eq. (1)
has the functional form of a $d_{xy}$ order parameter. For coupled layers, we
then expect the Cooper pair tunneling to induce a subdominant $d_{xy}$ order
parameter in the unrotated layer. By symmetry, the same should happen in the
other layer. Furthermore, if the two order parameters combine with a complex
phase, the result is a fully gapped topological superconductor with broken
${\cal T}$. The gapped nature of the phase follows from the fact that the
excitation gap is given by $|\Delta_{\bm{k}}|$ which, for a $d\pm id^{\prime}$
gap function, is non-vanishing everywhere at the Fermi level, Fig. 2(a,b).
${\cal T}$-breaking is evident because time reversal maps $d+id^{\prime}$ to
$d-id^{\prime}$. We demonstrate below that for a range of twist angles
$\theta$ between $\theta_{c}^{-}$ and $\theta_{c}^{+}$, such a ${\cal
T}$-broken phase becomes energetically favorable below a critical temperature
$\tilde{T}_{c}(\theta)$.

We note that ${\cal T}$-broken $d\pm id^{\prime}$ and $d\pm is$ phases have
been discussed theoretically in various contexts, including bulk high-$T_{c}$
cuprates in applied magnetic field [20, 21, 22], $c$-axis twist Josephson
junctions [23, 24, 25, 26], and highly doped single-layer graphene [27, 28].
However, conclusive experimental evidence for these states has never been
reported. A ${\cal T}$-broken topological phase with $p_{x}\pm ip_{y}$
symmetry has long been a leading candidate for superconducting order in
Sr2RuO4 [29, 30, 31]. Recent Knight shift measurements [32], however, appear
to rule out this possibility. By placing our emphasis on monolayer-thin
samples of an established $d_{x^{2}-y^{2}}$ superconductor, we argue that the
present work promises to finally break the logjam in the quest for a chiral
topological phase. As we show below, for $\simeq 45^{\rm o}$ twist, the ${\cal
T}$-breaking $d\pm id^{\prime}$ phase is nucleated robustly and reliably
through the simple and well understood mechanism of interlayer Cooper pair
tunneling that is largely insensitive to microscopic details.

##  I Ginzburg-Landau theory

The phenomenon outlined above is most easily quantified by examining the
phenomenological Ginzburg-Landau (GL) theory. For the sake of clarity and
simplicity, we begin by considering models with a single CuO2 plane per
monolayer. This is directly relevant to Bi2Sr2CuO6+δ (Bi2201), but, as we show
in the subsequent sections, our results are broadly applicable to cuprates
with multiple CuO2 plane structure, such as Bi2212, which has already been
exfoliated [1] to the monolayer limit.

The GL free energy density for two coupled monolayers can be written as

$\displaystyle{\cal F}[\psi_{1},\psi_{2}]$ $\displaystyle=$ $\displaystyle
f_{0}[\psi_{1}]+f_{0}[\psi_{2}]+A|\psi_{1}|^{2}|\psi_{2}|^{2}$
$\displaystyle+$ $\displaystyle B(\psi_{1}\psi_{2}^{\ast}+{\rm
c.c.})+C(\psi_{1}^{2}\psi_{2}^{\ast 2}+{\rm c.c.}),$

where $\psi_{a=1,2}$ are complex order parameters of the two layers and
$f_{0}[\psi]=\alpha|\psi|^{2}+{\frac{1}{2}}\beta|\psi|^{4}$ is the free energy
of a monolayer. We are interested in spatially uniform solutions, so the
gradient terms have been omitted. If the two layers are physically identical,
then the order parameter amplitudes must be the same and the most general
solution (up to an overall phase) is

$\psi_{1}=\psi,\ \ \ \psi_{2}=\psi e^{i\varphi},$ (3)

where we take $\psi$ to be real and positive.

Because $\psi_{a}$ describe order parameters with $d$-wave symmetry, they
transform as $\psi_{a}\to-\psi_{a}$ under a $90^{\rm o}$ rotation. Suppose we
increase the twist angle $\theta$ continuously from 0 to $\pi/2$. Since
$\psi_{2}$ changes sign under this operation, we conclude that for the free
energy (I) to remain invariant, the parameter $B$ must also change sign. The
simplest dependence consistent with this constraint is
$B=-B_{0}\cos(2\theta)$, where we take $B_{0}>0$. The negative sign is chosen
to allow the phase difference between the layers to vanish as the twist angle
approaches zero. This also implies that $B=0$ for $\theta=\pi/4$, consistent
with the intuition that a $d_{x^{2}-y^{2}}$ Cooper pair cannot tunnel between
the layers in this configuration. Parameters $A$ and $C$ are not required by
symmetry to depend on $\theta$, so we treat them as constants henceforth.
Incorporating these ingredients into Eq. (I), we find

${\cal F}(\varphi)={\cal
F}_{0}+2B_{0}\psi^{2}\left[-\cos(2\theta)\cos{\varphi}+{\cal
K}\cos({2\varphi)}\right],$ (4)

where ${\cal K}=C\psi^{2}/B_{0}$ and ${\cal F}_{0}$ contains all terms
independent of $\varphi$. It is now a simple matter to minimize ${\cal
F}(\varphi)$ for a given twist angle $\theta$. An interesting scenario arises
when ${\cal K}>0$, as indicated by Fig. 2(c). Although symmetry alone does not
constrain the sign of ${\cal K}$, a simple physical argument suggests that
${\cal K}$ should, in general, be positive. One can interpret the
$\psi_{1}^{2}\psi_{2}^{\ast 2}$ term in Eq. (I) as representing the coherent
tunneling of two Cooper pairs 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 pair tunneling.
Coefficient $C$ (and, hence, also ${\cal K}$) will then be positive,
regardless of the sign of $B$. Below, we demonstrate that microscopic models
fully support this argument, yielding ${\cal K}>0$.

At $\theta=45^{\rm o}$ and for ${\cal K}>0$, the free energy minimum clearly
occurs at $\varphi_{\rm min}=\pm\pi/2$ which indicates a $d\pm id^{\prime}$
superconducting state. For $\theta\neq 45^{\rm o}$, the free energy is
minimized by

$\varphi_{\rm min}=\arccos{\left(\frac{\cos{2\theta}}{4{\cal K}}\right)}.$ (5)

As illustrated in Fig. 2(f), a nontrivial phase difference occurs for a range
of angles $\theta_{c}^{-}<\theta<\theta_{c}^{+}$ with
$\theta_{c}^{\pm}=\frac{1}{2}\arccos{(\mp 4{\cal K})}$. This indicates a
${\cal T}$-broken phase with order parameter $d_{x^{2}-y^{2}}+e^{\pm
i\varphi_{\rm min}}d_{xy}$ which is fully gapped and topologically non-trivial
as long as $\varphi_{\rm min}\neq 0,\pi$, see also Fig. 2(a,b). As ${\cal K}$
depends on temperature through its dependence on $\psi$, the critical
temperature $\tilde{T}_{c}$ of the topological phase will be a function of the
twist angle $\theta$. If we adopt the standard GL temperature dependence for
the order parameter $\psi(T)=\psi_{0}\sqrt{1-T/T_{c}}$, it is straightforward
to deduce

$\tilde{T}_{c}(\theta)=T_{c}\left(1-\frac{|\cos{2\theta}|}{4{\cal
K}_{0}}\right),\ \ \ \theta_{c}^{-}<\theta<\theta_{c}^{+},$ (6)

where ${\cal K}_{0}=C\psi_{0}^{2}/B_{0}$. This defines the phase diagram in
Fig. 2(d). We observe that, remarkably, $\tilde{T}_{c}(45^{\rm o})$ coincides
with the $T_{c}$ of a single monolayer which can be as high as 90 K in
carefully prepared Bi2212 flakes. Away from the ideal $45^{\rm o}$ twist, the
critical temperature falls approximately as $\tilde{T}_{c}(\theta)\simeq
T_{c}(1-|\theta-\pi/4|/2{\cal K}_{0})$. Microscopic models with realistic
Bi2212 parameters discussed in the following give typical values ${\cal
K}_{0}\simeq 0.1-0.2$, implying a significant extent for the topological
phase, as indicated in Fig. 2(d).

##  II Microscopic models

We now turn to the microscopic theory. Although the pairing mechanism in
high-$T_{c}$ cuprates remains a subject of debate, it is widely accepted that
most physical properties of the superconducting state in the optimally doped
and overdoped regime are accurately described within the framework of the
standard BCS theory with a $d_{x^{2}-y^{2}}$ order parameter. Therefore, we
begin by modeling the twisted bilayer using a simple continuum model of
coupled $d$-wave superconductors, and complement this with a calculation based
on an attractive Hubbard model on the square lattice. The continuum
formulation has the advantage of being applicable to an arbitrary twist angle
$\theta$, similar to the Bistritzer-MacDonald model for graphene [8]. In
contrast with the continuum model above, the lattice calculation captures
various microscopic details of CuO2 planes more effectively, such as the shape
of the Fermi surface. It is, however, limited to commensurate twist angles
that produce relatively small Moiré unit cells due to practical constraints.

Figure 3: Lattice model results. Panels (a,b) show the temperature dependence
of the minimum gap $\Delta_{\rm min}$, the maximum gap $\Delta_{\rm max}$ and
phase $\varphi$, based on a fully self-consistent lattice calculation for
coupled layers with commensurate twist angles $\theta_{1,2}\simeq 53.13^{\rm
o}$, and $\theta_{2,5}\simeq 43.60^{\rm o}$ corresponding to a unit cell with
10 and 58 sites, respectively. Panels (c-e) show zero temperature phase
diagrams of the system for three twist angles as a function of chemical
potential $\mu$ and interlayer coupling $g_{0}$. The range of chemical
potentials $\mu\in(-1.6t,-1.0t)$ correspond to range of fillings
$n\in(0.033,0.04)$ near optimal doping. Each data point in the phase diagram
corresponds to an independent self-consistent solution and the radius of
markers is proportional to the size of the minimum gap $\Delta_{\rm min}$. The
color indicates the Chern number as shown in the legend. Panel (f) displays
tunneling conductance $\sigma(E)/\sigma_{N}$ calculated for $\theta\simeq
53.13^{\rm o}$, $\mu/t=-1.3$, $g_{0}=20$meV and temperatures ranging from 0 to
$1.2T_{c}$ as indicated by color scale. Curves for different temperatures have
been offset for clarity. Panels (g) and (h) illustrate the edge modes for
$C=4$ and $C=2$ topological phases, respectively, in $\theta_{2,5}$
configuration for parameters $\mu=-1.3t$ and $g_{0}=20,52$ meV. The energy
spectrum is shown for a bilayer system in the infinite strip geometry with
width of 90 unit cells. The color scale represents the normalized position
expectation value of the eigenstate along the direction of finite length.

In keeping with the strategy outlined above, we begin with a single CuO2 plane
per monolayer, as for Bi2201. The continuum theory for such a bilayer is
defined by the Hamiltonian

$\displaystyle{\cal H}$ $\displaystyle=$ $\displaystyle\sum_{{\bm{k}}\sigma
a}\xi_{{\bm{k}}a}c^{\dagger}_{{\bm{k}}\sigma a}c_{{\bm{k}}\sigma
a}+g\sum_{{\bm{k}}\sigma}\left(c^{\dagger}_{{\bm{k}}\sigma 1}c_{{\bm{k}}\sigma
2}+{\rm h.c.}\right)$ $\displaystyle+$
$\displaystyle\sum_{{\bm{k}}a}\left(\Delta_{{\bm{k}}a}c^{\dagger}_{{\bm{k}}\uparrow
a}c^{\dagger}_{-{\bm{k}}\downarrow a}+{\rm
h.c.}\right)-\sum_{{\bm{k}}a}\Delta_{{\bm{k}}a}\langle
c^{\dagger}_{{\bm{k}}\uparrow a}c^{\dagger}_{-{\bm{k}}\downarrow a}\rangle.$

Here, $c^{\dagger}_{{\bm{k}}\sigma a}$ creates an electron with crystal
momentum ${\bm{k}}$ and spin $\sigma$ in layer $a=1,2$ while $\xi_{{\bm{k}}a}$
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

$\Delta_{{\bm{k}}a}={\sum_{\bm{p}}}^{\prime}V_{{\bm{k}}{\bm{p}}}^{(a)}\langle
c_{-{\bm{p}}\downarrow a}c_{{\bm{p}}\uparrow a}\rangle.$ (8)

The prime on the summation indicates a restriction to momentum states with
energy within $\epsilon_{c}$ of the Fermi level. The Hamiltonian (II) should
be regarded as a mean-field approximation to the BCS pairing Hamiltonian with
an interaction term
$\sum_{{\bm{k}}{\bm{p}}}V_{{\bm{k}}{\bm{p}}}^{(a)}c^{\dagger}_{{\bm{k}}\uparrow
a}c^{\dagger}_{-{\bm{k}}\downarrow a}c_{-{\bm{p}}\downarrow
a}c_{{\bm{p}}\uparrow a}$, where $V_{{\bm{k}}{\bm{p}}}^{(a)}$ denotes the
interaction matrix element in layer $a$. We shall use a simple separable form

$V_{{\bm{k}}{\bm{p}}}^{(a)}=-2{\cal
V}\cos{(2\alpha_{\bm{k}})}\cos{(2\alpha_{\bm{p}})},$ (9)

where $\alpha_{\bm{k}}$ represents the polar angle of the vector ${\bm{k}}$.
This is known to yield a robust solution with $d_{x^{2}-y^{2}}$ symmetry for a
single CuO2 layer, namely
$\Delta_{{\bm{k}}a}=\Delta_{d}\cos{(2\alpha_{\bm{k}})}$.

In order to incorporate the twist, we take the interaction in layer 1 as in
Eq. (9), but we rotate the interaction in layer 2 by angle $\theta$ such that
$V_{{\bm{k}}{\bm{p}}}^{(2)}=-2{\cal
V}\cos{(2\alpha_{\bm{k}}-2\theta)}\cos{(2\alpha_{\bm{p}}-2\theta)}$. For the
sake of simplicity, we consider a circular Fermi surface generated by
$\xi_{{\bm{k}}a}=\hbar^{2}k^{2}/2m-\mu$ that remains invariant under rotation,
see Fig. 2(e). The problem posed by Hamiltonian (II) is then solved by
defining a four-component Nambu spinor $\Psi_{\bm{k}}=(c_{{\bm{k}}\uparrow
1},c^{\dagger}_{-{\bm{k}}\downarrow 1},c_{{\bm{k}}\uparrow
2},c^{\dagger}_{-{\bm{k}}\downarrow 2})^{T}$ in terms of which ${\cal
H}=\sum_{\bm{k}}\Psi^{\dagger}_{\bm{k}}h_{\bm{k}}\Psi_{\bm{k}}+E_{0}$. Here,
the Bogoliubov-de Gennes (BdG) Hamiltonian reads

$h_{\bm{k}}=\begin{pmatrix}\xi_{\bm{k}}&\Delta_{{\bm{k}}1}&g&0\\\
\Delta_{{\bm{k}}1}^{\ast}&-\xi_{\bm{k}}&0&-g\\\
g&0&\xi_{\bm{k}}&\Delta_{{\bm{k}}2}\\\
0&-g&\Delta_{{\bm{k}}2}^{\ast}&-\xi_{\bm{k}}\end{pmatrix}$ (10)

and
$E_{0}=\sum_{\bm{k}}2\xi_{\bm{k}}-\sum_{{\bm{k}}a}\Delta_{{\bm{k}}a}\langle
c^{\dagger}_{{\bm{k}}\uparrow a}c^{\dagger}_{-{\bm{k}}\downarrow a}\rangle$.
Diagonalizing $h_{\bm{k}}$ gives two pairs of energy eigenvalues $\pm
E_{{\bm{k}}\alpha}$ $(\alpha=1,2$) for each momentum ${\bm{k}}$. With the
Hamiltonian (II) expressed in diagonal form, the free energy of the system can
be calculated from the standard expression

${\cal F}_{\rm
BdG}=E_{0}-2k_{B}T\sum_{{\bm{k}}\alpha}\ln\left[2\cosh{(E_{{\bm{k}}\alpha}/2k_{B}T)}\right].$
(11)

By performing a systematic expansion of ${\cal F}_{\rm BdG}$ in powers of the
order parameter amplitudes $\Delta_{d}$, it is possible to ascertain various
coefficients entering the GL free energy (I). This confirms the form of the GL
coefficients $B$ and $C$ deduced previously on the basis of general symmetry
arguments, and that ${\cal K}$ in Eq. (4) is positive. This calculation is
summarized in Methods.

We proceed with a fully self-consistent solution which follows from minimizing
${\cal F}_{\rm BdG}$ with respect to $\Delta_{{\bm{k}}a}$ and can be performed
for any given ${\cal V}$, twist angle $\theta$ and temperature $T$. To
facilitate this calculation we assume that the order parameters in the two
layers have the same amplitude $\Delta_{d}$, but can differ in phase as in Eq.
(3). Free energy computed from Eq. (11) then shows the same qualitative
behaviour as the GL theory (4) with ${\cal K}>0$: for small twist angles,
${\cal F}_{\rm BdG}(\varphi)$ has a single minimum at $\varphi=0$, while for
twist angles close to $45^{\rm o}$ there are always two degenerate minima at
$\varphi=\pm\varphi_{\rm min}$, indicating formation of a ${\cal T}$-broken
topological phase. Fig. 2(f) shows the calculated $\varphi_{\rm min}$ as a
function of the twist angle $\theta$ for realistic Bi2212 parameters [33, 34,
35] $\Delta_{d}=40$ meV, $g=30$ meV, $\epsilon_{c}=60$ meV at $T=0$. We
observe an excellent agreement with the prediction of the GL theory when
${\cal K}=0.125$. The BdG theory also allows us to predict the minimum
excitation gap which, in this simple model, can be as large as 20 meV when
$\theta$ is close to $45^{\rm o}$.

The BdG theory can be formulated directly on the lattice by starting from the
Hubbard model with on-site repulsion and nearest-neighbour attraction. This is
known to produce a $d_{x^{2}-y^{2}}$ superconductor when applied to a single
CuO2 monolayer. The Hamiltonian is

$\displaystyle H$ $\displaystyle=$ $\displaystyle-\sum_{ij,\sigma
a}t_{ij}c^{\dagger}_{i\sigma a}c_{j\sigma a}-\mu\sum_{i\sigma a}n_{i\sigma a}$
(12) $\displaystyle+$
$\displaystyle\sum_{ij,a}V_{ij}n_{ia}n_{ja}-\sum_{ij\sigma}g_{ij}c^{\dagger}_{i\sigma
1}c_{j\sigma 2},$

where $t_{ij}$ encodes the normal-state band structure of the single layer,
$g_{ij}$ describes the interlayer tunneling and $V_{ij}$ denotes density-
density interactions. The mean-field calculations are performed in bilayer
geometries characterized by a twist vector ${\bm{v}}=(m,n)$ and commensurate
twist angle $\theta_{m,n}=2\arctan{(m/n)}$ as explained in Methods. Fig. 3
summarizes our main results. These lend further support to our conclusions
drawn on the basis of the GL and continuum BdG approaches and show additional
interesting features.

The lattice model confirms the onset of the ${\cal T}$-breaking phase below
critical temperature $\tilde{T}_{c}(\theta)$ whose dependence on $\theta$ is
consistent with the GL prediction, see Fig. 3(a,b). The results are in
reasonable quantitative agreement with the continuum BdG theory but the
lattice model gives a smaller spectral gap. This may be attributed to more
complicated Fermi surface geometry resulting from Brillouin zone folding that
accompanies the large real-space Moiré unit cell. The lattice model allows for
a direct evaluation of the Chern number $C$, as discussed in Methods. In
addition to the $C=4$ topological phase anticipated on the basis of the
continuum BdG theory, the $T=0$ phase diagrams in Fig. 3(c-e) reveal the
existence of $C=2$ and $C=0$ gapped phases that can be reached by varying the
chemical potential $\mu$ and the interlayer coupling strength $g_{0}$. In an
experiment, the former can be tuned over a wide range by oxygen annealing [1],
while the latter may depend on the twist angle $\theta$ and could be further
controlled by applying hydrostatic pressure as demonstrated in twisted
graphene [4]. As explained in Methods, the $C=0$ phase corresponds to a $d+is$
superconductor, while the $C=2$ phase can be effectively thought of as a
single-layer $d+id^{\prime}$ superconductor. For a strip geometry, the lattice
model predicts protected chiral edge modes traversing the bulk gap, Fig.
3(g,h), confirming the topological nature of the $d+id^{\prime}$ phase.

Additional interesting behaviour is observed at nonzero temperature. Fig. 10
in Methods shows a phase transition from $d+is$ to $d+id$ phase driven by
increasing $T$. As we discuss in more detail below, all these phenomena are
experimentally accessible using standard spectroscopic, thermodynamic and
transport techniques.

##  III Stability analysis and the Density functional theory

In order to establish the physical stability of the twisted double layer
structure and to estimate the interlayer coupling, we performed extensive
Density Functional Theory (DFT) simulations. We have employed the structure of
Bi2201; with a single CuO2 plane per monolayer, it has a direct connection to
the $N=1$ layer continuum and lattice models. Furthermore, the volume of the
unit cell is greatly reduced with respect to Bi2212, facilitating calculation
of the freestanding bilayer heterostructure. This analysis is focused on the
coupling between two CuO2 blocks; as for the previous section, we do not
expect the additional internal structure of each block in the $N>1$ case to
change our results substantively.

Figure 4: Crystal structure and results of the twisted double layer DFT
calculations. In panel (a), side view of the unit cell with unshifted
monolayers. Top view of the CuO2 planes are illustrated for the unshifted (b)
and shifted (c) unit cells. For the unit cell in (b), we plot the splitting at
the $X$ point as a function of the Cu - Cu distance in (d). A narrow window of
the twisted bilayer bandstructure near the $X$ point of the original Bi2201
Brillouin zone is provided in (e), at a Cu-Cu spacing of 11.97 Å . The pair of
bands used for defining the splitting in (d) are identified by the orange
arrow. A fit of the $X$-point splitting to $1/x$ is shown as a dashed red line
in (d). The equilibrium layer spacing is solved for as the minimum in the
system’s total energy, as illustrated by the inset. The total energy of the
bilayer unit cell is given for both unit cells in (b,c). The energy scale is
defined relative to the limit of two isolated monolayers.

Bi2201 supercells in the ${\bm{v}}=(1,2)$ geometry were constructed from a
bulk model defined with an in-plane Cu-Cu distance of $a=3.99$ Å, and unit-
cell height of $c=23.98$ Å. The tetragonal $I4/mmm$ space group was imposed to
minimize complications due to the small orthorhombicity of the actual
material. Two unit cells, each containing a single CuO2 plane and counter-
rotated by the twist angle were combined in a freestanding bilayer
construction, with a vacuum buffer of 15Å, see Fig. 4(a). We considered two
high-symmetry configurations of the bilayer unit cell, with the two layers
translated by (0,0) and ($\sqrt{5}a/2$,$\sqrt{5}a/2$), as illustrated in Fig.
4(b,c).The Kohn-Sham equations were solved self-consistently using the SCAN
meta-GGA exchange correlation potential [36], as implemented in VASP [37, 38].
As opposed to standard GGA, the recently developed SCAN functional
incorporates kinetic-energy corrections, supporting the study of mixed-bonding
environments [36]. In bilayer cuprates, weak van der Waals and ionic bonding
coexist, making this approach well-suited to the problem [39, 40]. This
circumvents the need for empirical van der Waals corrections to the GGA. We
repeated this process for various inter-bilayer distances to ascertain the
energetic stability of the bilayer compound. In the inset of panel (d), we
plot the total ground state energy, evaluated relative to the limit of two
independent monolayers. In both translational configurations, we find a
similar energy minimum near 12.87 Å. We note that this is the distance between
neighbouring CuO2 layers; the BiO layers which constitute the bilayer
interface are 9.45 Å closer together. These results are largely independent of
the exchange correlation potential used, as explored in Methods.

To estimate the strength of interlayer coupling between neighbouring CuO2
layers, we identify an avoided crossing at the $X$ point of the Bi2201
Brillouin zone, indicated in the bandstructure of Fig. 4(e). The size of this
splitting decays with the Cu-Cu distance, to a value of $\sim$ 5 meV at the
optimal spacing of 12.87 Å, as shown in Fig. 4(d). We emphasize that this
splitting does not provide a direct measure of the parameter $g_{0}$, but
establishes evidence for finite Cu-Cu interlayer coupling within the energy
scale surveyed in the text. It is important to note that there is precedent
for a renormalization of the DFT-derived interlayer coupling strength in
cuprates, as measured via photoemission. These renormalization factors vary
over a range of at least 2-5 in different materials [41]. Additional details
of the calculations presented here can be found in the Methods section.

##  IV Outlook

Our results establish a new avenue for engineering high-temperature
topological superconductors by combining monolayers of van der Waals-bonded
$d$-wave superconducting cuprates, such as Bi2Sr2CuO6+δ or Bi2Sr2CaCu2O8+δ
[1]. The topological phase is fully gapped, breaks time reversal symmetry and
supports protected chiral Majorana modes at the sample boundary. The chiral
$d\pm id^{\prime}$ phase arises robustly, driven by Cooper pair tunneling
between the layers when they are twisted by an angle close to 45o. Unlike in
graphene, where the twist angle must be tuned very accurately, we predict that
a wide range of angles (45${}^{\rm o}\pm 10^{\rm o}$) will produce nontrivial
physics in these systems.

The resulting fully gapped topological phase can be probed by conventional
spectroscopic or transport techniques. For example, on lowering the sample
temperature, angle resolved photoemission or tunneling spectroscopies [34, 35]
will detect the transition at $T=\tilde{T}_{c}(\theta)$ from the gapless
$d_{x^{2}-y^{2}}$ superconductor with Dirac nodes to a fully gapped excitation
spectrum. This change should also be observable via thermodynamic probes such
as those measuring the specific heat [42] and superfluid density [43, 44].
Signatures of the broken ${\cal T}$-symmetry characteristic of the $d\pm
id^{\prime}$ and $d+is$ phases may be identified by circular birefringence and
dichroism experiments [45]. In transport, the non-zero Chern number manifests
as quantized electronic contribution to the thermal Hall conductance
$\kappa_{xy}=C(\pi^{2}k_{B}^{2}T/3h)$ [22]. In mesoscopic samples,
longitudinal thermal conductance $\kappa_{xx}$ will also be quantized due to
the ballistic contribution of the edge states. The $d\pm id^{\prime}$ phase
exhibits anomalous interlayer Josephson effect due to the two-minima structure
of the the free energy depicted in Fig. 2(c); the Josephson current
$I_{J}=(2e/\hbar)\partial{\cal F}/\partial\varphi$ will develop pronounced
deviations from the usual sinusoidal dependence on the phase bias
$\Delta\varphi=\varphi-\varphi_{\rm min}$.

It is worth noting that existing experimental results on $c$-axis twist
Josephson tunneling in cuprates already offer some support for the ideas
advanced here. According to symmetry analysis, which assumes unbroken time
reversal and pure $d_{x^{2}-y^{2}}$ order parameter, the Josephson current
$I_{J}$ must vanish when $\theta=45^{\rm o}$ [23, 24, 25]. In Methods, we show
that for a $d+id^{\prime}$ state $I_{J}$, however, remains nonzero at
$\theta=45^{\rm o}$. Non-vanishing $I_{J}$ has often been interpreted as
evidence for conventional $s$-wave order parameter in experiments [46, 47].
Given what we know about the order parameter symmetry in cuprates today, one
may argue that a more plausible explanation for non-vanishing $I_{J}$ is the
spontaneously ${\cal T}$-broken phase advocated in this work.

The authors are indebted to D.A. Bonn, D.M. Broun, J.A. Folk, C. Kallin, C.
Li, É. Lantagne-Hurtubise, S. Plugge, S. Sahoo, O. Vafek and Z. Ye for
valuable discussions and correspondence. The work described here was supported
by NSERC. This research was undertaken thanks in part to funding from the Max
Planck-UBC-UTokyo Centre for Quantum Materials and the Canada First Research
Excellence Fund, Quantum Materials and Future Technologies Program. OC is
supported by International Doctoral Fellowship from UBC.

###  IV.1 Data availability

This manuscript contains no experimental data.

###  IV.2 Code availability

The complete code used to obtain results shown in Fig. 3 will be made publicly
available at https://github.com/ocanphys/tbcuprate/. DFT results shown in Fig.
4 were obtained VASP [37, 38] . Input files are also available through our
repository.

## References

  * Yu  _et al._ [2019] Y. Yu, L. Ma, P. Cai, R. Zhong, C. Ye, J. Shen, G. D. Gu, X. H. Chen, and Y. Zhang, High-temperature superconductivity in monolayer bi2sr2cacu2o8+$\delta$, Nature 575, 156 (2019).
  * Cao  _et al._ [2018a] 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, Correlated insulator behaviour at half-filling in magic-angle graphene superlattices, Nature (London) 556, 80 (2018a).
  * Cao  _et al._ [2018b] Y. Cao, V. Fatemi, S. Fang, K. Watanabe, T. Taniguchi, E. Kaxiras, and P. Jarillo-Herrero, Unconventional superconductivity in magic-angle graphene superlattices, Nature (London) 556, 43 (2018b).
  * Yankowitz  _et al._ [2019] M. Yankowitz, S. Chen, H. Polshyn, Y. Zhang, K. Watanabe, T. Taniguchi, D. Graf, A. F. Young, and C. R. Dean, Tuning superconductivity in twisted bilayer graphene, Science 363, 1059 (2019).
  * Sharpe  _et al._ [2019] A. L. Sharpe, E. J. Fox, A. W. Barnard, J. Finney, K. Watanabe, T. Taniguchi, M. A. Kastner, and D. Goldhaber-Gordon, Emergent ferromagnetism near three-quarters filling in twisted bilayer graphene, Science 365, 605 (2019).
  * Lu  _et al._ [2019] 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, Superconductors, orbital magnets and correlated states in magic-angle bilayer graphene, Nature 574, 653 (2019).
  * Wang  _et al._ [2019] 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. Pasupathy, and C. R. Dean, Magic continuum in twisted bilayer WSe2, arXiv e-prints , arXiv:1910.12147 (2019), arXiv:1910.12147 [cond-mat.mes-hall] .
  * Bistritzer and MacDonald [2011] R. Bistritzer and A. H. MacDonald, Moiré bands in twisted double-layer graphene, Proceedings of the National Academy of Science 108, 12233 (2011).
  * Nam and Koshino [2017] N. N. T. Nam and M. Koshino, Lattice relaxation and energy band modulation in twisted bilayer graphene, Phys. Rev. B 96, 075311 (2017).
  * Kang and Vafek [2018] J. Kang and O. Vafek, Symmetry, maximally localized wannier states, and a low-energy model for twisted bilayer graphene narrow bands, Phys. Rev. X 8, 031088 (2018).
  * Zou  _et al._ [2018] L. Zou, H. C. Po, A. Vishwanath, and T. Senthil, Band structure of twisted bilayer graphene: Emergent symmetries, commensurate approximants, and Wannier obstructions, Phys. Rev. B 98, 085435 (2018).
  * Po  _et al._ [2018] H. C. Po, L. Zou, A. Vishwanath, and T. Senthil, Origin of mott insulating behavior and superconductivity in twisted bilayer graphene, Phys. Rev. X 8, 031089 (2018).
  * Xie and MacDonald [2020] M. Xie and A. H. MacDonald, Nature of the correlated insulator states in twisted bilayer graphene, Phys. Rev. Lett. 124, 097601 (2020).
  * Koshino  _et al._ [2018] M. Koshino, N. F. Q. Yuan, T. Koretsune, M. Ochi, K. Kuroki, and L. Fu, Maximally localized wannier orbitals and the extended hubbard model for twisted bilayer graphene, Phys. Rev. X 8, 031087 (2018).
  * Guinea and Walet [2018] F. Guinea and N. R. Walet, Electrostatic effects, band distortions, and superconductivity in twisted graphene bilayers, Proceedings of the National Academy of Science 115, 13174 (2018).
  * Kang and Vafek [2019] J. Kang and O. Vafek, Strong Coupling Phases of Partially Filled Twisted Bilayer Graphene Narrow Bands, Phys. Rev. Lett. 122, 246401 (2019).
  * Hejazi  _et al._ [2019] K. Hejazi, C. Liu, H. Shapourian, X. Chen, and L. Balents, Multiple topological transitions in twisted bilayer graphene near the first magic angle, Phys. Rev. B 99, 035111 (2019).
  * Andrei and MacDonald [2020] E. Y. Andrei and A. H. MacDonald, Graphene Bilayers with a Twist, arXiv e-prints , arXiv:2008.08129 (2020), arXiv:2008.08129 [cond-mat.mes-hall] .
  * Elliott and Franz [2015] S. R. Elliott and M. Franz, Colloquium: Majorana fermions in nuclear, particle, and solid-state physics, Rev. Mod. Phys. 87, 137 (2015).
  * Laughlin [1998] R. B. Laughlin, Magnetic induction of ${\mathit{d}}_{{\mathit{x}}^{2}-{\mathit{y}}^{2}}+{\mathrm{id}}_{\mathrm{xy}}$ order in high- ${T}_{c}$ superconductors, Phys. Rev. Lett. 80, 5188 (1998).
  * Franz and Tešanović [1998] M. Franz and Z. Tešanović, Self-consistent electronic structure of a ${\mathit{d}}_{{\mathit{x}}^{2}-{\mathit{y}}^{2}}$ and a ${\mathit{d}}_{{\mathit{x}}^{2}-{\mathit{y}}^{2}}+{\mathit{id}}_{\mathit{xy}}$ vortex, Phys. Rev. Lett. 80, 4763 (1998).
  * Vishwanath [2001] A. Vishwanath, Quantized thermal hall effect in the mixed state of $\mathit{d}$-wave superconductors, Phys. Rev. Lett. 87, 217004 (2001).
  * Kuboki and Sigrist [1996] K. Kuboki and M. Sigrist, Proximity-induced time-reversal symmetry breaking at josephson junctions between unconventional superconductors, Journal of the Physical Society of Japan 65, 361 (1996).
  * Bille  _et al._ [2001] A. Bille, R. A. Klemm, and K. Scharnberg, Models of c-axis twist josephson tunneling, Phys. Rev. B 64, 174507 (2001).
  * Yokoyama  _et al._ [2007a] T. Yokoyama, S. Kawabata, T. Kato, and Y. Tanaka, Theory of macroscopic quantum tunneling in high-${T}_{c}$ $c$-axis josephson junctions, Phys. Rev. B 76, 134501 (2007a).
  * Yang  _et al._ [2018] Z. Yang, S. Qin, Q. Zhang, C. Fang, and J. Hu, $\pi$/2-josephson junction as a topological superconductor, Phys. Rev. B 98, 104515 (2018).
  * Pathak  _et al._ [2010] S. Pathak, V. B. Shenoy, and G. Baskaran, Possible high-temperature superconducting state with a $d+id$ pairing symmetry in doped graphene, Phys. Rev. B 81, 085431 (2010).
  * Nandkishore  _et al._ [2012] R. Nandkishore, L. S. Levitov, and A. V. Chubukov, Chiral superconductivity from repulsive interactions in doped graphene, Nature Physics 8, 158 (2012).
  * Rice and Sigrist [1995] T. M. Rice and M. Sigrist, Sr2ruo4: an electronic analogue of 3he?, Journal of Physics: Condensed Matter 7, L643 (1995).
  * Ishida  _et al._ [1998] K. Ishida, H. Mukuda, Y. Kitaoka, K. Asayama, Z. Q. Mao, Y. Mori, and Y. Maeno, Spin-triplet superconductivity in sr2ruo4 identified by 17o knight shift, Nature 396, 658 (1998).
  * Kallin and Berlinsky [2009] C. Kallin and A. J. Berlinsky, Is sr2ruo4a chiral p-wave superconductor?, Journal of Physics: Condensed Matter 21, 164210 (2009).
  * Pustogow  _et al._ [2019] 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, Constraints on the superconducting order parameter in sr2ruo4 from oxygen-17 nuclear magnetic resonance, Nature 574, 72 (2019).
  * Andersen  _et al._ [1995] O. Andersen, A. Liechtenstein, O. Jepsen, and F. Paulsen, Lda energy bands, low-energy hamiltonians, $t^{\prime}$, $t^{\prime\prime}$, $t_{\perp}(k)$, and $j_{\perp}$, arXiv preprint cond-mat/9509044 (1995).
  * Damascelli  _et al._ [2003] A. Damascelli, Z. Hussain, and Z.-X. Shen, Angle-resolved photoemission studies of the cuprate superconductors, Rev. Mod. Phys. 75, 473 (2003).
  * Fischer  _et al._ [2007] O. Fischer, M. Kugler, I. Maggio-Aprile, C. Berthod, and C. Renner, Scanning tunneling spectroscopy of high-temperature superconductors, Rev. Mod. Phys. 79, 353 (2007).
  * Sun  _et al._ [2015] J. Sun, A. Ruzsinszky, and J. P. Perdew, Strongly constrained and appropriately normed semilocal density functional, Phys. Rev. Lett. 115, 036402 (2015).
  * Kresse and Hafner [1993] G. Kresse and J. Hafner, Ab initio molecular dynamics for open-shell transition metals, Phys. Rev. B 48, 13115 (1993).
  * Kresse and Furthmüller [1996] G. Kresse and J. Furthmüller, Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set, Phys. Rev. B 54, 11169 (1996).
  * Furness  _et al._ [2018] J. W. Furness, Y. Zhang, C. Lane, I. G. Buda, B. Barbiellini, R. S. Markiewicz, A. Bansil, and J. Sun, An accurate first-principles treatment of doping-dependent electronic structure of high-temperature cuprate superconductors, Communications Physics 1, 11 (2018).
  * Zhang  _et al._ [2020] Y. Zhang, C. Lane, J. W. Furness, B. Barbiellini, J. P. Perdew, R. S. Markiewicz, A. Bansil, and J. Sun, Competing stripe and magnetic phases in the cuprates from first principles, Proceedings of the National Academy of Sciences 117, 68 (2020).
  * Markiewicz  _et al._ [2005] R. S. Markiewicz, S. Sahrakorpi, M. Lindroos, H. Lin, and A. Bansil, One-band tight-binding model parametrization of the high-${T}_{c}$ cuprates including the effect of ${k}_{z}$ dispersion, Phys. Rev. B 72, 054519 (2005).
  * Riggs  _et al._ [2011] 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, Heat capacity through the magnetic-field-induced resistive transition in an underdoped high-temperature superconductor, Nature Physics 7, 332 (2011).
  * Tallon  _et al._ [2003] J. L. Tallon, J. W. Loram, J. R. Cooper, C. Panagopoulos, and C. Bernhard, Superfluid density in cuprate high-${T}_{c}$ superconductors: A new paradigm, Phys. Rev. B 68, 180501 (2003).
  * Hosseini  _et al._ [2004] A. Hosseini, D. M. Broun, D. E. Sheehy, T. P. Davis, M. Franz, W. N. Hardy, R. Liang, and D. A. Bonn, Survival of the $d$-wave superconducting state near the edge of antiferromagnetism in the cuprate phase diagram, Phys. Rev. Lett. 93, 107003 (2004).
  * Yip and Sauls [1992] S. K. Yip and J. A. Sauls, Circular dichroism and birefringence in unconventional superconductors, Journal of Low Temperature Physics 86, 257 (1992).
  * Li  _et al._ [1997] Q. Li, Y. Tsay, M. Suenaga, G. Gu, and N. Koshizuka, Superconducting coupling along the c-axis of [001] twist grain-boundaries in bi2sr2cacu2o${}_{8}+\delta$ bicrystals, Physica C: Superconductivity 282-287, 1495 (1997).
  * Zhu  _et al._ [2019] Y. Zhu, M. Liao, Q. Zhang, F. Meng, R. Zhong, J. Schneeloch, G. Gu, L. Gu, X. Ma, D. Zhang, _et al._ , Isotropic josephson tunneling in c-axis twist bicrystals of bi2sr2cacu2o8+$\delta$, arXiv preprint arXiv:1903.07965 (2019).
  * Harrison [2012] W. Harrison,  _Electronic Structure and the Properties of Solids: The Physics of the Chemical Bond_ , Dover Books on Physics (Dover Publications, 2012).
  * Bansil  _et al._ [2005] A. Bansil, M. Lindroos, S. Sahrakorpi, and R. Markiewicz, Influence of the third dimension of quasi-two-dimensional cuprate superconductors on angle-resolved photoemission spectra, Physical Review B 71, 012503 (2005).
  * Slezak  _et al._ [2008] J. Slezak, J. Lee, M. Wang, K. McElroy, K. Fujita, B. Andersen, P. Hirschfeld, H. Eisaki, S. Uchida, and J. Davis, Imaging the impact on cuprate superconductivity of varying the interatomic distances within individual crystal unit cells, Proceedings of the National Academy of Sciences 105, 3203 (2008).
  * Coleman [2015] P. Coleman,  _Introduction to many-body physics_ (Cambridge University Press, 2015).
  * Bernevig and Hughes [2013] B. A. Bernevig and T. L. Hughes,  _Topological insulators and topological superconductors_ (Princeton university press, 2013).
  * Feng  _et al._ [2001] D. L. Feng, N. P. Armitage, D. H. Lu, A. Damascelli, J. P. Hu, P. Bogdanov, A. Lanzara, F. Ronning, K. M. Shen, H. Eisaki, C. Kim, Z.-X. Shen, J.-i. Shimoyama, and K. Kishio, Bilayer splitting in the electronic structure of heavily overdoped ${\mathrm{bi}}_{2}{\mathrm{sr}}_{2}{\mathrm{cacu}}_{2}{O}_{8+\delta}$, Phys. Rev. Lett. 86, 5550 (2001).
  * Yokoyama  _et al._ [2007b] T. Yokoyama, S. Kawabata, T. Kato, and Y. Tanaka, Theory of macroscopic quantum tunneling in high-t c c-axis josephson junctions, Physical Review B 76, 134501 (2007b).
  * Perdew  _et al._ [1996] J. P. Perdew, K. Burke, and M. Ernzerhof, Generalized gradient approximation made simple, Phys. Rev. Lett. 77, 3865 (1996).
  * Grimme [2006] S. Grimme, Semiempirical gga-type density functional constructed with a long-range dispersion correction, Journal of Computational Chemistry 27, 1787 (2006), https://onlinelibrary.wiley.com/doi/pdf/10.1002/jcc.20495 .

##  V Methods

###  V.1 Microscopic derivation of GL coefficients

From the continuum BdG model, defined by the matrix Hamiltonian (10), it is
possible to derive the corresponding GL theory by directly expanding the free
energy ${\cal F}_{\rm BdG}$ given in Eq. (11) in the powers of the
superconducting order parameters. This procedure is made feasible by the fact
that explicit expressions for the two positive-energy eigenvalues of
$h_{\bm{k}}$ can be obtained and are given by

$E_{{\bm{k}}\alpha}=\sqrt{(\Delta_{{\bm{k}}1}^{2}+\Delta_{{\bm{k}}2}^{2})/2+\xi_{\bm{k}}^{2}+g^{2}+(-)^{\alpha}D_{\bm{k}}}$
(13)

where
$D_{k}^{2}=(\Delta_{{\bm{k}}1}^{2}-\Delta_{{\bm{k}}2}^{2})^{2}/4+g^{2}(\Delta_{{\bm{k}}1}^{2}+\Delta_{{\bm{k}}2}^{2}+4\xi_{\bm{k}}^{2})-2g^{2}\Delta_{{\bm{k}}1}\Delta_{{\bm{k}}2}\cos\varphi$
and $\alpha=1,2$. We now express the two $d$-wave gap functions as

$\Delta_{{\bm{k}}1}=\psi\cos(2\alpha_{\bm{k}}),\ \ \
\Delta_{{\bm{k}}2}=\psi\cos(2\alpha_{\bm{k}}-2\theta),$ (14)

substitute into (13) and expand the free energy (11) in the powers of the
order parameter amplitude $\psi$ to fourth order. This expansions results in
many terms and is best carried out with the help of Mathematica or similar
software capable of symbolic manipulations.

We are primarily interested in terms governing the phase difference between
the two layers, that is, coefficients $B$ and $C$ in the GL free energy (I).
These can be isolated by focusing on terms containing powers of $\cos\varphi$.
Coefficient $B$ is thus found as a prefactor of $\psi^{2}\cos{\varphi}$ in the
above expansion and reads

$B=\cos{(2\theta)}N_{F}\sum_{s=\pm}\int_{0}^{\epsilon_{c}}d\xi\frac{-\pi
sg}{\xi|\xi-sg|}\frac{\tanh{|\xi-sg|}}{2k_{B}T}$ (15)

In deriving this expression we replaced the crystal momentum sum using the
usual prescription $\sum_{\bm{k}}\to
N_{F}\int_{-\epsilon_{c}}^{\epsilon_{c}}d\xi\int_{0}^{2\pi}d\phi$ where
$N_{F}=L^{2}m/\hbar^{2}$ is the normal-state density of states at the Fermi
level, $L$ is the system size and the angular integral has been performed. The
remaining integral over the energy variable cannot be completed in the closed
form but one can show that the integrand $f_{B}(\xi)$ is negative for all
values of $\xi$ and thus the integral itself is also negative, as expected on
physical grounds. This is illustrated in Fig. 5(a).

Figure 5: Analysis of the energy integrands. (a) Integrand $f_{B}(\xi)$
entering the expression Eq. (15) for the GL theory coefficient $B$. (b)
Integrand $f_{C}(\xi)$ entering the expression for the GL coefficient $C$
compared to $f_{B}(\xi)^{2}$. Note that the amplitude of the latter is scaled
by a constant factor of 1/30 for easier comparison. The temperature is chosen
such that $k_{B}T=0.2g$ but the qualitative properties of the integrands
remain the same in a wide range of temperatures.

Coefficient $C$ can likewise be identified as a prefactor of the
$\psi^{4}\cos{2\varphi}$ term in the expansion and has the form
$C=(2+\cos{4\theta})I_{C}$. Here
$I_{C}=N_{F}\int_{0}^{\epsilon_{c}}f_{C}(\xi)d\xi$ is the corresponding energy
integral similar to Eq. (15). Its integrand $f_{C}(\xi)$ has a more
complicated structure and for the sake of brevity we do not show it
explicitly. Importantly $f_{C}(\xi)$ is positive for all $\xi$ and, as
anticipated in the main text, is approximately equal to the integrand of Eq.
(15) squared, as shown in Fig. 5(b). We thus conclude that the continuum BdG
model indeed supports the GL phenomenology deduced in the main text on the
basis of general physical arguments. Specifically the calculation explicitly
shows that prefactor ${\cal K}$ of the $\cos{2\varphi}$ term in Eq. (4) is
positive which we argued was the key condition underlying the emergence of the
${\cal T}$-broken topological phase in the twisted $d$SC bilayers.

###  V.2 Self-consistent mean field theory on the lattice

To describe the two-layer system microscopically, we employ the Hubbard
Hamiltonian (12) and treat the interaction term at the mean field level to
account for $d$-wave superconductivity in each layer. The MF Hamiltonian
becomes

$\displaystyle H$ $\displaystyle=$ $\displaystyle-t\sum_{\langle
ij\rangle\sigma a}c^{\dagger}_{i\sigma a}c_{j\sigma
a}-t^{\prime}\sum_{\langle\langle ij\rangle\rangle\sigma
a}c^{\dagger}_{i\sigma a}c_{j\sigma a}-\mu\sum_{i\sigma a}n_{i\sigma a}$ (16)
$\displaystyle+$ $\displaystyle\sum_{\langle ij\rangle
a}\left(\Delta_{ij,a}c^{\dagger}_{i\uparrow a}c^{\dagger}_{j\downarrow a}+{\rm
h.c.}\right)-\sum_{ij\sigma}g_{ij}c^{\dagger}_{i\sigma 1}c_{j\sigma 2},$

where $t$ and $t^{\prime}$ are the hopping amplitudes between first and second
neighbour sites on the square lattice, $\mu$ is the chemical potential that
controls on-site particle density $n_{i\sigma a}$ and $\Delta_{ij,a}=V\langle
c_{i\uparrow a}c_{j\downarrow a}\rangle$ is the complex order parameter. It is
believed that the interlayer tunneling processes in cuprates are mediated, to
leading order, by the spherically symmetric Cu $4s$ orbitals [33], 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 phenomenological form for the
interlayer tunneling amplitude that decays exponentially with distance:

$g_{ij}=g_{0}e^{-(r_{ij}-c)/\rho},$ (17)

where $r_{ij}^{2}=c^{2}+d_{ij}^{2}$, with $d_{ij}$ denoting the in-plane
separation between two sites $i$ and $j$ on different layers and $c$ being the
interlayer distance. In practice, we limit these terms to a few unit cells,
beyond which the amplitudes are negligible. The characteristic radial extent
of the Cu $4s$ orbital has been denoted by $\rho$. We have also considered
interlayer couplings of the form $g_{ij}=g_{0}c^{2}/r_{ij}^{2}$ [48]. Self
consistent calculations still result in topological phases, but this form of
coupling, when compared to (17) for same $g_{0}$, effectively corresponds to
stronger interlayer coupling as the decay of $1/r^{2}$ is slower compared to
the exponential.

While the above real space Hamiltonian is valid for any twist angle, the
bilayer forms a crystal only for a discrete set of commensurate rotation
angles. To carry out the lattice mean-field calculations we focus on
commensurate twist geometries that are amenable to numerical study using the
standard Bloch representation. We note that any commensurate twist can be
described by an integer-valued ‘twist’ vector ${\bm{v}}=(m,n)$ as shown in
Fig. 6. The corresponding moire unit cell comprises of $2q=2(m^{2}+n^{2})$
sites and the twist angle is determined by the relation

$\theta_{m,n}=2\arctan{(m/n)}.$ (18)

This can be understood as a rotation of two perfectly aligned square lattices
in opposite directions by $\theta_{m,n}/2$, which leads to the sites
originally at the locations ${\bm{v}}$ and $\bar{{\bm{v}}}=(-m,n)$ to lie atop
each other. The unit cells for a selection of commensurate twist angles close
to $45^{\rm o}$ are illustrated in Fig. 6.

Figure 6: Geometry of the commensurate unit cells for bilayers with relative
twists of (a) $\theta_{1,2}\simeq 53.13^{\rm o}$(b) $\theta_{2,5}\simeq
43.60^{\rm o}$ and (c) $\theta_{5,12}\simeq 45.24^{\rm o}$. Twisting procedure
is illustrated in panel (a) where orange and blue arrows indicate integer
valued vectors ${\bm{v}}=(m,n)$ and $\bar{{\bm{v}}}=(-m,n)$ introduced in the
text. Moire unit cell is shown by the shaded area.

To model the Fermi surface of Bi2212 we follow Ref. [24] and choose $t=153$
meV, $t^{\prime}=-0.45t$ with $\mu=-1.35t$ near optimal doping. Other
parameters that determine the interlayer tunneling amplitudes are as follows:
the average spatial extent of electrons in the $4s$ orbitals is $\rho\approx
4\text{a.u.}=2.11$\mathrm{\SIUnitSymbolAngstrom}$$, the typical inter-layer
distance is $c\approx 12$\mathrm{\SIUnitSymbolAngstrom}$$ [49] and the in-
plane separation between the Cu atoms is $a\approx
5.4$\mathrm{\SIUnitSymbolAngstrom}$$ [50].

In the remainder of this section we outline the procedure for obtaining gap
equations for the order parameters $\Delta_{ij,a}$ using the imaginary time
path integral formalism. After decoupling the attractive pairing potential in
the Cooper channel with the Hubbard-Stratonovich transformation, we obtain the
action corresponding to the Hamiltonian (16) (cf. Ref. [51] for details)

$S=\sum_{{\bm{k}}n}\Psi^{\dagger}_{{\bm{k}}n}[-i\omega_{n}+h_{\bm{k}}]\Psi_{{\bm{k}}n}+\frac{\beta\cal
N}{V}\sum_{\begin{subarray}{c}i,j\in\text{u.c.}\\\
a\end{subarray}}\Delta_{ij,a}\Delta_{ij,a}^{\ast}.$ (19)

Here $h_{\bm{k}}$ is the BdG Hamiltonian that follows from Fourier
transforming Eq. (16), $\beta=1/k_{B}T$ is the inverse temperature,
$\omega_{n}=(2n+1)\pi/\beta$ are the Matsubara frequencies, $V$ is the
strength of attractive nearest neighbour interactions, ${\cal N}$ is the
number of unit cells and
$\Psi_{\bm{k}}=(\Phi_{{\bm{k}}1},\Phi_{{\bm{k}}2})^{T}$ is the Nambu spinor
such that $\Phi_{{\bm{k}}a}=(c^{(1)}_{{\bm{k}}\uparrow
a},\ldots,c^{(q)}_{{\bm{k}}\uparrow a},c^{(1){\dagger}}_{-{\bm{k}}\downarrow
a},\ldots,c^{(q){\dagger}}_{-{\bm{k}}\downarrow a})^{T}$, with the
superscripts $1,\ldots,q$ labelling the sublattice degrees of freedom in layer
$a$.

The effective action for $\Delta_{ij,a}$ can be determined by integrating out
fermionic degrees of freedom

$S_{\text{eff}}=-\sum_{{\bm{k}}n}\text{Tr}\log[-\mathcal{G}({\bm{k}},i\omega_{n})^{-1}]+\frac{\beta{\cal
N}}{V}\sum_{\begin{subarray}{c}i,j\in\text{u.c.}\\\
a\end{subarray}}\Delta_{ij,a}\Delta_{ij,a}^{\ast},$ (20)

where we defined the Matsubara Green’s function
$\mathcal{G}({\bm{k}},i\omega_{n})=-(-i\omega_{n}+h_{\bm{k}})^{-1}$. To
evaluate the saddle point condition for this effective action, $\partial
S_{\text{eff}}/\partial\Delta_{ij,a}^{\ast}=0$, we employ the identity
$\partial_{\Delta}(\text{Tr}\log A)=\text{Tr}[\partial_{\Delta}AA^{-1}]$.
After performing the Matsubara sum, this simplifies to

$\Delta_{ij,a}=-\frac{V}{{\cal N}}\sum_{{\bm{k}}}\text{Tr}\left[\frac{\partial
h_{\bm{k}}}{\partial\Delta_{ij,a}^{\ast}}U_{\bm{k}}n_{F}(E_{\bm{k}})U_{\bm{k}}^{\dagger}\right].$
(21)

In the above, $U_{\bm{k}}$ is the unitary matrix that diagonalizes the Bloch
Hamiltonian as $U_{\bm{k}}^{\dagger}h_{\bm{k}}U_{\bm{k}}=E_{\bm{k}}$ and
$n_{F}(E_{\bm{k}})$ is a diagonal matrix where Fermi function is applied to
the eigenvalues $E_{\bm{k}}$. For decoupled layers ($g_{0}=0$), with the tight
binding parameters given above, a self-consistent solution converges to a
$d$-wave form with maximal gap $\Delta_{\rm max}=40$ meV, when $V=146$ meV. We
use this value for calculations where two layers are coupled ($g_{0}\neq 0$).
Fig. 3 in the main text summarizes the results obtained from these
calculations.

###  V.3 Topological phases

Here, we discuss in greater detail the various phases found in the self-
consistent treatment of the lattice model.

To reliably determine the topology of the band structure, we turn to the Chern
number, which is defined as an integral of the Berry curvature over the full
Brillouin zone: $C=\frac{1}{2\pi}\int_{\rm BZ}\Omega({\bm{k}})~{}d^{2}k$.
Computation of the Berry curvature directly from the eigenstates is
complicated by the arbitrary phases that accompany numerical diagonalization.
To counter this, we rely on a gauge independent formulation where [52]

$\Omega({\bm{k}})=\sum_{m\neq n}\frac{\langle
n|\nabla_{\bm{k}}h_{\bm{k}}|m\rangle\times\langle
m|\nabla_{\bm{k}}h_{\bm{k}}|n\rangle}{(E_{m}-E_{n})^{2}}.$ (22)

Here $|n\rangle$ and $|m\rangle$ are eigenstates of the Bloch Hamiltonian with
energies $E_{n}$ and $E_{m}$ respectively.

The $C=\pm 4$ and $\pm 2$ phases (indicated in Fig. 3(e) in red and green
respectively) can be understood by considering the continuum model with a
circular Fermi surface in the normal state, which we defined in Eq. (II).
Similar arguments hold for the more realistic Fermi surface with hole pockets.
When two layers are decoupled ($g=0$), Fermi surfaces of each layer overlap,
Fig. 7(a), and each layer hosts four Dirac cones in the presence of a
$d_{x^{2}-y^{2}}$ order parameter. Upon breaking of ${\cal T}$, each Dirac
point is gapped out and can contribute a Berry flux of $\pm\pi$, thereby
allowing for the maximal value $C=\pm 4$. As we turn on the interlayer
coupling, bands from each monolayer split, Fig. 7(b). If the interlayer
coupling $g$ is strong enough to push one of the bands above the Fermi level,
half of the Dirac cones disappear, leaving four Dirac cones behind. In this
case, the maximum Chern number is $\pm 2$, a situation realized in the phase
diagram Fig. 3(e) for large $g_{0}$.

Figure 7: Schematic picture for the normal state band structure of the
continuum model defined in Eq. (II). Both electron and hole bands that enter
the BdG Hamiltonian are shown while the Fermi level is indicated by a dashed
red line. Panel (a) shows decoupled layers with doubly degenerate bands
yielding a circular Fermi surface. In the presence of a $d$-wave order
parameter we have eight Dirac cones. (b) Interlayer coupling $g$ splits the
bands and Fermi surface now consists of a pair of concentric circles. (c) The
inner Fermi circle shrinks to a point and half of the Dirac cones disappear.
(d) Four Dirac cones survive in the strong interlayer coupling ($g>\mu$) case.

It is interesting to note that the spectral gap tends to be larger in the
$C=2$ phase. We attribute this to the fact that the $C=2$ phase is stabilized
at stronger values of the interlayer tunneling $g_{0}$ which enhances the
Cooper pair tunneling amplitude. Moreover, in second order perturbation
theory, the gap is expected to scale as $g_{0}^{2}/\Delta_{0}$. Temperature
evolution of the corresponding tunneling density of states is shown in Fig. 8.

Figure 8: Evolution of the density of states with temperature in the $C=2$
phase with parameters $g_{0}=40$ meV and $\mu=-1.35t$. Temperature $T/T_{c}$
shown as color scale. The superconducting gap persists up until $T=T_{c}$.

The gapped, topologically trivial $C=0$ phase in the phase diagram can be
understood as follows. As a representative of this phase, we consider the
point $g_{0}=20$ meV and $\mu/t=-1.3$ for ${\bm{v}}=(1,3)$ in the phase
diagram Fig. 3(c). A self-consistent solution of the gap equation (21) yields
the following forms for order parameters on the two layers

$\displaystyle\Delta_{{\bm{k}}1}=$
$\displaystyle\Delta_{0}(\cos{k_{x}}-e^{i\phi}\cos{k_{y}})$ (23)
$\displaystyle\Delta_{{\bm{k}}2}=$
$\displaystyle\Delta_{0}(\cos{k_{y}}-e^{i\phi}\cos{k_{x}})$ (24)

where $\Delta_{1}$ and $\Delta_{2}$ correspond to top and bottom layers
respectively. The numerically obtained phase $\phi$ tends to be close to
$-\pi/4$. If we define $\Delta^{d}=\Delta_{0}e^{i\phi/2}\cos{\phi/2}$ and
$\Delta^{s}=\Delta_{0}e^{i\phi/2}\sin{\phi/2}$, order parameters can be cast
as

$\displaystyle\Delta_{{\bm{k}}1}=$
$\displaystyle\Delta^{d}(\cos{k_{x}}-\cos{k_{y}})-i\Delta^{s}(\cos{k_{x}}+\cos{k_{y}})$
$\displaystyle\Delta_{{\bm{k}}2}=$
$\displaystyle-\Delta^{d}(\cos{k_{x}}-\cos{k_{y}})-i\Delta^{s}(\cos{k_{x}}+\cos{k_{y}}),$

which shows that both layers have $d+is$ type superconducting order parameter.
The transition between the $d+is$ and the $d+id$ phases can be seen by tuning
the chemical potential $\mu$ for a fixed interlayer coupling $g_{0}$, as shown
in Fig. 9.

Figure 9: The topological phase transition between $C=0$ ($\mu<-1.34t$) and
$C=4$ ($\mu>-1.34t$) phases is marked by a closing of the gap as the chemical
potential is tuned while keeping the interlayer coupling fixed at $g_{0}=20$
meV and $\theta_{2,5}\simeq 43.6^{\circ}$

###  V.4 Temperature dependence

The self-consistent gap equations can be also solved at non-zero temperature
$T$. Although the Chern number is well defined only at $T=0$ we may still
compute $C$ from Eq. (22) using the Bloch eigenstates obtained at $T>0$. The
resulting integer then describes the underlying mean-field band structure of
the system at non-zero temperature. In this case, there will be no strictly
quantized physical response (such as the thermal Hall conductance), but a
nonzero $C$ indicates the existence of protected chiral edge states that can
be observed experimentally. In this sense the $T>0$ phases with non-zero $C$
can be regarded as topological and their properties related are to the $T=0$
counterparts.

In Fig. 10 we compare phase diagrams at $T=T_{c}/2$ and $T=0$. We observe
that, remarkably, increasing temperature tends to suppress the topologically
trivial $d+is$ phase in favor of topological $d+id^{\prime}$ phases.
Specifically at some points in the phase diagram increasing temperature can
trigger a transition from the $d+is$ phase to the $d+id^{\prime}$. This is
illustrated in Fig. 10(c). As a result, we find that at elevated temperatures
(but still well below $T_{c}$) essentially all of the phase diagram is taken
over by the $C=4$ and $C=2$ topological phases.

Figure 10: Temperature dependence of the phase diagram is shown for
$\theta_{2,5}\simeq 43.6^{\circ}$ at $T=0$ and $T=T_{c}/2$ in panels (a) and
(b) respectively. Panel (c) shows the minimum gap as a function of temperature
for $g_{0}=32$meV and $\mu=-1.35t$ (the parameters are indicated in panels (a)
and (b) by a dashed circle)

###  V.5 Multiple ${\rm CuO}_{2}$ planes

We now briefly discuss the extension of our theory to systems with $N>1$ CuO2
planes per monolayer, focusing on compounds with $N=2$ such as Bi2212. The GL
theory can be generalized by including complex order parameters $\chi_{a}$
describing the extra layers as illustrated in Fig. 11(a). For $N=2$, the
corresponding free energy density reads

$\displaystyle{\cal F}$ $\displaystyle=$ $\displaystyle{\cal
F}_{0}+B(\psi_{1}\psi_{2}^{\ast}+{\rm c.c.})+C(\psi_{1}^{2}\psi_{2}^{\ast
2}+{\rm c.c.})$ (25) $\displaystyle+$
$\displaystyle\sum_{a=1}^{2}[B^{\prime}(\psi_{a}\chi_{a}^{\ast}+{\rm
c.c.})+C^{\prime}(\psi_{a}^{2}\chi_{a}^{\ast 2}+{\rm c.c.})],$

where ${\cal F}_{0}$ collects all the terms that are independent of the
relative phases between the layers and we only include coupling between
adjacent layers. The most general ansatz for order parameters which respect
the requisite symmetries would be

$\displaystyle\psi_{1}=\psi e^{-i\varphi/2},\ \ \psi_{2}=\psi e^{i\varphi/2},$
(26) $\displaystyle\chi_{1}=\chi e^{-i\varphi^{\prime}/2},\ \ \chi_{2}=\chi
e^{i\varphi^{\prime}/2},$

with the amplitudes $\psi$ and $\chi$ real and positive. The free energy (25)
then becomes

$\displaystyle{\cal F}$ $\displaystyle=$ $\displaystyle{\cal
F}_{0}+2B_{0}\psi^{2}\left[-\cos(2\theta)\cos{\varphi}+{\cal
K}\cos({2\varphi)}\right]$ (27) $\displaystyle+$ $\displaystyle
2B^{\prime}_{0}\psi\chi\left[-\cos(\delta\varphi)+{\cal
K}^{\prime}\cos{(2\delta\varphi)}\right],$

where $\delta\varphi=(\varphi-\varphi^{\prime})/2$ and ${\cal
K}^{\prime}=C^{\prime}\psi\chi/B_{0}^{\prime}$. Observe that terms describing
interlayer and intralayer phase differences, $\varphi$ and $\delta\varphi$
respectively, can be minimized independently. The intralayer free energy,
represented by the second line of Eq. (27), will be dominated by the single
pair tunneling term $-\cos(\delta\varphi)$ on the account of the two CuO2
planes being in a natural untwisted configuration. This implies
$\delta\varphi=0$: the two order parameters within each monolayer will be in
phase, as expected. The remaining interlayer component of the free energy then
coincides with Eq. (I) and its analysis proceeds exactly as before. This
confirms that at the level of the GL theory, the presence of the extra layers
has no effect on our prior conclusions.

Figure 11: Physics of systems with multiple CuO2 planes per monolayer. a) The
structure and notation used for the case with $N=2$ CuO2 planes. Panels b) and
c) show phase diagrams of the lattice BdG model for bilayers with $N=2$,
relevant to the Bi2212 crystal structure with the intra-bilayer coupling set
to $t_{z}=40$meV.

Microscopic models can likewise be extended to the $N>1$ case by adding
appropriate CuO2 planes to each monolayer. We illustrate this by considering a
straightforward generalization of our lattice Hamiltonian Eq. (12) to $N=2$.
This is achieved by adding another square lattice of sites, which is described
by the same Hamiltonian, to each monolayer and coupling the partner lattice
sites by an intralayer tunneling term with amplitude $t_{z}$. Results of the
$N=2$ self-consistent calculation are shown in Fig. 11(b,c) where we chose
$t_{z}=40$meV as appropriate for Bi2212 [53]. We observe physics very similar
to the $N=1$ case, except that topological phases now exhibit Chern numbers 8
and 6, consistent with the notion that in the effective $d+id^{\prime}$ state
each CuO2 plane contributes $C=2$ to the aggregate Chern number of the system.

###  V.6 Josephson tunneling

In the past, several works have studied the problem of Cooper pair tunneling
across $c$-axis twist junctions in order to understand the order parameter
symmetry of cuprates [24, 54]. Assuming a purely $d_{x^{2}-y^{2}}$ order
parameter, it can be shown that the critical Josephson current must vanish
when the twist angle is $45^{\circ}$. Recent experiments in twisted ultra-thin
layers of Bi2212, however, observed that the critical current is essentially
independent of the twist [47]. Here, we show the current is indeed non-
vanishing when ${\cal T}$ is broken.

For the Josephson junction formed by the bilayer, the phase dependent
supercurrent for a given $c$ axis twist can be determined within the framework
of the phenomenological GL theory via the relation

$I(\varphi,\theta)=\frac{2e}{\hbar}\frac{\partial{\cal
F}(\varphi,\theta)}{\partial\varphi},$ (28)

where the free energy has the form given in Eq. (4). A typical current-phase
relation thus obtained is plotted in Fig. 12(a). The critical Josephson
current across the two layers for a given twist configuration,
$I_{c}(\theta)$, may now be obtained by maximizing (28) with respect to
$\varphi$.

As discussed in the main text, for twists in the vicinity of $45^{\circ}$, a
${\cal K}>0$ ushers in a nonzero phase difference between the order parameters
in the static case. This non-trivial character is manifest also in the
behaviour of $I_{c}(\theta)$, as shown in Fig. 12(b). To elucidate this, we
note that setting ${\cal K}=0$ reproduces the behaviour of a junction formed
by superconductors with pure $d$-wave character [24, 54]. Note that
$I_{c}(45^{\circ})=0$. For a positive ${\cal K}$, however, the critical
current at $45^{\circ}$ is non-vanishing. Given that rotations by $\theta$ and
$\pi/2-\theta$ are equivalent from the perspective of the bilayer, the
functional form of the current respects this symmetry and so does its global
maximum. The cusps at $\theta=45^{\circ}$ can be seen as a consequence of
this.

Figure 12: (a) The current-phase relation in units of the critical current at
twist $\theta=41.4^{\circ}$ when ${\cal K}=0.125$. (b) Twist angle dependence
of the critical Josephson current normalized by the current in the untwisted
case. While the two curves show similar behaviour close to
$\theta=0^{\circ}\text{ and }90^{\circ}$, a non-negative ${\cal K}$ results in
a finite value of the current at $\theta=45^{\circ}$.

###  V.7 Density functional theory

We solved for the total energy of the twisted bilayer unit in DFT using the
meta-GGA SCAN exchange correlation potential [36], as implemented in VASP [37,
38]. Projector-Augmented-Wave pseudopotentials were used for the calculations.
Results were converged for a plane-wave energy cutoff of 650 eV on a
(5$\times$5$\times$1) mesh. Similar results to those plotted in Fig. 4 of the
main text were found in the generalized-gradient approximation [55], with
modestly a increased optimal interlayer spacing. For the GGA calculations, we
added Grimme’s semi-empirical DFT-D2 approximation [56] to better account for
the long-range van der Waals binding between the bilayers.

In the main text, we made an effort to parameterize the interlayer coupling by
evaluating the strength of hybridization between Cu $3d$ bands at the unfolded
zone edge, or $X$ point. This momentum point was chosen as previous studies of
interlayer coupling in cuprates have indicated its vanishing strength along
the $\Gamma M$ line [41]. The $X$ point is then a suitable high-symmetry point
where an interlayer-coupling driven hybridization may be studied. Despite the
high density of bands, we have identified a suitable band splitting which is
stable and monotonic over a sufficiently wide range of interlayer spacings, so
that an approximate interlayer coupling evolution can be extracted. The
results, for the unshifted-unit cell, are summarized in Fig. 4(d) of the main
text, and reproduced in Fig. 13(b).

The recently developed SCAN meta-GGA presents the opportunity to perform DFT
on this system without introducing the empirical vdW corrections required for
standard GGA. It was developed with the intention of modelling mixed-bonding
systems [36]. SCAN has been successfully applied to calculations of both LSCO
[39] and YBCO [40]. To confirm the reliability of our results, and their
independence of the exchange-correlation potential, we compare the results of
our interlayer coupling calculations, as computed using standard GGA (with the
Grimme correction), against the SCAN meta-GGA. For simplicity, we use the
unshifted unit cell of Fig. 4(b). A summary of the cohesive energy, as well as
the interlayer coupling, is provided in Fig. 13. Both SCAN and the vdW-
corrected GGA predict a total energy minimum at 12.87 A, albeit with different
cohesive energies.

Figure 13: Bilayer in Density Functional Theory: In panel (a), comparison of
cohesive energy of twisted bilayer Bi2201, as computed in both the GGA and
SCAN meta-GGA. Both exchange correlation potentials recover the same
equilibrium spacing. (b) shows the $X$-point splitting, computed as in LDA
described in the main text.

The similarities are even closer when comparing the X-point splitting,
presented in Fig. 13(b). The splitting decays inversely with the Cu-Cu
spacing, and the equilibrium splitting are nearly indistinguishable. The
agreement between these results indicates the reliability of our $X$-point
metric as an indicator of the interlayer coupling.