Twisted multilayer nodal superconductors

Twisted bilayers of nodal superconductors were recently proposed as a promising platform to host superconducting phases that spontaneously break time-reversal symmetry. Here we extend this analysis to twisted multilayers, focusing on two high-symmetry stackings with alternating ($\pm \theta$) and constant ($\theta$) twist angles. In analogy to alternating-twist multilayer graphene, the former can be mapped to twisted bilayers with renormalized interlayer couplings, along with a remnant gapless monolayer when the number of layers $L$ is odd. In contrast, the latter exhibits physics beyond twisted bilayers, including the occurrence of “magic angles” characterized by cubic band crossings when $L\mathrm{mod}4=3$. Due to their power-law divergent density of states, such multilayers are highly susceptible to secondary instabilities. Within a BCS mean-field theory, defined in the continuum and on a lattice, we find that both stackings host chiral topological superconductivity in extended regions of their phase diagrams.

Figure 1

Spontaneous $\mathcal{T}$ breaking in twisted trilayer nodal superconductors. The top panels (a)–(d) depict data for the alternating twist (AT) stacking, and the bottom panels (e)–(h) depict data for the chiral twist (CT) stacking. (a),(e) Lattice geometries. (a) The AT stacking is obtained by twisting successive layers by $\pm \theta$, and it is invariant under an out-of-plane mirror reflection ${\mathcal{M}}_z$ with respect to the middle layer. (e) The CT stacking is obtained by twisting successive layers by the same angle $\theta$ with respect to a common origin in the plane, and it is invariant under the ${\mathcal{C}}_2$ symmetry that implements a $\pi$ rotation about the diagonals of the middle layers. (b), (f) Contour plot of the BdG free energy (in arbitrary units) for $\theta ={45}^{\circ }$ in the space of the phases ${\phi }_1, {\phi }_3$ of the top and bottom layer order parameters (setting ${\phi }_2=0$). (c), (g) Low-energy spectrum in one quadrant of the BZ for large twist angles. The CT stacking is fully gapped, whereas the AT stacking is gapless due to the Dirac cone dispersion of the decoupled monolayer. (d), (h) Phase differences ${\phi }_1$ and ${\phi }_3$ that minimize the BdG free energy at zero temperature, and the spectral gap (gray curve) as a function of the twist angle. In (d), the physics is similar to that of a standalone bilayer, expect that the BdG spectrum is gapless due to the decoupled Dirac cone. In (h), the nontrivial phases break $\mathcal{T}$ but preserve the product ${\mathcal{C}}_2\mathcal{T}$. The corresponding dispersion is gapped except at ${45}^{\mathrm{o}}$. The parameters are chosen as ${\epsilon }_c=60\mathrm{meV}, |{\mathrm{\Delta }}_l|\approx 40\mathrm{meV}$, and $g=20\mathrm{meV}$.

Figure 2

Spectrum of the AT trilayer (top) and the CT trilayer (bottom) on an infinite cylinder geometry for a commensurate twist angle ${\theta }_{1,2}=36.9^{\circ }$ (see Appendix pp5 for modeling details). The color scale shows the expectation value $\langle \widehat{y}\rangle$ of the eigenstates, with the length $y$ along the open direction normalized to unity. The AT trilayer shows four chiral edge 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 Dirac cones. In contrast, the CT trilayer shows a fully gapped bulk with six chiral edge modes. The unit cells in the two configurations have 15 and 75 sites, respectively. Both calculations use 90 unit cells along the open direction with parameters $t=1, \mu =-1.3, \mathrm{\Delta }=0.6$, and $g=0.5$.

Figure 3

Dispersion of BdG quasiparticles in one quadrant of the Brillouin zone, absent any gap-opening instabilities, as the system is tuned across the magic angle ${\theta }_{\mathrm{M}}$ for trilayers in the AT (top row) and CT (bottom row) configurations. (a) The decoupling of AT trilayers in two sectors is manifest by the merging of two Dirac cones in a quadratic band 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 the diagonals at $\mathit{q}=\mathit{0}$ is protected by the ${\mathcal{C}}_2$ symmetry of the CT stacking and leads to a cubic band touching (CBT) with a divergent density of states $D\left(E\right)\sim E^{-1/3}$ at low energies, as shown in (d).

Figure 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=12\mathrm{meV}$, showing the layer-resolved order-parameter amplitudes $|{\mathrm{\Delta }}_l|$, phase factor ${\phi }_1$, and the spectral gap $E_{\mathrm{gap}}$ as a function of $\theta$. The $\mathcal{T}$-breaking phase is nucleated near the magic angles ${\theta }_{\mathrm{M}}$ (vertical dotted gray lines), where the density of states is maximal, and persists to large twist angles up to $\pi /2-{\theta }_{\mathrm{M}}$. For AT trilayers the spectrum remains gapless throughout due to the decoupled monolayer sector. (b), (e) Phase diagrams showing ${\phi }_1/\pi$ in the $g\text{−}\theta$ plane. The analytical magic angle conditions are indicated by dashed gray lines, which track the onset of the topological regime for small twist angles. The orange and purple lines show the $g=12\mathrm{meV}$ cut used in (c), (f) and the $g=20\mathrm{meV}$ cut used in Fig. 1, respectively. (c), (f) Phase diagrams showing ${\phi }_1/\pi$ in the $T\text{−}\theta$ plane, with the temperature normalized by the critical temperature $T_c$ for zero twist angle. The $\mathcal{T}$-broken domes where topological SC is stabilized around the magic angles are delineated in (f).

Figure 5

Self-consistent mean-field results for quadrilayers in the AT (top) and CT (bottom) configuration, using $g=12\mathrm{meV}, |{\mathrm{\Delta }}_l|\sim 40\mathrm{meV}, {\epsilon }_c=60\mathrm{meV}$ and zero temperature. For twist angles around $\theta =\pi /4$, the SC phase breaks $\mathcal{T}$ for both stackings. Top: The magic angles are determined by the effective tunneling in the bilayer blocks ${\tilde{g}}_{1,2}=(\sqrt{5}\mp 1)g/2$ and are denoted by dashed gray lines. For these parameters, $\mathcal{T}$-breaking sets in only at the second magic angle. Bottom: For CT quadrilayers, chiral SC which preserves the ${\mathcal{C}}_2\mathcal{T}$ symmetry is stabilized for a large range of twist angles.

Figure 6

Phase diagram comparison of chirally twisted multilayers. The color scale indicates the $\mathcal{T}$-breaking gap in (a) bilayers, (b) trilayers, (c) quadrilayers, and (d) pentalayers. The extent of the topological region of the phase diagram increases with $L$, and the maximal gaps are obtained away from $\theta =\pi /4$. We choose parameters ${\epsilon }_c=60\mathrm{meV}, |{\mathrm{\Delta }}_l|\approx 40\mathrm{meV}$, and $T=0$. The gray areas in (c) and (d) denote regions where numerical convergence is difficult due to the small energy scales involved.

Figure 7

Evolution of the BdG quasiparticle dispersion of chirally twisted quadrilayers at small angles. There are two magic angles ${\theta }_{\mathrm{M}}^{\left(1\right)}$ and ${\theta }_{\mathrm{M}}^{\left(2\right)}$ where two quadratic band touchings occur simultaneously at momenta $\pm {\mathit{q}}_{\mathrm{M}}^{\left(1\right)}$ and $\pm {\mathit{q}}_{\mathrm{M}}^{\left(2\right)}$, respectively.

Figure 8

Self-consistent mean-field results for the $s$-wave order parameter amplitude of the middle layer ${\mathrm{\Delta }}_2^s$ (black lines) and the $d$-wave phase factor ${\phi }_1$ (red lines) as a function of twist angle $\theta$ for trilayers in the AT (top panel) and CT (bottom panel) stacking. We consider a secondary $s$-wave pairing channel with relative strength $\mathcal{U}/\mathcal{V}=0.28$, 0.32, and 0.36, and we use parameters $g=12\mathrm{meV}$ and ${\epsilon }_c=60\mathrm{meV}$, while the primary $d$-wave order parameter $\mathrm{\Delta }\sim 40\mathrm{meV}$. For $\mathcal{U}/\mathcal{V}\lesssim 0.30$ we find that the $d+is$ solution is not nucleated, and the $d+i{d}^{\prime }$ solution first develops around the magic angles. For large values of $\mathcal{U}$ we observe a competition between the two gapped $\mathcal{T}$-broken phases, with the $d+i{d}^{\prime }$ solution preferred at large angles and the $d+is$ solution preferred near the magic angles.