Prediction of time-reversal-symmetry breaking fermionic quadrupling condensate in twisted bilayer graphene

Recent mean-field calculations suggest that the superconducting state of twisted bilayer graphene exhibits either a nematic order or a spontaneous breakdown of the time-reversal symmetry. The two-dimensional character of the material and the large critical temperature relative to the Fermi energy dictate that the material should have significant fluctuations. We study the effects of these fluctuations using Monte Carlo simulations. We show that in a model proposed earlier for twisted bilayer graphene there is a fluctuation-induced phase with quadrupling fermionic order for all considered parameters. This four-electron condensate, instead of superconductivity, shows a spontaneous breaking of time-reversal symmetry. Our results suggest that twisted bilayer graphene is an especially promising platform to study different types of condensates, beyond the pair-condensate paradigm.


INTRODUCTION
The recently discovered superconducting state which emerges in magic-angle twisted bilayer graphene exhibits a critical temperature that is exceptionally high compared to the Fermi energy [1][2][3][4][5]. This, and the fact that the system is two-dimensional, implies the presence of strong pairing fluctuations.
While superconductivity is a more than century-old state of matter, which results from electron pairing, the presence of strong fluctuations suggest the tantalizing possibility that magic-angle twisted bilayer graphene can be an especially promising system to realize different states of matter in the form of condensates of electronic quadruplets. In principle, the standard Bardeen-Cooper-Schrieffer theory does not allow fermionic quadrupling condensates. However, if the low-temperature regime of twisted bilayer graphene exhibits a superconducting ground state that breaks multiple symmetries, then, as we show below, it has the ideal ingredients for the formation of fluctuationinduced electron quadrupling states.
Multiple broken symmetries imply a multicomponent order parameter. Hence, it is described by multiple complex fields of the form |∆ i |e iφi . Consider a system that is a two-dimensional multicomponent superconductor: at finite temperature, and for a finite magnetic-field penetration length, the only nonvanishing order parameter in the thermodynamic limit has to be constructed out of at least four fermionic fields [6][7][8]. This is based on the observation that composite superconducting vortices, which have identical phase winding in all components, have finite energy due to supercurrents screening effects. Therefore, a fluctuating two-component system is unstable to the proliferation of composite vortices that disorder the superconducting phase, while preserving the relative density or the phase difference between the components of the order parameter. The phase difference φ i − φ j ∝ arccos Re∆ i ∆ * j is an order parameter proportional to the product of two complex fields and hence represents four-fermion correla-tions. Various other realizations of four-fermion order were discussed in two-dimensional systems that exhibit multi-component superconductivity at zero temperature [6,[9][10][11][12][13][14][15].
The recent microscopic study [16] derived an effective mean-field description for twisted bilayer graphene (TBG) near half filling of the valence band (n = −2). Upon particle doping, six Van Hove singularities give rise to as many Fermi patches, which are the leading contribution to the density of states. There are two interaction types that are permitted by symmetry, and which contribute to pairing-intrapatch and inter-patch coupling. In this scenario, the resulting mean-field theory features two complex order parameters ∆ 1 = |∆ 1 |e iφ1 and ∆ 2 = |∆ 2 |e iφ2 , and a free-energy potential of the form: where α 1 ∝ (T − T c0 ), with T c0 being the mean field critical temperature, β 1 > 0 and β 1 + β 2 > 0 for stability. The free-energy potential Eq.(1) permits two different ground state manifolds, that are determined by the sign of the coupling β 2 . For β 2 > 0, the ground state is a chiral superconductor that breaks time-reversal symmetry, while for β 2 < 0 the superconducting state develops a nematic order. Finally, for β 2 = 0, the potential exhibits an SU (2) symmetry. Note that the potential terms of the model are similar to two-component models that appear in many other instances of superconductors that break time-reversal symmetry [17], hence our results apply to other models as well.
According to the Mermin-Wagner theorem [18], two-dimensional systems with short-range interactions cannot spontaneously break a continuous symmetry at finite temperatures. However, while for the SU (2) symmetry case (β 2 = 0), the system does not exhibit any phase transition, the presence of a biquadratic term β 2 > 0, that explicitly breaks the arXiv:2206.02698v2 [cond-mat.str-el] 9 Mar 2023 SU (2) symmetry into a U (1) × Z 2 symmetry, allows for the emergence at low temperatures of an algebraic-ordered superconducting state that additionally breaks a Z 2 symmetry. When the magnetic field screening is negligible, a two-dimensional system preserves a U (1) symmetry at any finite temperature [18], while it exhibits a SC phase transition belonging to the Berezinskii-Kosterlitz-Thouless universality class [19][20][21].
In the limit of strong symmetry breaking, one may consider the London limit and the model can then be mapped onto effective models considered in [11,[22][23][24]. However, for small K, i.e. in the vicinity of the SU (2) point, fluctuation effects in the density sector may impact the resulting phase diagram.
In this letter, we focus on the effect of fluctuations in the microscopic model [16] in the scenario β 2 > 0. Starting with the free-energy functional proposed in [16], we employ large-scale Monte Carlo simulations to obtain the phase diagram of the system beyond the mean-field approximation.

THE MODEL
The Ginzburg-Landau free-energy density of the system reads: When coupled to a gauge field, this model only exhibits a quartic order in the thermodynamic limit [6]. However, for the case of TBG, the screening is negligible. Hence, we consider the problem of computing the phase diagram in the extreme type-II limit. The resulting description is characterized by an SU (2) symmetry that is explicitly broken down to U (1) × Z 2 by the β 2 term. The symmetry-breaking term renders fluctuations of the relative density massive. These fluctuations can be important in this model for the statistical problem of assessing the SC and the Z 2 critical temperatures. Correspondingly, we retain them as part of our description, while taking the total density to be constant, Rescaling the free energy by the total density |∆ 0 | 2 = |α 1 |/(β 1 + β 2 ), one can express Eq.(2) as a function of a single parameter K: given by K = β2 β1+β2 > 0. Next, by collecting the phase-difference gradient terms we obtain the free en-ergy: with For finite values of K, at low but finite temperatures, the system exhibits an algebraic-ordered SC state, that is destroyed at higher temperatures. However, a spontaneous symmetry breaking does occur in the Z 2 sector, which is associated with the two-fold degeneracy of the phase difference φ 1,2 = φ 1 − φ 2 = ±π/2, resulting from the presence of the biquadratic Josephson term.
To obtain the phase diagram of the model (4), and in particular identify the presence of a fermionic quadrupling condensate, it is necessary to assess, as a function of the parameter K, the two critical temperatures T BKT and T Z2 c . For (i) T BKT > T Z2 c , there arises a superconducting phase that preserves time-reversal symmetry; while for (ii) T Z2 c > T BKT , a metallic state that breaks the time-reversal symmetry forms as a result of the condensation of fermion quadruplets [11,22,25,26]. The observation of a quadrupling-fermionic condensate was recently reported in the three-dimensional material Ba 1−x K x Fe 2 As 2 [25].
The problem of whether a multicomponent system has a single transition or a four-fermion order is very complicated to assess, and most of the progress on such systems to date comes from large-scale numerical simulations [11,22,27,28]. Indeed, these nonsuperconducting phases are large and directly amenable for analytical arguments only in a few cases, such as two-dimensional superconductors with finite magnetic field penetration length [29] or systems where such order can be induced and tuned by an external magnetic field [7].
The reason why the problem is that challenging is that beyond the mean-field approximation, the physics of the system, and therefore its phase diagram, is governed by the proliferation of topological phase excitations that mutually interact with each other. These can be elementary vortex excitations, resulting from a phase winding in each condensate individually; composite vortices, resulting from the phase winding of both condensates around the same core; and domain walls separating regions with opposite phase differences. The elementary vortices (∆φ 1 = ±2π, ∆φ 2 = 0) ≡ (±1, 0) or (∆φ 1 = 0, ∆φ 2 = ±2π) ≡ (0, ±1) have a phase winding in the intercomponent phase difference and hence emit a domain wall. Consequently, their proliferation restores the Z 2 symmetry and simultaneously destroys the superconducting state leading to the BKT superfluid-stiffness jump to zero at the critical point. On the other hand, the proliferation of composite vortices of the kind ±(1, 1) can only affect the superconducting sector, leaving the Z 2 symmetry broken. Likewise, the proliferation of domain-wall excitations alone can only restore the Z 2 symmetry leaving the superfluid stiffness associated with the SC phase finite. The key problem is that the defects in the U (1) and Z 2 sectors are not decoupled. First, in contrast to the ordinary vortices [30], the composite vortices in this model consist of two spatially separated fractional vortices because of the condition ρ 2 = |∆ 1 | 2 +|∆ 2 | 2 = 1. Such defects carry a skyrmionic topological charge and exist also when one softens the ρ 2 = 1 constrain [31]. That implies that in contrast to conventional vortices thermal excitations in the U (1) sector also generate local defects in the phase-difference sector. Likewise, a thermally induced Ising domain wall interacts with vortices by splitting them into two half-quanta vortices [31], hence disorder in the Z 2 sector may induce disorder in the U (1) sector as well. As a consequence, it would in principle be incorrect to assess the critical temperatures of these two sectors by treating them separately, there are several correlation lengths and their interplay is highly nontrivial.
The BKT superconducting transition is associated with the emergence of a finite stiffness of the phasesum. Within the Ginzburg-Landau model Eq.(3), this can be assessed by computing the helicity-modulus sum Υ µ + , defined as the linear response of the system to an infinitesimal twist of the two phase condensates along the direction µ: where: Here, δ µ,i denotes the phase-twist parameter with respect to the i-th phase component. Here, L is the linear size of the two-dimensional system. The expectation value . . . is the thermal average, evaluated stochastically by the Monte-Carlo Metropolis algorithm. In our simulations, we compute the helicitymodulus sum along µ =x. In what follows, we will simply write: Υ + ≡ Υx + .
Ordinary U (1) systems in two dimensions exhibit a topological phase transition driven by the unbinding of vortex-antivortex pairs [20,21,33], which becomes entropically favorable at a finite temperature T BKT . The proliferation of free vortices leads to a discontinuous vanishing of the phase stiffness, that drops to zero at T BKT according to the Kosterlitz-Nelson universal relation [30].
When a system undergoes a BKT phase transition [20,21,33], the critical point can be located by finitesize scaling of the quantity [34]: where L 0 is a free parameter giving the best crossing point at finite temperature (see also Supplementary information [32] and Supplementary Fig. 1). For K = 5, the best crossing point is obtained for L 0 = 3, as shown in Fig. 1(a). Varying K, the value of L 0 varies as well. In particular, we find that L 0 increases with decreasing K (see Supplementary Fig. 2 and Supplementary Fig. 3), leading to very pronounced finite-size effects at small K. This finding stems from the multi-component nature of the system. Indeed, in contrast to the single-component case, the BKT transition is in this case driven by the proliferation of free composite vortices, resulting from the unbinding of a pair formed by a (1, 1) and a (−1, −1) vortex.
For large values of K, the superconducting phases of the two condensates are essentially locked, and the model (3) can effectively be described by a single component. In this limit, the two elementary vortices ±(1, 0) and ±(0, 1) that constitute ±(1, 1) composite vortices are tightly bound. However, for smaller values of K, this is no longer the case. Indeed, alongside the density-density interaction that promotes the separation of the composite vortices into their elementary constituents, in the limit K → 0 the model approaches the SU (2) symmetry where the composite vortices are unstable in a conventional sense. The finite size of these composite vortices results in an increase of the finite-size effects of the whole system leading to a larger value of L 0 . That also suggests that standard conventional-vortex-based estimate for the BKT transition are not accurate in this limit. To asses the Z 2 phase transition, we define an effective Ising order parameter m, related to the two possible values of φ 1,2 ∈ [−π; π) via: Finally, we extract the Z 2 critical temperature by means of a finite-size crossing analysis of the Binder cumulant U of m: In the thermodynamic limit, U tends to 1 in the hightemperature phase and to 1/3 in the low-temperature limit. The resulting crossing point for the case K = 5 is shown in Fig. 1 (b). The phase diagram obtained via this numerical study is shown in Fig. 2. Details on the finite-size scaling of the two critical temperatures can be found  Fig. 4-6). Our results reveal that for any finite value of K we considered, the system has a fermionic quadrupling state that breaks time-reversal symmetry. The range of temperatures where this phase appears (see the inset of Fig. 2 ) is larger for large values of K and saturates to a finite value in the limit K → ∞. The presence of a lattice provides a minimum size for the domain wall between two different chiralities. Consequently, the energy cost of such topological defects saturates to a finite value in the limit K → ∞, resulting in a saturation of the critical temperature associated with the Z 2 transition. In contrast to the previously studied conventional multiband models, in the limit of a very small intercomponent coupling K, the two transitions do not merge and we observe a relative increase of ∆T c for K < 1. We argue that this increase is related to the symmetry of the model in the limit K → 0 that, for the model derived in [16], is SU (2), rather than U (1) × U (1) symmetry as for the case of s + is superconductors. In two dimensions, SU (2)-symmetric systems exhibit no long-range or quasi-long-range order. Obtaining a significant fermion quadrupling phase in the case when U (1) × U (1) symmetry is explicitly broken to U (1) × Z 2 generally requires a very strong symmetry-breaking Josephson term [11]. By contrast, in the TBG model considered in this work, a notice-able quadrupling phase remains even if the term that breaks the Z 2 symmetry is small. In all the considered cases however the positions of the critical points are very correlated, and the difference in critical temperatures is of the order of 1% signalling that fluctuations in the U (1) and Z 2 sectors are nontrivially coupled. Therefore, the phase diagram of this or similar models cannot be accurately determined by treating the U (1) and Z 2 sectors separately.

CONCLUSIONS
In conclusion, we studied a Ginzburg-Landau model derived in connection with twisted bilayer graphene [16] at low temperatures. We have shown that this class of models can host a fermion quadrupling phase above the critical temperature of the superconducting phase. This phase is more robust than its counterpart in the previously studied class of models [11,[22][23][24] and extends to all finite values of the coupling parameter K considered but is still relatively small. Nonetheless these findings indicate that magic-angle twisted bilayer graphene can be an especially promising platform for realizing and observing fermion quadrupling order.
While the temperature range associated with fermion quadrupling is small, it is likely to be larger in real systems. The coupling to a vector potentialomitted in this work-reduces the energy of composite vortices, thus reducing the temperature of the onset of superconductivity so that a finite diamagnetism could lead to a significant fermion quadrupling phase. More importantly, the robustness of this phase in the model we studied suggests that its size can be amplified by applying a transverse magnetic field [7,25].
The fermion quadrupling state can be identified via a combination of thermal and electrical transport measurements, analogous to those performed in [25]. The effective model of the Z 2 quadrupling state [35] suggests that signatures of a Z 2 broken symmetry above the critical temperature can be detected via magnetic probes. Skyrmion excitations [35] or spontaneous magnetic fields can indeed appear in the presence of local strain, obtained by imposing local pressure or by local heating [36,37] in combination with local magnetic probes. Another route to probe this state in twisted bilayer graphene is through collective modes [38][39][40][41]. The inter-component collective modes, indeed, only depend on the relative phases and relative densities of the two components and thus they should survive in the non-superconducting state with broken time-reversal symmetry.