Introduction

The charge-4e/6e superconductivities (SCs) are exotic SCs characterized by \(\frac{1}{2}\)/\(\frac{1}{3}\) flux quantization. These novel SCs are formed by condensation of electron quartets/sextets1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25, which is beyond the conventional Bardeen–Cooper–Schrieffer mechanism26. Recently, it was proposed that these intriguing SCs can emerge as the high-temperature vestigial phases of the charge-2e SC in systems hosting multiple coexisting pairing order parameters (ODPs). Typical proposals for such multi-component pairings include the incommensurate pair-density-wave (PDW)10,11,14, the nematic pairing17,18, and the bilayer pairing system16,19. However, each proposal is still waiting for the experiment’s realization.

One proposal is through melting of incommensurate PDW10,11,14. The PDW has been reported in such materials as the cuprates27,28, the CsV3Sb529, and the transition-metal dichalcogenide30. This proposal, however, suffers from the difficulty that the PDWs observed in these experiments are always accompanied by a dominant uniform SC part. Another proposal is through the melting of nematic pairing17,18. Such a pairing state is formed through the real mixing of the two basis functions of a two-dimensional (2D) irreducible representation (IRRP) of the point group. More recently, a group-theory based classification of the vestigial phases generated by melting of the pairing states belonging to the 2D IRRPs was performed31, wherein such interesting phase as d-wave charge-4e SC was proposed. However, the experiment verification of these proposals is still on the way. Alternatively, a bilayer approach was recently proposed16,19 in which, two monolayers hosting SCs with different phase stiffness are coupled. Consequently, in an intermediate-temperature vestigial phase, one layer carries charge-2e SC while the other layer carries charge-4e SC19. The drawback of this proposal lies in that, in an out-of-plane magnetic field, while the charge-4e-SC layer allows for integer times of half magnetic flux, the charge-2e-SC layer only allows for integer flux. As the two layers experience the same magnetic flux, only the integer flux is allowed, and the hallmark of the charge-4e SC, i.e., the half flux quantization, cannot be experimentally detected in this proposal. Finally, the melting of the multi-component hexatic chiral superconductor leading to vestigial charge-6e SC was proposed in the context of kagome superconductors22. Presently, the material realization of the charge-4e/6e SC is still a big challenge.

Here in this work, we take advantage of the rapid development of the twistronics32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,61 and utilize it to design the intriguing charge-4e SC. Here we shall study materials made through stacking two identical monolayers with the largest twist angle, which host Moireless quasi-crystal (QC) structures62,63,64 and are dubbed as the twist-bilayer QC (TB-QC)65, exampled by the recently synthesized 30∘-twisted bilayer graphene66,67,68,69,70 and 45°-twisted bilayer cuprates71,72. Prominently, the TB-QC hosts a doubly enlarged fold of rotation axis relative to its monolayer. Previous study65,73 suggests that when each monolayer hosts a pairing state carrying the largest pairing angular momentum for the lattice, the second-order interlayer Josephson coupling (IJC) between the pairing ODPs from the two layers in the TB-QC makes them mix as 1: ± i, leading to time-reversal symmetry (TRS) breaking chiral topological SC (TSC). For example, as the monolayer cuprate carries the d-wave pairing, the 45°-twisted bilayer cuprates will host the d + id chiral TSC65,74,75,76,77. It’s interesting to investigate possible vestigial secondary orders above the Tc of these chiral TSC phases, driven by the second-order IJC between the pairing ODPs from the two layers.

In this paper, we study the secondary orders in the superconducting TB-QC. Its unique symmetry leads to a simplified low-energy effective Hamiltonian including decoupled total- and relative-phase fields between the bilayer. Significantly, the second-order IJC allows the relative phase to fluctuate between its two saddle points to restore the TRS. Consequently, while the unilateral order of the relative-phase field leads to the TRS-breaking chiral-metal phase, the unilateral quasi-order of the total-phase field leads to the charge-4e SC phase, in which two Cooper pairs from different layers pair to form quartets. These two vestigial phases occupy different regimes in the phase diagram obtained by our combined renormalization group (RG) and Monte-Carlo (MC) studies, which are unambiguously identified by various temperature-dependent quantities including the specific heat, the secondary ODPs, and their susceptibilities, as well as the spatial-dependent correlation functions.

Results

Model and symmetry

Taking two Dn-symmetric monolayers, let’s stack them by the twist angle π/n to form a TB-QC, as shown in Fig. 1 for n = 4 (e.g., the cuprates) and n = 6 (e.g., the graphene). Obviously, the point group is Dnd, isomorphic to D2n. There is an additional symmetry generator in the TB-QC which is absent in its monolayer, i.e., the \({C}_{2n}^{1}\) rotation accompanied by a succeeding layer exchange, renamed as \({\tilde{C}}_{2n}^{1}\) here.

Fig. 1: Schematic illustration of a TB-QC formed by two Dn-symmetric monolayers, with each monolayer carrying SC with pairing angular momentum \(L=\frac{n}{2}\).
figure 1

As examples, n = 4 (cuprates) and n = 6 (graphene) are shown in (a) and (b), respectively.

Full size image

Suppose that driven by some pairing mechanism, the monolayer μ = t/b (top/bottom) can host a pairing state with pairing angular momentum L = n/2. While the cuprate monolayer hosting the d-wave SC synthesized recently78 provides a good example for n = 4, some members in the graphene family which were predicted to host the f-wave SC73,79,80 set an example for n = 6. The pairing gap function in the μ layer is

$${{{\Delta }}}^{(\mu )}({{{{{{{\bf{k}}}}}}}})={\psi }_{\mu }{{{\Gamma }}}^{(\mu )}({{{{{{{\bf{k}}}}}}}}).$$
(1)

Here Γ(μ)(k) is the normalized real form factor, and ψμ is the “complex pairing amplitude”. Prominently, the Γ(μ)(k) for L = n/2 changes sign with every \({C}_{n}^{1}\) rotation. As shown in Fig. 1, we choose a gauge so that

$${{{\Gamma }}}^{({{{{{{{\rm{b}}}}}}}})}({{{{{{{\bf{k}}}}}}}})={\hat{P}}_{\frac{\pi }{n}}{{{\Gamma }}}^{({{{{{{{\rm{t}}}}}}}})}({{{{{{{\bf{k}}}}}}}}),\,\,{\hat{P}}_{\frac{2\pi }{n}}{{{\Gamma }}}^{(\mu )}({{{{{{{\bf{k}}}}}}}})=-{{{\Gamma }}}^{(\mu )}({{{{{{{\bf{k}}}}}}}}).$$
(2)

Here \({\hat{P}}_{\phi }\) indicates the rotation by the angle ϕ. As the interlayer coupling in the TB-QC is weak62,63,65, we can only consider the dominant intralayer pairing, but the two intralayer pairing ODPs can couple through the IJC73,74,75,76,77,80. We shall investigate the ground state and the vestigial secondary orders induced by this IJC.

Firstly, let’s make a saddle-point analysis for the Ginzburg–Landau (G–L) free energy F as functional of ψt/b. For the saddle-point solution, the ψt/b are spatially uniform constant numbers. F is decomposed as,

$$F\left({\psi }_{{{{{{{{\rm{t}}}}}}}}},\;{\psi }_{{{{{{{{\rm{b}}}}}}}}}\right)={F}_{0}\left({\left|{\psi }_{{{{{{{{\rm{t}}}}}}}}}\right|}^{2}\right)+{F}_{0}\left({\left|{\psi }_{{{{{{{{\rm{b}}}}}}}}}\right|}^{2}\right)+{F}_{J}\left({\psi }_{{{{{{{{\rm{t}}}}}}}}},\;{\psi }_{{{{{{{{\rm{b}}}}}}}}}\right),$$
(3)

where \({F}_{0}({|{\psi }_{\mu }|}^{2})\) are the monolayers terms and FJ is the IJC. The TRS-allowed first-order IJC takes the form,

$${F}_{J}^{(1)}\left({\psi }_{{{{{{{{\rm{t}}}}}}}}},\;{\psi }_{{{{{{{{\rm{b}}}}}}}}}\right)=-\alpha \left({\psi }_{{{{{{{{\rm{t}}}}}}}}}{\psi }_{{{{{{{{\rm{b}}}}}}}}}^{*}+c.c\right).$$
(4)

Under \({\tilde{C}}_{2n}^{1}\), the gap function on the μ layer changes from \({{{\Delta }}}^{(\mu )}({{{{{{{\bf{k}}}}}}}})={\psi }_{\mu }{{{\Gamma }}}^{\left(\mu \right)}({{{{{{{\bf{k}}}}}}}})\) to \({\tilde{{{\Delta }}}}^{(\mu )}({{{{{{{\bf{k}}}}}}}})={\psi }_{\bar{\mu }}{\hat{P}}_{\frac{\pi }{n}}{{{\Gamma }}}^{\left(\bar{\mu }\right)}({{{{{{{\bf{k}}}}}}}})\) which, under Eq. (2), can be rewritten as \({\tilde{\psi }}_{\mu }{{{\Gamma }}}^{\left(\mu \right)}({{{{{{{\bf{k}}}}}}}})\) with

$$\tilde{{\psi }_{{{{{{{{\rm{b}}}}}}}}}}={\psi }_{{{{{{{{\rm{t}}}}}}}}},\;\tilde{{\psi }_{{{{{{{{\rm{t}}}}}}}}}}=-{\psi }_{{{{{{{{\rm{b}}}}}}}}}.$$
(5)

The invariance of F under \({\tilde{C}}_{2n}^{1}\) requires α = 0. Thus, the following second-order IJC should be considered,

$${F}_{J}\left({\psi }_{{{{{{{{\rm{t}}}}}}}}},\;{\psi }_{{{{{{{{\rm{b}}}}}}}}}\right)={A}_{0}\left({\psi }_{{{{{{{{\rm{t}}}}}}}}}^{2}{\psi }_{{{{{{{{\rm{b}}}}}}}}}^{2*}+{{{{{{{\rm{c.c.}}}}}}}}\right)+O\left({\psi }^{6}\right).$$
(6)

Equation (6) is minimized at ψb = ±iψt for A0 > 0 or ψb = ±ψt for A0 < 0. Previous microscopic calculations favor the former for the 45°-twisted bilayer cuprates65,74 and 30°-twisted bilayer of the graphene family73,80, leading to d + id or f + if chiral TSCs ground state.

Secondly, let us provide the low-energy effective Hamiltonian for the pairing-phase fluctuations. In this study, we fix Γ(μ) and set \({\psi }_{\mu }\to {\psi }_{\mu }\left({{{{{{{\bf{r}}}}}}}}\right)\) as a slowly varying “envelope” function to describe the spatial fluctuation of the complex pairing amplitude. Focusing on the phase fluctuation, ψt/b are written as \({\psi }_{{{{{{{{\rm{t/b}}}}}}}}}={\psi }_{0}{e}^{i{\theta }_{{{{{{{{\rm{t/b}}}}}}}}}\left({{{{{{{\bf{r}}}}}}}}\right)}\) where ψ0 > 0 is a constant. The \({\theta }_{{{{{{{{\rm{t/b}}}}}}}}}\left({{{{{{{\bf{r}}}}}}}}\right)\) are further written as

$${\theta }_{{{{{{{{\rm{t}}}}}}}}}\left({{{{{{{\bf{r}}}}}}}}\right)={\theta }_{+}\left({{{{{{{\bf{r}}}}}}}}\right)+{\theta }_{-}\left({{{{{{{\bf{r}}}}}}}}\right),\;{\theta }_{{{{{{{{\rm{b}}}}}}}}}\left({{{{{{{\bf{r}}}}}}}}\right)={\theta }_{+}\left({{{{{{{\bf{r}}}}}}}}\right)-{\theta }_{-}\left({{{{{{{\bf{r}}}}}}}}\right).$$
(7)

Here \({\theta }_{+}\left({{{{{{{\bf{r}}}}}}}}\right)\) and \({\theta }_{-}\left({{{{{{{\bf{r}}}}}}}}\right)\) denote the total and relative pairing phases. The low-energy effective Hamiltonian reads

$$H={H}_{0}\left[{\partial }_{\pm }{\theta }_{+},\;{\partial }_{\pm }{\theta }_{-}\right]+{A}_{0}{\psi }_{0}^{4}\int\cos 4{\theta }_{-}\left({{{{{{{\bf{r}}}}}}}}\right){d}^{2}{{{{{{{\bf{r}}}}}}}},$$
(8)

with ∂± ≡ ∂x ± i∂y. Up to the lowest-order expansion, the H0 takes the following explicit form in the k-space,

$${H}_{0}= \frac{1}{2}\int{d}^{2}{{{{{{{\bf{k}}}}}}}}\left[{\theta }_{+}\left({{{{{{{\bf{k}}}}}}}}\right){\theta }_{+}\left({{{{{{{\bf{-k}}}}}}}}\right)\left(\alpha {k}_{+}^{2}+\beta {k}_{-}^{2}+\rho {k}_{+}{k}_{-}\right)\right.\\ +{\theta }_{+}\left({{{{{{{\bf{k}}}}}}}}\right){\theta }_{-}\left({{{{{{{\bf{-k}}}}}}}}\right)\left(\omega {k}_{+}^{2}+\delta {k}_{-}^{2}+\eta {k}_{+}{k}_{-}\right)\\ +\left.{\theta }_{-}\left({{{{{{{\bf{k}}}}}}}}\right){\theta }_{-}\left({{{{{{{\bf{-k}}}}}}}}\right)\left(\epsilon {k}_{+}^{2}+\xi {k}_{-}^{2}+\kappa {k}_{+}{k}_{-}\right)\right].$$
(9)

Under \({\tilde{C}}_{2n}^{1}\), the gap function on the μ layer changes from \({{{\Delta }}}^{(\mu )}={\psi }_{\mu }{{{\Gamma }}}^{\left(\mu \right)}\) to \({\tilde{{{\Delta }}}}^{(\mu )}={\psi }_{\bar{\mu }}{\hat{P}}_{\frac{\pi }{n}}{{{\Gamma }}}^{\left(\bar{\mu }\right)}\) which, under Eq. (2), can be rewritten as \({\tilde{\psi }}_{\mu }{{{\Gamma }}}^{\left(\mu \right)}\) with

$${\tilde{\psi }}_{{{{{{{{\rm{b}}}}}}}}}\left({{{{{{{\bf{r}}}}}}}}\right)={\psi }_{{{{{{{{\rm{t}}}}}}}}}\left({\hat{P}}_{\frac{\pi }{n}}^{-1}{{{{{{{\bf{r}}}}}}}}\right),\;{\tilde{\psi }}_{{{{{{{{\rm{t}}}}}}}}}\left({{{{{{{\bf{r}}}}}}}}\right)=-{\psi }_{{{{{{{{\rm{b}}}}}}}}}\left({\hat{P}}_{\frac{\pi }{n}}^{-1}{{{{{{{\bf{r}}}}}}}}\right).$$
(10)

Consequently, we have

$${\theta }_{+}\left({{{{{{{\bf{k}}}}}}}}\right)\to {\theta }_{+}\left({\hat{P}}_{\frac{\pi }{n}}^{-1}{{{{{{{\bf{k}}}}}}}}\right),\,{\theta }_{-}\left({{{{{{{\bf{k}}}}}}}}\right)\to -{\theta }_{-}\left({\hat{P}}_{\frac{\pi }{n}}^{-1}{{{{{{{\bf{k}}}}}}}}\right).$$
(11)

The invariance of Eq. (9) under (11) only allows for nonzero ρ and κ, leading to the real-space Hamiltonian

$$H=\int{d}^{2}{{{{{{{\bf{r}}}}}}}}\left(\frac{\rho }{2}{\left|\nabla {\theta }_{+}\right|}^{2}+\frac{\kappa }{2}{\left|\nabla {\theta }_{-}\right|}^{2}+{A}_{0}{\psi }_{0}^{4}\cos 4{\theta }_{-}\right).$$
(12)

Equation (12) shows two important features. Firstly, the θ+ and θ− fields are dynamically decoupled, with each hosting different stiffness parameter ρ or κ derived by the G-L expansion in the Sec. I of Supplementary Information (SI). Secondly, the second-order IJC allows θt − θb = 2θ− to fluctuate between its two saddle points, i.e., ±π/2, to restore the TRS. Note that although the term \(\cos (4{\theta }_{-})\) in Eq. (12) leads to four different values of θ−: ±π/4 and ± 3π/4 for the ground state, θ− = π/4 (−π/4) leads to gauge equivalent state with θ− = −3π/4 (3π/4). So the system only possesses two-fold Ising anisotropy. Here the unilateral quasi-ordering of the θ+ field leads to the ODP Δ(t)(k) ⋅ Δ(b)(−k) characterizing the charge-4e SC in which two Cooper pairs from different layers pair. The unilateral ordering of the θ− field leads to the ODP Δ(t)*(k) ⋅ Δ(b)(k) characterizing the TRS breaking chiral metal81,82. Note that while θ+ and θ− each can host either integer or half-integer vortices, Eq. (7) requires that they can only simultaneously host integer or half-integer vortices to ensure the single-valuedness of ψt/b10,18. This sets the “kinematic constraint” in the low-energy “classical Hilbert space” for allowed vortices of the two fields.

RG study: To perform the RG study, we start with the following effective action at the temperature T,

$$S=\int{d}^{2}{{{{{{{\bf{r}}}}}}}}\left(\frac{\rho }{2T}{\left|\nabla {\theta }_{+}\right|}^{2}+\frac{\kappa }{2T}{\left|\nabla {\theta }_{-}\right|}^{2}+{g}_{4}\cos 4{\theta }_{-}\right)$$
(13)

Here g4 > 0 is proportional to A0. This action can be mapped to a two-component Sine-Gordon model,

$${S}_{{{{{{{{\rm{SG}}}}}}}}}= \int{d}^{2}{{{{{{{\bf{x}}}}}}}}\left(\frac{T}{2\rho }{\left|\nabla {\tilde{\theta }}_{+}\right|}^{2}+\frac{T}{2\kappa }{\left|\nabla {\tilde{\theta }}_{-}\right|}^{2}+{g}_{4}\cos 4{\theta }_{-}-{g}_{2,0}\right.\\ \left.\times \cos 2\pi {\tilde{\theta }}_{+}-{g}_{0,2}\cos 2\pi {\tilde{\theta }}_{-}-{g}_{1,1}\cos \pi {\tilde{\theta }}_{+}\cos \pi {\tilde{\theta }}_{-}\right)$$
(14)

The dual bosonic fields \({\tilde{\theta }}_{+}\) and \({\tilde{\theta }}_{-}\) describe the vortices of the fields θ+ and θ−. g2,0, g0,2, and g1,1 are coupling parameters proportional to the fugacities of different types of vortices (g2,0/g0,2: integer vortices; g1,1: half vortices).

The phase diagram obtained by the one-loop RG analysis provided in “Methods” is shown in Fig. 2a. Variation of the initial coupling parameters does not change the topology of the phase diagram, which always includes the chiral TSC, charge-4e SC, chiral metal, and normal metal phases, see the Section IV of SI. At low enough T, the vortex fugacities g2,0, g0,2, and g1,1 are all irrelevant while the IJC parameter g4 is relevant, suggesting that both the θ± fields are locked, leading to the TRS breaking chiral SC. With the enhancement of T, in the low κ/ρ regime, g0,2 first gets relevant (and suppresses g4) suggesting that the θ− vortices proliferate to restore the TRS, to form the charge-4e SC. In the high κ/ρ regime instead, the g2,0 first gets relevant suggesting the θ+ vortices proliferate to kill the SC, to form the chiral metal. In both regimes, at high enough T, g2,0, and g0,2 are both relevant, forming the normal metal phase. In the regime κ ≈ ρ, with the enhancement of T, the system transits into a phase wherein the coupling g1,1 is relevant and the half vortices involving both fields proliferate to kill both (quasi) orders, suggesting that the system directly transit to the normal state.

Fig. 2: The phase diagram.
figure 2

Phase diagram provided by a the RG study and b the MC study. The initial values of the coupling parameters in a are g2,0 = g0,2 = 0.1, g1,1 = g4 = 0.01 in Eq. (14) and in (b) are A = 0.025ρ and \(\gamma=\frac{1}{4}\rho \kappa /(\rho+\kappa )\) in Eq. (15). The three dotted lines and four dots (A–D) in (b) are used for subsequent interpretation of the phase diagram.

Full size image

In the charge-4e SC, the Josephson-coupling phase, i.e., θ−, is disordered. However, this phase should not be understood as a layer-decoupled charge-2e SC from each layer, as in this phase the pairing phase of each layer is also disordered. To remind you, the charge-4e SC proposed here only lives in the intermediate temperature above the Tc of the pairing state, wherein each layer is no longer superconducting. In the chiral metal phase, the time-reversal symmetry breaking can be verified by the polar Kerr effect. Furthermore, there can be spontaneously generated inner magnetic field in the material, which can be detected by the muon spin resonance experiment.

MC study

To perform the MC study, we discretize the Hamiltonian (12) on the square lattice to obtain

$$H= -\alpha \mathop{\sum}\limits_{\langle ij\rangle }\cos [{\theta }_{{{{{{{{\rm{t}}}}}}}}}({{{{{{{{\bf{r}}}}}}}}}_{i})+{\theta }_{{{{{{{{\rm{b}}}}}}}}}({{{{{{{{\bf{r}}}}}}}}}_{i})-{\theta }_{{{{{{{{\rm{t}}}}}}}}}({{{{{{{{\bf{r}}}}}}}}}_{j})-{\theta }_{{{{{{{{\rm{b}}}}}}}}}({{{{{{{{\bf{r}}}}}}}}}_{j})]\\ -\lambda \mathop{\sum}\limits_{\langle ij\rangle }\cos [{\theta }_{{{{{{{{\rm{t}}}}}}}}}({{{{{{{{\bf{r}}}}}}}}}_{i})-{\theta }_{{{{{{{{\rm{b}}}}}}}}}({{{{{{{{\bf{r}}}}}}}}}_{i})-{\theta }_{{{{{{{{\rm{t}}}}}}}}}({{{{{{{{\bf{r}}}}}}}}}_{j})+{\theta }_{{{{{{{{\rm{b}}}}}}}}}({{{{{{{{\bf{r}}}}}}}}}_{j})]\\ -\gamma \mathop{\sum}\limits_{\langle ij\rangle }\cos [{\theta }_{{{{{{{{\rm{t}}}}}}}}}({{{{{{{{\bf{r}}}}}}}}}_{i})-{\theta }_{{{{{{{{\rm{t}}}}}}}}}({{{{{{{{\bf{r}}}}}}}}}_{j})]+\cos [{\theta }_{{{{{{{{\rm{b}}}}}}}}}({{{{{{{{\bf{r}}}}}}}}}_{i})-{\theta }_{{{{{{{{\rm{b}}}}}}}}}({{{{{{{{\bf{r}}}}}}}}}_{j})]\\ +A\mathop{\sum}\limits_{i}\cos [2{\theta }_{{{{{{{{\rm{t}}}}}}}}}({{{{{{{{\bf{r}}}}}}}}}_{i})-2{\theta }_{{{{{{{{\rm{b}}}}}}}}}({{{{{{{{\bf{r}}}}}}}}}_{i})].$$
(15)

Here 〈ij〉 represents nearest-neighbor bonding, and the positive coefficients α, λ, and γ satisfy,

$$\alpha=\frac{\rho -2\gamma }{4},\;\lambda=\frac{\kappa -2\gamma }{4}.$$
(16)

Note that although different α, λ, and γ satisfying Eq. (16) leads to the same continuous Hamiltonian (12) in the continuum limit, it is required that all of them should be positive so as to reproduce the correct low-energy “classical Hilbert space” for allowed vortices. The reason is as follows. Here the α > 0 and λ > 0 terms energetically allow for integer or half-integer θ+ and θ− vortices, while the γ > 0 term energetically only allows for integer θt or θb vortices and hence imposes the “kinematic constraint” between the θ+ and θ− vortices. Note that although the γ term does not naturally emerge from Eq. (12), the singlevaluedness of the ψt/b field dictates it. This term is crucial to yield the correct topology of the phase diagram. As shown in the Sec. VI of SI, if we turn off the γ term, θ+ and θ− are decoupled, leading to a topologically wrong phase diagram. A comparison between the correct phase diagram and the wrong one shows that the kinematic correlation makes the vestigial phase regimes largely shrink. For thermodynamic limit, even an infinitesimal γ can energetically guarantee the ”kinematic constraint”. Here in the discrete lattice, we set \(\gamma=\frac{1}{4}\rho \kappa /(\rho+\kappa ),\,A=0.025\rho\), and their other values lead to similar results, see Section VI of SI.

The MC phase diagram shown in Fig. 2b is qualitatively consistent with the RG one shown in Fig. 2a. Various T-dependent quantities are shown in Fig. 3 for κ/ρ = 0.3, 1, 2.2 marked in Fig. 2b, with the formulas adopted in the MC calculations provided in “Methods”. For κ/ρ = 0.3, the specific heat Cv is shown in Fig. 3a, where the high-T broad hump characterizes the Kosterlitz–Thouless (K–T) phase transition between the normal state and the charge-4e SC and the low-T sharp peak characterizes the Ising phase transition between the charge-4e SC and the chiral SC. For this κ/ρ, Fig. 3b shows the phase stiffness S characterizing the SC and the Ising ODP I characterizing the relative-phase order16, which emerge at the critical temperatures corresponding to the broad hump and sharp peak in Fig. 3a, respectively. Furthermore, the total- (χ+) and relative- (χ−) phase susceptibilities16 shown in Fig. 3c diverge at the same critical temperatures. For κ/ρ = 1, the specific heat shown in Fig. 3d exhibits only one peak, suggesting a direct phase transition from the normal state to the chiral SC. Such a result is also reflected in Fig. 3e, f which shows that the total- and relative-phase (quasi) orders emerge at the same temperature. For κ/ρ = 2.2, the corresponding results shown in Fig. 3g–i reveal that following the decrease of T, the system will successively experience the normal state, the chiral metal, and the chiral TSC phases. The results presented in Fig. 3 are well consistent with the phase diagram shown in Fig. 2b.

Fig. 3: Various T-dependent quantities.
figure 3

Various T-dependent quantities for κ = 0.3 (a–c), κ = 1 (d–f), and κ = 2.2 (g–i). a, d, g The specific heat Cv. b, e, h The phase stiffness S (blue) and Ising ODP I (red). c, f, i The susceptibilities χ+ (blue) and χ− (red). The ρ is set as the unit of κ and T.

Full size image

The total- (+) and relative- (−) phase correlation functions η± are shown in Fig. 4. See their formulas in “Methods”. Figure 4a, b shows that for the representative point A marked in Fig. 2b, while η+(Δr) power-law decays with Δr suggesting quasi-long-range order of the total phase, η−(Δr) decays exponentially with Δr, suggesting disorder of the relative phase. Obviously, these electron correlations are consistent with the charge-4e SC phase. Figure 4c, d shows that for the point D, while η+(Δr) decays exponentially with Δr suggesting disorder of the total phase, η−(Δr) saturates to a constant number for large enough Δr suggesting long-range order of the relative phase, consistent with the chiral-metal phase. For comparison, the η± for points B and C provided in Section V of SI is also consistent with the normal-metal and chiral-SC phases.

Fig. 4: The correlation functions.
figure 4

The correlation function η± for a and b for A(κ = 0.2, T = 0.2), and for c and d for D(κ = 2.2, T = 0.5) marked in Fig. 2b. Insets of a the log–log plot, and b and c only the y-axes are logarithmic.

Full size image

Discussions

In comparison with previous proposals for the charge-4e/6e SC based on melting of the PDW10,11,14 or the nematic pairing17,18, our proposal is based on a more definite and easily realized start point: here we only need to start from non-topological d-wave SC (or f-wave SC) in any fourfold (or sixfold) symmetric monolayers. Particularly, we have provided concrete synthesized materials to realize our proposal, i.e., the 45o-twisted bilayer cuprates and the 30o-twisted bilayer of some graphene family. Furthermore, superior to the previous bilayer approach, here a Cooper pair from the top layer pairs with a Cooper pair from the bottom layer to form the charge-4e SC between the layers. Consequently, the half flux quantization can be experimentally detected as a hallmark of the charge-4e SC in our proposal.

The TB-QC provides a better platform to realize the vestigial phases than conventional chiral superconductors such as the p + ip or d + id ones on the square or honeycomb lattices. The latter also hosts two degenerate pairing ODPs and hence can accommodate both total and relative-phase fluctuations of the two ODPs. However, the rotational symmetry of the monolayer system is not as high as that of the TB-QC studied here. Consequently, for chiral TSC in monolayer systems, there can be many nonzero coefficients in Eq. (9). Particularly, the two-phase fields are generally dynamically coupled as the symmetries in these systems allow for extra terms such as ∇±θ+ ⋅ ∇±θ− in the Hamiltonian density in Eq. (12). See more details in Section II of SI. As shown in Fig. 2 and Fig. S5a, the kinematic correlation between θ+ and θ− has already made the vestigial phase regimes largely shrink, their extra dynamic coupling might make them further shrink or even vanish.

In conclusion, we have predicted the realization of the charge-4e SC or the chiral metal in the TB-QC, emerging as the unilateral (quasi) ordering of the total- or relative- pairing phase of the two layers, above the chiral-TSC ground state. The TB-QC provides a better platform to realize these vestigial phases than previous proposals as here we can start from a more definite and easily realized start point.

Methods

The RG approach

Here we provide some technique details for the RG study. With standard RG analysis, the flow equations at the one-loop level are given by:

$$\frac{d{g}_{2,0}}{d\ln b} =(2-\pi {\rho }^{{\prime} }){g}_{2,0}\\ \frac{d{g}_{0,2}}{d\ln b} =(2-\pi {\kappa }^{{\prime} }){g}_{0,2}\\ \frac{d{g}_{1,1}}{d\ln b} =\left(2-\frac{\pi }{4}({\rho }^{{\prime} }+{\kappa }^{{\prime} })\right){g}_{1,1}\\ \frac{d{g}_{4}}{d\ln b} =\left(2-\frac{4}{\pi {\kappa }^{{\prime} }}\right){g}_{4}\\ \frac{d{\rho }^{{\prime} }}{d\ln b} =-16{g}_{2,0}^{2}{\rho }^{{\prime} 3}-\frac{{g}_{1,1}^{2}}{2}{\rho }^{{\prime} 2}({\rho }^{{\prime} }+{\kappa }^{{\prime} })\\ \frac{d{\kappa }^{{\prime} }}{d\ln b} =\frac{256{g}_{4}^{2}}{{\pi }^{4}{\kappa }^{{\prime} }}-16{g}_{0,2}^{2}{\kappa }^{{\prime} 3}-\frac{{g}_{1,1}^{2}}{2}{\kappa }^{{\prime} 2}({\rho }^{{\prime} }+{\kappa }^{{\prime} }),$$
(17)

Here b represents the renormalization scale, g2,0, g0,2, and g1,1 represent the coupling strength of different types of topological defects, \({\rho }^{{\prime} }=\rho /T\) and \({\kappa }^{{\prime} }=\kappa /T\) represent two kinds of stiffness parameters.

The fixed points of N general RG flow equation \(\frac{d{{{{{{{\bf{g}}}}}}}}}{d\ell }=R({{{{{{{\bf{g}}}}}}}})\) is obtained by R(g*) = 0. In Table 1, we present four fixed points of the RG flow equation (17) and the corresponding phases in our calculation. We find the renormalized values of the stiffness parameters \({\rho }^{{\prime} }\) and \({\kappa }^{{\prime} }\) are consistent with the phase revealed by the RG flow result of the g-couplings. Specifically, the \({\rho }^{{\prime} }\) flows to a finite positive value if the U(1)-gauge symmetry is (quasi) broken, otherwise it flows to zero; the \({\kappa }^{{\prime} }\) flows to infinity if the time- reverse symmetry is broken, otherwise it flows to zero. Furthermore, the stability analysis of the fixed points can be provided following the standard process83. We only outline here and see more details in the Section III of SI. The β function of the coupling constant which is very close to the fixed point g* can be replaced by a linear mapping:

$$R({{{{{{{\bf{g}}}}}}}})=R\left(({{{{{{{\bf{g}}}}}}}}-{{{{{{{{\bf{g}}}}}}}}}^{*})+{{{{{{{{\bf{g}}}}}}}}}^{*}\right)\simeq M({{{{{{{\bf{g}}}}}}}}-{{{{{{{{\bf{g}}}}}}}}}^{*})$$
(18)

where we have used R(g*) = 0, and \({M}_{\alpha \beta }=\frac{\partial {R}_{\alpha }}{\partial {g}_{\beta }}{| }_{{{{{{{{\bf{g}}}}}}}}={{{{{{{{\bf{g}}}}}}}}}^{*}}\). To get the stability properties of the flow, we have to diagonalize the matrix MN×N. The eigenvalues denoted by λα, α = 1, 2, . . . , N. If the real parts of all the eigenvalues are negative or, at worst, zero, i.e., the scaling fields are all irrelevant or marginal. There are stable fixed points corresponding to the “stable phases". Complementary to the stable fixed points, if all the eigenvalues are positive and the scaling fields are all relevant, there are unstable fixed points. Additionally, there is a generic class of fixed points with both relevant and irrelevant scaling fields. These points are associated with the boundary of the phase transition.

Table 1 Fixed points of the coupling parameters under RG, and the corresponding phases
Full size table

The Monte-Carlo approach

Here we provide some formulas for the MC calculations.

The phase stiffness characterizes the quasi-long-range order of the total phase and hence the SC is16

$$S=\frac{1}{N}\left( \left \langle {H}_{x} \right \rangle -\beta \left \langle {I}_{x}^{2} \right \rangle \right)$$
(19)

with

$${H}_{x}= 4\alpha \mathop{\sum}\limits_{ < ij{ > }_{x}}\cos [{\theta }_{{{{{{{{\rm{t}}}}}}}}}({{{{{{{{\bf{r}}}}}}}}}_{i})+{\theta }_{{{{{{{{\rm{b}}}}}}}}}({{{{{{{{\bf{r}}}}}}}}}_{i})-{\theta }_{{{{{{{{\rm{t}}}}}}}}}({{{{{{{{\bf{r}}}}}}}}}_{j})+{\theta }_{{{{{{{{\rm{b}}}}}}}}}({{{{{{{{\bf{r}}}}}}}}}_{j})]\\ +\gamma \mathop{\sum}\limits_{ < ij{ > }_{x}}\cos [{\theta }_{{{{{{{{\rm{t}}}}}}}}}({{{{{{{{\bf{r}}}}}}}}}_{i})-{\theta }_{{{{{{{{\rm{t}}}}}}}}}({{{{{{{{\bf{r}}}}}}}}}_{j})]+\cos [{\theta }_{{{{{{{{\rm{b}}}}}}}}}({{{{{{{{\bf{r}}}}}}}}}_{i})-{\theta }_{{{{{{{{\rm{b}}}}}}}}}({{{{{{{{\bf{r}}}}}}}}}_{j})]\\ {I}_{x} = 2\alpha \mathop{\sum}\limits_{ < ij{ > }_{x}}\sin [{\theta }_{{{{{{{{\rm{t}}}}}}}}}({{{{{{{{\bf{r}}}}}}}}}_{i})+{\theta }_{{{{{{{{\rm{b}}}}}}}}}({{{{{{{{\bf{r}}}}}}}}}_{i})-{\theta }_{{{{{{{{\rm{t}}}}}}}}}({{{{{{{{\bf{r}}}}}}}}}_{j})+{\theta }_{{{{{{{{\rm{b}}}}}}}}}({{{{{{{{\bf{r}}}}}}}}}_{j})]\\ +\gamma \mathop{\sum}\limits_{ < ij{ > }_{x}}\sin [{\theta }_{{{{{{{{\rm{t}}}}}}}}}({{{{{{{{\bf{r}}}}}}}}}_{i})-{\theta }_{{{{{{{{\rm{t}}}}}}}}}({{{{{{{{\bf{r}}}}}}}}}_{j})]+\sin [{\theta }_{{{{{{{{\rm{b}}}}}}}}}({{{{{{{{\bf{r}}}}}}}}}_{i})-{\theta }_{{{{{{{{\rm{b}}}}}}}}}({{{{{{{{\bf{r}}}}}}}}}_{j})],$$
(20)

where N is the site number, and β = 1/kBT.

The Ising order parameter characterizing the relative-phase ordering breaking the time-reversal symmetry is,

$$I\equiv \frac{1}{{N}^{2}}\mathop{\sum}\limits_{ij}\left\langle \sin [{\theta }_{{{{{{{{\rm{t}}}}}}}}}({{{{{{{{\bf{r}}}}}}}}}_{i})-{\theta }_{{{{{{{{\rm{b}}}}}}}}}({{{{{{{{\bf{r}}}}}}}}}_{i})]\cdot \sin [{\theta }_{{{{{{{{\rm{t}}}}}}}}}({{{{{{{{\bf{r}}}}}}}}}_{j})-{\theta }_{{{{{{{{\rm{b}}}}}}}}}({{{{{{{{\bf{r}}}}}}}}}_{j})]\right\rangle .$$

The total- (+) and relative- (−) phase susceptibilities for temperatures above the Tc of the corresponding orders are defined by

$${\chi }_{\pm }\equiv \frac{1}{NT}\mathop{\sum}\limits_{i}\left\langle {\left| {e}^{i\left[{\theta }_{{{{{{{{\rm{t}}}}}}}}}\left({{{{{{{{\bf{r}}}}}}}}}_{i}\right)\pm {\theta }_{{{{{{{{\rm{b}}}}}}}}}\left({{{{{{{{\bf{r}}}}}}}}}_{i}\right)\right]}\right| }^{2}\right\rangle .$$
(21)

The total- (+) and relative- (−) phase correlation functions are defined as

$${\eta }_{\pm }({{\Delta }}{{{{{{{\bf{r}}}}}}}})=\frac{1}{N}\mathop{\sum}\limits_{{{{{{{{\bf{r}}}}}}}}}\left\langle {e}^{i[{\theta }_{t}({{{{{{{\bf{r}}}}}}}})\pm {\theta }_{b}({{{{{{{\bf{r}}}}}}}})-{\theta }_{t}({{{{{{{\bf{r}}}}}}}}+{{\Delta }}{{{{{{{\bf{r}}}}}}}})\mp {\theta }_{b}({{{{{{{\bf{r}}}}}}}}+{{\Delta }}{{{{{{{\bf{r}}}}}}}})]}\right\rangle .$$
(22)