Effective model for superconductivity in magic-angle graphene

We carry out large-scale quantum Monte Carlo simulations of a candidate field theory for the onset of superconductivity in magic-angle twisted bilayer graphene. The correlated insulating state at charge neutrality spontaneously breaks U(1) Moir\'e valley symmetry. Owing to the topological nature of the bands, skyrmion defects of the order parameter carry charge $2e$ and condense upon doping. In our calculations we encode the U(1) symmetry by an internal degree of freedom such that it is not broken upon lattice regularization. Furthermore, the skyrmion carries the same charge. The nature of the doping-induced phase transitions depends on the strength of the easy-plane anisotropy that reduces the SU(2) valley symmetry to U(1) $\times \mathbb{Z}_2 $. For large anisotropy, we observe two distinct transitions separated by phase coexistence. While the insulator to superconducting transition is of mean-field character, the U(1) transition is consistent with three-dimensional XY criticality. Hence, the coupling between the gapless charge excitations of the superconducting phase and the XY order parameter is irrelevant. At small anisotropy, we observe a first-order transition characterized by phase separation.

Introduction.-Magic-angletwisted bilayer graphene (MATBG) provides a new platform to study correlationinduced phenomena.Aside from correlated insulating states at commensurate fillings, it is of great present interest to understand the physics when the electron filling factor is close to charge neutrality [1].In particular, extensive attention has been paid in understanding how the correlated insulator gives way to a superconducting state via doping.Aside from transport properties [2,3], there is a lack of experimental tools [4] to study these two phases.A theoretical hypothesis from Ref. [5,6] attempts to explain the two states within a unifying framework.In this framework the SU(2) valley symmetry is weakly broken down to U(1) ×Z 2 .The so-called Kramer inter-valley coherent insulator (K-IVC) [7] spontaneously breaks the U(1) charge conservation within each Moiré valley; skyrmion defects of the order parameter carry electron charge 2e, which, upon doping condense and trigger superconductivity [5,6,8].The success of this continuum limit picture relies crucially on the conservation of valley quantum numbers corresponding to the chiral symmetry.
In constructing an effective lattice Hamiltonian for the above, a major difficulty arises since the valley (or chiral) symmetry will invariably be broken by the regularization.On the lattice, Dirac cones cannot be rotated independently.Furthermore, within the theory proposed For each value of ∆ we consider a value of λ that maximizes the K-IVC order at charge neutrality (µ = 0).Red bullets denote the superconducting transition points while the blue dot corresponds to the K-IVC transition point.The box with dashed line indicates the uncertainty in identifying the position of critical end point, which depends on the specific choice of λ at each value of ∆.The inset illustrates the spin-orbit interactions inside a plaquette.
in Ref. [5] the spin does not play a key role in the pairing mechanism, and in fact spinless versions of the theory were put forward in Ref. [8].In this letter we propose to swap chiral and spin symmetries.We use the spin degree of freedom to encode the (Moiré) valley degrees of freedom and the two Moiré valley bands are reformulated arXiv:2304.02428v2[cond-mat.str-el]3 May 2023 in our case as physical spin (up and down).It follows naturally [9] that, the U(1)×Z 2 symmetry-broken states (K-IVC and Valley Hall, [7]) map onto a dynamically generated quantum spin Hall (QSH) insulator with easyplane anisotropy.From previous works, Refs.[10][11][12], it is known that skyrmion excitations of the QSH order parameter carry charge 2e.In fact in the absence of easyplane anisotropy, it was argued in Ref. [12] that the proliferation of skyrmions simultaneously destroys the QSH order and generates superconductivity (SC) as the system is doped.Hence, the low energy physics of MATBG can be realized by including easy-plane anisotropy, such that merons, as opposed to skyrmions, become the low energy charged textures.
Although there is no essential difference between (pairs of) merons and skyrmions in terms of the defining topology and their associated electron charges, the symmetry difference between U(1) ×Z 2 and SU(2) can lead to dramatic effects.Depending upon the easy axis anisotropy, meron (pairs) may have a large excitation gap and eventually result in a doping-induced phase transition of mean-field character.This letter aims at revealing the Cooper-pair condensation when doping a U(1) broken symmetry K-IVC insulator and at understanding the interplay between the two order parameter fluctuations.
Model and Method.-Weconsider a model of Dirac fermions in 2 + 1 dimensions on the honeycomb lattice with Hamiltonian Ĥ = Ĥt + Ĥλ .Here, The spinor ĉ † i = ĉ † i,+ , ĉ † i,− where ĉ † i,τ creates an electron at lattice site i with z-component of the internal degree of freedom τ , and the sum runs over the nearest neighbors of the honeycomb lattice.
The SU(2) invariant version (∆ = 1) of this Hamiltonian has been well studied in Ref. 11 and 12. Associating τ to the spin degree of freedom, a dynamically generated QSH insulator that spontaneously breaks SU(2) τrotational symmetry is found at intermediate interacting strength (λ) and at half-filling.An SSC can be realized by increasing λ or doping.For ∆ < 1, we realized a U(1) symmetry-broken QSH state by reducing the full τ -rotational symmetry of the Hamiltonian [14].
The Hamiltonian we considered here captures the key ingredients of MATBG around charge neutrality.In MATBG, the Chern number of the flat bands is associated with the (Moiré) valley quantum number [7].As documented in the supplemental material, the spin degrees of freedom in our toy model play the role of the (Moiré) valley in their case.Therefore, the Kramers inter-valley coherent phase, which breaks spontaneously the U(1) valley symmetry is not different than our QSH insulator from a symmetry point of view.The following consequence from this symmetry argument is that in both cases, electron pairing arises from Kramers doublet based on pairs of meron configurations in the U(1) broken insulator.Upon doping, these pre-formed pairs condensate and superconductivity is formed.Crucially, although charge conservation within each valley (e.g., chiral symmetry) is an exact symmetry in the continuum limit, it is generically not possible realize it in a lattice Hamiltonian.A simple way of performing this regularization is to substitute the valley degrees of freedom with spin ones [15].
Here we focus on the case where ∆ ∈ [0, 1) such that the SU(2) spin rotational symmetry is reduced to U(1)×Z 2 .We investigate three values of the anisotropy and only consider values of λ deep inside the chargeneutral K-IVC state [14]: λ = 0.043 for ∆ = 0.1, λ = 0.035 for ∆ = 0.5, and λ = 0.03 for ∆ = 0.75.Our unit of energy is set by t = 1.We used the ALF (Algorithms for Lattice Fermions) implementation [16] of the auxiliary-field quantum Monte Carlo (QMC) method [17][18][19].Because λ > 0 and ∆ > 0, a real Hubbard-Stratonovich decomposition for the perfect square term does not break the time-reversal symmetry.Since charge is conserved, the eigenvalues of the fermion determinant come in complex conjugate pairs such that no sign problem occurs [20].We simulate periodic systems with size of L × L. The imaginary time interval is ∆τ = 0.2 and a symmetric Trotter decomposition is chosen to ensure the hermiticity of the Trotterized imaginary time evolution [11].Additionally, we apply a checkerboard decomposition to the exponential of hopping matrix Ĥt .Following our previous work [12], we used a projective QMC algorithm (PQMC) [19,21,22].For the trial wave function, we take the ground state of the hopping Hamiltonian in Eq. (1) with spin-dependent twisted boundary conditions.This Slater determinant is then uniquely defined and also time-reversal symmetric.
QMC results.-Tocapture the K-IVC order, we define the local operators Ôr,n ≡ Ĵr+δn,r+ηn .Here, r denotes a unit cell and n = 1, 2, ...6 are the six next-nearest neighbor bonds of the corresponding hexagon with legs r + δ n and r+η n .The K-IVC order with broken U(1) rotational symmetry is detected by computing: As for SSC we consider: where S O is the largest eigenvalue of the corresponding correlation matrix (O=K-IVC,SSC).q 0 is the ordering wave vector and q 0 + δq is the neighboring wave vector (|δq| = 4π √ 3L ).We define the doping factor δ by the density of doped holes relative to charge neutrality, δ ≡ 2L FIG. 3. Same as Fig. 2, for ∆ = 0.5.
N e counts the number of electrons.Due to large autocorrelation times related to particle fluctuations, it is convenient to adopt a canonical ensemble.Within this ensemble, the chemical potential is evaluated from the slope of the 'tower of states' in the charge sector: Here, ∆ η − (N e ) is the pairing (η − ) gap extrapolated from the time-displaced correlation function.
Our results are summarized in the ground state phase diagram in the ∆ versus µ plane, as shown in Fig. 1.For each value of ∆ we consider a coupling λ that places us well within the correlated insulating phase at charge neutrality [14].Below we will document that at small anisotropy we observe a first-order transition characterized by phase separation between the K-IVC and SSC.At larger values of ∆ phase coexistence sets in.
Figures 2 and 3 show the correlation ratios for the SSC and K-IVC instabilities at ∆ = 0.1 and ∆ = 0.5.As apparent, the K-IVC order survives finite doping, with the critical doping, δ K-IVC c , being given by the crossing point of the curves.On the other hand, superconductivity sets in at any finite doping for both considered values of the anisotropy, ∆.Since we are working in a canonical ensemble, we have to check for the stability against phase separation as signaled by an infinite compressibility, ∂δ ∂µ .
Figure 4(b) plots this quantity at, ∆ = 0.5.As apparent, in the range 0 < δ < δ K-IVC c the data is consistent with a diverging compressibility.Hence, in this doping range and in the thermodynamic limit, we expect to observe puddles of K-IVC insulating phases with a total density set by 1 − δ δ K-IVC c and regions of SSC with a total density set by δ δ K-IVC c [23].On the other hand, at strong anisotropy, ∆ = 0.1, Fig. 4(a) shows no sign of diverging compressibility, thus indicating phase coexistence.However we observe two non-analytical points, at µ c1 ≈ 0.31 and µ c2 ≈ 0.35 corresponds to the insulator to superconductor transition at δ = 0 + and to the K-IVC transition at δ = δ K-IVC c .Both transitions show no features of first-order transitions, as can be directly seen from the continuous behavior of δ −µ dependence in Fig. 4(a).We note that at the superconducting transition at µ c1 , the hyper-scaling law: holds.This is due to the fact that the generator of the SSC order parameter couples to the chemical potential, µ, corresponding to the tuning parameter of this quantum phase transition [24].The linear δ-µ dependence at µ c1 ≈ 0.31 in Fig. 4(a) suggests that z = 2 at the superconducting phase transition.Assuming that the background of K-IVC ordering does not couple to any critical fluctuations at the superconducting critical point, a mean field phase transition with z = 2, ν = 0.5 at µ c1 is expected.This transition is in the very same universality class as that of the Mott-insulator-superfluid transition in doped Bose-Hubbard system [24].
To understand the nature of the K-IVC transition in the background of superconducting order, we fit our data to the following form: (q 0 )/L 2 is the K-IVC order parameter.Our results, Fig. 5, suggest that the transition is consistent with Lorentz invariance, z = 1, and that it falls into the 3D XY universality class.In particular we considered ν ≈ 0.67169 and η ≈ 0.03810 from previous Monte Carlo simulations of the 3D XY model [25] as well as δ K-IVC c = 0.0448 for our data collapse.As shown in Fig. 5, both two quantities show nice collapse for system sizes from L = 6 to L = 18.We also performed a collective polynomial fit using L = 15 and 18 based on Eq. ( 8).The results are listed in Tab.I.Even when a wide fitting range is taken into consideration, the value of ν that we obtain is consistent with 3D XY universality class.
An interesting question to ask is why symmetry allowed coupling terms between the K-IVC order parameter and the critical charge fluctuations of the superconductor play no role at the 3D XY phase transition.In Ref.
[9] we provide a power counting argument to show that the Goldstone modes of the SSC order parameter are irrelevant at the 3D XY fixed point.
Discussion and summary.-Wehave considered an effective lattice model that captures the physics of a candidate theory of MATBG [5], which unifies the K-IVC insulator and superconducting phases: skyrmions defect of the K-IVC order parameter carry charge 2e and condense upon doping.The key insight to obtaining lattice regularization of this physics is to encode the valley symmetry with the spin degree of freedom.In Ref. [9] we show that the single-particle gap does not vanish.This allows us to integrate out the electronic degrees of freedom and precisely obtain the low energy topological field theory presented in Ref. [5].Large-scale QMC simula- FIG. 5. Data collapse at ∆ = 0.1 with δ K−IVC c = 0.0448 and 3D XY exponent whose ν = 0.67169 (7), η = 0.03810( 8) for (a) the correlation ratio and (b) the K-IVC order parameter.
tions reveal two different types of transitions depending upon the strength of the easy-plane anisotropy.At large anisotropy we observe phase coexistence.As a function of doping the insulating state gives way to a superconducting phase, with the universality class being identical to that of the Bose-Hubbard model [24].At larger doping we observe the vanishing of the U(1) order in the background of the superconducting phase.Our results show that this phase transition belongs to the 3D XY universality class, such that coupling to gapless charge fluctuations are irrelevant at this critical point.
For small anisotropic case, we observe phase separation.In conjunction with our previous results for the SU(2) case [12], we observe that the doping range in which phase separation occurs decreases as one approaches the SU(2) symmetric point, ∆ → 1.This challenges the conclusion of Ref. [12] where a seemingly continuous transition with large dynamical exponent is observed.Alternatively, if symmetry plays a crucial role, an interesting possibility is that the large z continuous transition only exists in the SU(2) symmetric case.Our parameters are such that the single-particle gap at charge neutrality is independent of the anisotropy and the pairing gap decreases as the anisotropy grows [9].We inter-pret this in terms of merons pairs that become energetically more expensive due to a smaller core size at strong anisotropy.It is hence tempting to interpret the phase diagram of Fig. 1 as dominated by topology in the vicinity of the SU(2) symmetric point and from the point of view of a Ginzburg-Landau theory at strong anisotropy.This statement is substantiated by calculations in Ref.
[9] at small anisotropy, that show the locking of the charge density and curvature of the K-IVC order parameter.
Let us now return to MATBG, where the easy-plane anisotropy is expected to be small [7].The numerical study in Ref. [8] suggests that the mechanism of skyrmion superconductivity is stable to the Coulomb repulsion, such that our model can very well capture the low energy physics of MATBG.In this case, our results suggest that the doping-induced transition is first-order and entails phase separation.Here we present data for ∆ = 0.75.Upon inspection of Fig. 6(a)-(c) the data bears very similarities to the case of ∆ = 0.5 presented in the main text.In particular, we observe a first-order phase transition between the K-IVC and SSC.However in comparison to the case ∆ = 0.5 we see that the doping range where phase separation occurs 0 < δ < δ K-IVC c is smaller.In this section, we show the single-particle spectrum at finite chemical potential.At zero temperature we can sharply distinguish between the particle addition and removal spectra:

Single-particle spectrum
We obtain the spectral function by analytic continuation in particle and hole channels: Here, |0 in Eq. ( 10) is the ground state at finite doping and n| is an eigenstate of the Hamiltonian with energy E n and an additional particle (hole) relative to the ground state.We have used the stochastic Maxent [26] implementation of the ALF [16] library.
In Fig. 7, we plot the spectral functions for the L = 18 lattice and dopings, δ = 0, 6 648 , 14 648 , and 28 648 .As apparent, we always observe a single-particle gap in the spectra and an approximate particle-hole symmetry around the Fermi level.At finite doping, this stems from the superconducting nature of the ground state and the energy cost of breaking a Cooper pair.Importantly, since fermion excitations are gapped, one can integrate them out, to derive the purely bosonic field theory discussed in Ref. [5].

Preformed pairs
Here, we would like to clarify two points: i) pairs of merons are the lowest energy excitations and ii) the gap of preformed pairs increases as a function of anisotropy.
We extrapolate the excitation energy of the pairing gap ∆ SC as obtained from the SSC imaginary-time correlations, as well as the fermionic single-particle gap from the Green's function: We note that the minimal gap is at the momentum M point in the Brillouin zone (see Fig. 7).
Overall, the single-particle gap ∆ sp is larger than half of the pairing gap ∆ SC /2 for all three values of ∆, thus indicating that the lowest charge excitation is a pair.The pairing energy is defined as ∆ Pairing = 2∆ sp − ∆ SC .Let us now consider only the even lattice sizes for which the single-particle gap shows little size effects.In this case, we observe that the pairing energy decreases as the easyplane anisotropy grows.This is consistent with our understanding based on the meron picture: although the norm of the K-IVC U(1) order parameter (corresponding to the single-particle gap) has no significant dependence on the anisotropy, pairs of merons require less excitation energy as the anisotropy grows.This stems from our understanding that as the anisotropy grows the meron core becomes more energetically expensive.

Spin current texture
The fact that skyrmions in S 2 space carry charge 2e can be understood locally: electron densities are associated with the curvature of the local spin-orbit fluctuations.Correspondingly, spin current textures will also be observed when the electron (hole) distribution is not uniform.In this section, we will demonstrate numerically the one-to-one relation between these two quantities.
We trap a pair of holes at two separate regions of the honeycomb lattice, by adding a site-dependent pinning potential: such that the electron density relative to half-filling is reduced around the center of two potential wells r c1 and r c2 , as shown by Fig. 10(a).
Correspondingly, the curvature of the spin-orbit order parameter is given by where r denotes the position of unit cell and m labels the six next nearest neighbor bonds.We use m = 1, 2, 3 respectively on hexagons at site r, r + a 1 and r + a 2 , as shown in Fig. 9.Here A is the area of a hexagon.This way of defining curvature is to avoid explicit density operators via the commutator of spin current operators.From Fig. 10(b), it is clear that the curvature distribution in real space follows exactly the spatial pattern of the electron density.

Scaling analysis of the coupling between two order parameters
The aim of this section is to show that at the 3D XY critical point coupling to charge fluctuations is irrelevant.To this end, we write down the ϕ 4 theory describing the K-IVC phase transition within the superconducting state: where S K-IVC is the bare action for K-IVC field, and S SSC is the one for SSC order.Here we denote the two U(1) order parameters for K-IVC and SSC as n and φ, respectively.In the SSC phase, we can decompose the field as  S coupling describes the charge-or spin-neutral coupling terms between two fields.The leading term reads where D = d + 1 and d = 2.Note that terms like d D r|n| 2 |φ| 2 will not be present at this critical point, since φ has a finite mass along the longitudinal direction in the SSC phase.
Let ∆ n be the scaling dimension of the field n and ∆ θ be the scaling dimension of the Goldstone boson θ.According to standard scaling theory, the scaling dimension of λ should be ∆ Assume that in our system, the Wilson-Fisher fixed point of the field n belongs to the 2+1D universality class ( ∆ n = D−2+η 2 = 0.51905, with η = 0.03810).On the other hand, the Goldstone mode θ of a superconducting state in 2 + 1D is associated with a power-law decay correlation function along the transverse direction: This is based on the saddle point expansion of ϕ 4 theory along the massless direction.Therefore, ∆ θ = 0.5, and as a consequence, ∆ λ = D−2∆ n −2∆ θ −2 < 0. This implies that the coupling between two U(1) order parameters is irrelevant.This is consistent with our numerical evidence of the 2+1D O(2) universality class at the K-IVC critical point.
Relation to magic-angle twisted bilayer graphene (MATBG) Our model is closely analogous to those described in Refs.[5,7] that aim at accounting for superconductivity resulting from the condensation of skyrmions in magicangle twisted bilayer graphene (MATBG).The starting point is the Bistritzer-MacDonald continuum model [27] in the chiral limit, in which interlayer hopping between the same sub-lattice is set to zero.In this case, the lowlying bands are flat and can be labeled by sub-lattice polarization, which is complete in this limit, as well as a valley and spin index.In this basis, the bands carry non-trivial topology characterized by a Chern number.Let ĉ † τ,σ,s (k) create an electron with momentum k + K τ , in valley τ , sub-lattice polarization σ, and physical spin s.The Chern number of the bands is given by σ z τ z .In this chiral limit, the form factor of the density fluctuations depends solely on σ z τ z , such that the complete model possesses a U(4) × U(4) symmetry that rotates the bands in a given Chern sector.This model captures the dominant energy scales.Perturbations beyond the chiral limit break this symmetry.In fact, the phases and phase transitions discussed in [5,7] are captured by the effective model: Ĥ= + ĉ † (x)τ y σ y ĉ (x) 2 + ∆ ĉ † (x)σ z ĉ (x) 2 P .
Here, P reflects the projection onto the low energy Hilbert space.
At ∆ = 1 the model has an SU V (2) valley symmetry with generators i 2 M VH i , M VH j , corresponding to (τ x σ x , τ y σ x , τ z ), as well as an SU S (2) spin symmetry and a U C (1) charge symmetry.For ∆ = 1, the SU V (2) symmetry is reduced to a U V (1) × Z 2 and for ∆ < 1 the K-IVC state is favored.Since the bands have a nontrivial Chern index, skyrmion excitations of the threecomponent K-IVC and VH order parameters carry charge 2e.We note that, since the SU V (2) is reduced to U V (1) × Z 2 , skyrmions have to be seen in terms of a pair of merons, each carrying charge e.Hence, all in all, the model has U V (1) × Z 2 × U C (1) × SU s (2) symmetry.The doping-induced transition between the K-IVC and SSC put forward in [5] and based on the proliferation of skyrmions, does not involve the spin degrees of freedom.Spinless versions of the model have been put forward to capture the relevant physics [5,8].
The model of Eq. 18 does not support lattice regularization since it will break the valley (or chiral) symmetry.Furthermore, the topology of the bands does not allow for a local Wannier basis.As a result, simulations of MATBG are carried out in the continuum [8,28,29].Our model provides a possibility to avoid this by encoding the SU V (2) symmetry as SU S (2).A continuum limit of our model reads: −λ V d 2 x ĉ † (x)s x τ z σ z ĉ (x) 2 + ĉ † (x)s y τ z σ z ĉ (x) 2 + ∆ ĉ † (x)s z τ z σ z ĉ (x) 2 , where M QSH = (s x , s y , s z )τ z σ z correspond to the three QSH mass terms.The parallel now becomes apparent: the M QSH z mass corresponds to the Valley Hall insulator, and the first two components to the K-IVC insulator.Topologically, both models are equivalent since the skyrmion of the three-component QSH order parameter carries charge 2e.The symmetry of the Hamiltonian is given by U s (1) × Z 2 × U C (1) × SU V (2).Here the SU V (2) symmetry is again spanned by the generators: (τ x σ x , τ y σ x , τ z ).

FIG. 1 .
FIG.1.Schematic phase diagram in the anisotropy ∆ and chemical potential µ plane.Blue line indicates first-order phase transition, while black lines continuous ones.For each value of ∆ we consider a value of λ that maximizes the K-IVC order at charge neutrality (µ = 0).Red bullets denote the superconducting transition points while the blue dot corresponds to the K-IVC transition point.The box with dashed line indicates the uncertainty in identifying the position of critical end point, which depends on the specific choice of λ at each value of ∆.The inset illustrates the spin-orbit interactions inside a plaquette.

FIG. 2 .
FIG. 2. Equal-time correlation ratio for (a) R K−IVC and (b) R SSC as a function of doping factor δ for ∆ = 0.1.The vertical dashed line is a guide to the eye, fitted from the crossing point of K-IVC order parameter between L = 15 and 18.

5 FIG. 4 .
FIG. 4. Doping factor δ as a function of chemical potential µ for (a) ∆ = 0.1 and (b) ∆ = 0.5.The horizontal line is a guide to the eye, which is the same as the vertical one in Figs. 2 and 3.

FIG. 6 .
FIG. 6.Data for ∆ = 0.75.Correlation ratios for (a) R K−IVC and (b) R SSC as a function of doping factor δ. (c) δ as a function of chemical potential µ.The vertical and horizontal dashed lines are a guide to the eye, fitted from the crossing point analysis of the K-IVC order parameter.Instead of choosing the largest two sizes, we fit the crossing points between lines of L and L + 3 from L = 6 to L = 15 and extrapolate to the thermodynamic limit.

FIG. 10
FIG.10.(a) Electron particle number distribution in each unit cell, and (b) the curvature in S 2 space of spin current in the presence of a finite pinning potential based on Eq. (13).The simulation is performed on a 12 × 12 lattice with δ = 2/288, ∆ = 0.75, and λ = 0.03.