Chiral $p$-wave superconductivity in twisted bilayer graphene from dynamical mean field theory

We apply cluster dynamical mean field theory with an exact-diagonalization impurity solver to a Hubbard model for magic-angle twisted bilayer graphene, built on the tight-binding model proposed by Kang and Vafek (2018), which applies to the magic angle $1.30^\circ$. We find that triplet superconductivity with $p+ip$ symmetry is stabilized by CDMFT, as well as a subdominant singlet $d+id$ state. A minimum of the order parameter exists close to quarter-filling and three-quarter filling, as observed in experiments.


Introduction
Twisted bilayer graphene (TBG) consists of two layers of graphene deposited on top of each other with a slight rotation, or twist. At commensurate twist angles, the bilayer forms a moiré pattern with a period that depends closely on the twist angle. It has been predicted that for some "magic angles", the resulting band structure has a few relatively flat bands at low energy, separated from 2 Low-energy model There have been a few proposals for an effective tight-binding Hamiltonian describing the lowenergy bands of TBG [1,[18][19][20]. We adopt in this work the model described in Ref. [1] and inspired by Ref. [19]. It is based on four Wannier orbitals per unit cell, with maximal symmetry, on an effective honeycomb lattice and is appropriate for a twist angle θ = 1.30 • .
It is customary to derive effective models for TBG directly from continuum models. In that framework a valley symmetry emerges and the model is endowed with a fragile topology. It can be shown that in a model with nontrivial topology, time-reversal symmetry (TRS) cannot be represented simply by a set of localized Wannier states: its action is not strictly local [21]. However, as shown in [22], the error committed by using a localized Wannier basis is exponentially small. Since we are going to truncate the hopping matrix to a few terms and introduce strong interactions that would likely destroy any existing topology, this issue should not be of concern here. Fig. 1 offers a schematic view of the orbitals w 1 and w 3 . Orbitals w 2 = w * 1 and w 4 = w * 3 are not shown. Ref. [1] computes a large number of hopping integrals, of which we will only retain (green) on which our model Hamiltonian is built. The charge is maximal at the AA superposition points (blue circles) forming a triangular lattice. The Wannier functions are centered on the triangular plaquettes that form a graphene-like lattice (black dots), whose unit cell is shaded in red. The underlying moiré pattern illustrated corresponds to (m, n) = (9,8), but the functions used in this work correspond to (m, n) = (26,25). The basis vectors E 1,2 of the moiré lattice are shown (they are also basis vectors of the graphene-like lattice of Wannier functions), as well as the elementary nearest-neighbor vectors a 1,2,3 .
the largest, as listed in Table 1. The notation used is that of Ref. [1].
Remarkably, the most important hopping terms are between w 1 and w 4 (and between w 2 and w 3 ), i.e., between graphene sublattices. It therefore makes sense physically to picture the system as made of two layers and to assign w 1 and w 4 to the first layer, whereas w 2 and w 3 are assigned to the second layer. The rather small t 13 [0, 0] hopping (and its equivalents) is the only term that couples the two layers. The concept of layer is useful when visualizing the model in space and when arranging local clusters of sites in CDMFT, since it is preferable to have the more important hopping terms within a cluster; it is merely a book-keeping device. The drawing next to Table 1 illustrates the range and multiplicity of the intra-layer hopping terms retained.
To this tight-binding model we will add a local interaction term U. This is a rather approximate description of the interactions in this system, but has the merit of simplicity and tractability in the context of dynamical mean field theory. A more refined description of the interactions would not only contain extended interactions (see, e.g., [23,24]) but would include terms not of the densitydensity form [25]. We will defer the study of extended interactions to future work. The values of U in our calculations range from 0.5 meV to 5 meV. Fig. 3d of Ref. [5] leads us to expect a wide range of values of U depending on twist angle, and a rather large U ∼ 20 meV at an angle of 1.30 • . However, Ref. [11] predicts a value U ∼ 5 meV for this angle and the range of U values symbol value (meV) Table 1: Hopping amplitudes used in this work. They are the most important amplitudes computed in Ref. [1].
Here ω = e 2πi/3 and the vector [a, b] following the symbol represents the bond vectors in the (E 1 , E 2 ) basis shown on Fig. 1. Note that t 23 = t * 14 and t 24 = t * 13 . On the right: schematic view of the hopping terms t 14 within a given layer (the unit cell is the blue shaded area). Lines 2, 3, and 4 of the table correspond to the red, blue and green links, respectively. Dashed and full lines are for t 14 and t 23 , respectively. predicted in Fig. 9 of Ref. [26] is largely compatible with the range we have selected.
The model is invariant under a rotation C 3 by 2π/3 about the AA site, and under a π-rotation C 2 about an axis in the plane of the bilayer (the vertical axis on Fig. 1). These transformations generate the point group D 3 and affect the Wannier orbitals as follows [1]: where ω = e 2πi/3 andω = e −2πi/3 . In other words, the orbitals w 1 and w 3 transform between themselves, and so do w 2 and w 4 . The model also has time-reversal symmetry (TRS), under which w 1 ↔ w 2 and w 3 ↔ w 4 . Possible superconducting pairings are either singlet or triplet (there is no spin orbit coupling). It is reasonable to assume that pairing will be more important between sites that also correspond to the most important hopping integrals. Let us therefore concentrate on pairing states involving nearest neighbors on a given layer, i.e., between orbitals w 1 and w 4 (or w 2 and w 3 ). Because of the strong local repulsion in our model, we ignore on-site pairing. Let us then define the pairing operators where c r,σ annihilates an electron at graphene site r of the first layer (in orbital w 1 or w 4 depending on the sublattice). The elementary vectors a i are defined on Fig. 1, but apply to the layer in the current context. Likewise, we define operators S i,r and T i,r on the second layer, in terms of orbitals w 2 and w 3 ). Under the transformations C 3 and C 2 , the six singlet (triplet) pairing operators transform amongst themselves and may be organized into irreducible representations of D 3 , as Irrep singlet pairing triplet pairing and likewise for the combinations s , d ± id , etc. for the second layer. A similar analysis could be carried out with longer-range pairing, with the same classification: This would simply add harmonics to the basic pairing functions. This organization into representations of D 3 is contingent on the importance of the inter-layer hopping t 13 , which is an order of magnitude smaller than the intra-layer hopping. If t 13 were zero, the two layers would be independent, the symmetry would be upgraded to C 6v and the classification of pairing states would be the same as in Ref. [27], Since t 13 is small, we expect that the different pairing states of Table 2 (for a given total spin) will be nearly impossible to differentiate from an energetics point of view, except for the difference between s and d ± id (or between f and p ± i p).

Cluster dynamical mean field theory
In order to probe the possible existence of superconductivity in this model, we use cluster dynamical mean-field theory (CDMFT) [28][29][30] with an exact diagonalization solver at zero temperature (or ED-CDMFT). Let us summarize this method.

General description
The infinite lattice is tiled into identical, repeated units; this defines a superlattice, and an associated reduced Brillouin zone, smaller than the original Brillouin zone. In the present study the unit cell of the superlattice (or supercell) is made of four clusters of four sites each: Two clusters tile each of the two layers (Fig. 2c). Ref. [31] explains the particulars of CDMFT when the supercell contains more than one cluster. Each cluster is coupled to a bath of uncorrelated, auxiliary orbitals, and is governed by an Anderson impurity model (AIM): where H c is the infinite-lattice Hamiltonian, but restricted to the cluster, c i annihilates an electron on orbital i of the cluster (i labels both site and spin) and a r annihilates an electron on orbital r of the bath. The bath parameters (ε r , θ ir ) are found by imposing a self-consistency condition, as explained below.
Hamiltonian (3) is solved by exact diagonalization. Without taking into account any symmetry of the Hamiltonian, the dimension of the Hilbert space for an impurity of 4 cluster sites and 6 bath sites would be d = 4 4+6 ∼ 10 6 . Because we are investigating a broken symmetry state where the number of particles is not conserved, the only Abelian symmetry that can be used is the conservation of the z-component of the spin (we cannot use point group symmetries in general). Assuming a S z = 0 state (singlet or triplet), this reduces the dimension of the Hilbert space to 184,756.
The electron Green function on the cluster, G c (ω), is needed by CDMFT. We use the band Lanczos method to compute it; for details, please see Refs [32,33]. This method provides a Lehmann representation for the Green function. This is a L c × L c matrix, L c being the number of orbitals on the cluster (including spin). It may be expressed in terms of the electron self-energy on cluster c, Σ c (ω), and the associated hybridization function Γ c (ω): where and t c is the matrix of one-body terms of H c (including the chemical potential µ).
The fundamental approximation of CDMFT is to replace the exact electron self-energy by the self-energy obtained by assembling the various cluster self-energies: where the direct sum is carried over the various clusters forming the supercell. The Green function on the infinite lattice is then approximated by wherek is a wave vector in the reduced Brillouin zone and t(k) is the noninteracting dispersion relation expressed in real space within the supercell and in reciprocal space within the reduced Brillouin zone. If L tot = c L c is the total number of orbitals in the supercell, then G(k, ω), t(k) and Σ(ω) are L tot × L tot matrices. We further define the projected Green function This is the Fourier transform of the infinite-lattice Green function (7) to a single supercell around the origin. The CDMFT self-consistency condition requires that the L c × L c diagonal blocks of G(ω) (notedḠ c (ω)) should coincide with the corresponding cluster Green functions G c (ω). This cannot be satisfied exactly with a finite number of bath orbitals, because it should hold for all frequencies and only a finite number of bath parameters are at hand. Therefore this condition is replaced by the optimization of a distance function: where the weights W (iω n ) are chosen in some appropriate way along a grid a Matsubara frequencies associated with some fictitious temperature β −1 . This is where some arbitrariness arises in the method, as will be commented below. Let us then quickly summarize the actual CDMFT algorithm: 1. A trial value of the bath parameters (ε r , θ ir ) is chosen. When looping over an external parameter, the previous converged value or an extrapolation thereof is chosen.
2. The cluster Green functions G c (ω) are computed, with the help of an impurity solver (here an exact diagonalization method).

4.
A new set of bath parameters is found by minimizing the distance function (9) with respect to the bath parameters entering G c (ω) through Eqs (4,5) for a fixed value of Σ c .

5.
We go back to step 2 until the bath parameters or the hybridization functions Γ c converge.
Once the converged solution is found, various quantities may be computed either from the impurity model ground state (averages, etc.) or from the associated lattice Green function G(k, ω).

Cluster-bath system
The cluster-bath system for the current problem is illustrated on Fig. 2. The supercell contains four 4-site clusters; one layer is illustrated on Panel (c). Note that the only hopping term included in the impurity model is t 14 [0, 0] and its equivalents, represented by red lines on Fig. 1. The other hopping terms have an effect through the self-consistent CDMFT procedure. Each cluster contains four sites and six bath orbitals and the various bath parameters are illustrated on panels (a) and (b). The four black, numbered circles are the cluster sites per se. The six red squares are the bath orbitals. Even though their positions have no meaning, they are, on this diagram, assigned a virtual position that makes them look as if they were physical sites on neighboring clusters. They are then given "nearest-neighbor" hybridizations θ 1,2 and "secondneighbor" hybridizations η 1,2 . In order to probe superconductivity, we add pairing amplitudes within the bath itself, as shown on bath orbitals, and two others p 1,2 between "second neighbor" bath orbitals. In the context of Eq. (3), these pairing amplitudes must be understood in the restricted Nambu formalism, in which a particle-hole transformation is applied to the spin-down orbitals, giving the pairing operators the looks of hopping amplitudes. Specifically, in terms of the multiplet (C ↑ , C † ↓ , A ↑ , A † ↓ ), where C σ = (c 1,σ , c 2,σ , c 3,σ , c 4,σ ) and A σ = (a 1,σ , · · · , a 6,σ ) (σ =↑, ↓), the noninteracting part of the impurity Hamiltonian takes the form Here t c is the hopping matrix restricted to the cluster, θ is a 4×4 matrix containing the parameters θ 1,2 and η 1,2 , ε is a diagonal matrix containing the bath energies 1,2 and ∆ is a 6 × 6 matrix containing the parameters d 1,2 and p 1,2 .
In total, the AIM contains 10 bath parameters, some real, some complex. The impurity Hamiltonian does not contain pairing operators on the cluster sites themselves. However, the operators defined in Eqs (1) may develop a nonzero expectation value on the impurity through the selfconsistent bath.
The hybridization pattern shown in the figure is appropriate for triplet pairing (it is directional, as indicated by the arrows) in a p + i p state (because of the phases ω and ω 2 =ω appearing in the bath pairing amplitudes as one circles around). This may be readily adapted to probing a p − i p state (by replacing ω ↔ω) or a f state (by replacing ω,ω → 1). Likewise, singlet states are probed by introducing singlet pairing between bath sites. In principle, we could leave all pairings free, at the price of tripling the number of bath parameters, but CDMFT convergence has proven problematic when this was tested. It is easier, and no less general, to separately probe the p ± i p and f states (and likewise for the singlet states).
One could also treat the bath parameters of all four clusters of the supercell as independent. In practice, this is not necessary as they are related. The two clusters belonging to the same layer have identical bath parameters by symmetry, except for the triplet pairings which must change sign between the two clusters because the second cluster is obtained from the first by a spatial inversion. According to Table 2, we expect the complex-valued bath parameters of the second layer to be the complex conjugates of those of the first layer. These constraints effectively reduce the total number of variational parameters to the equivalent of 13 real parameters.
Minimizing the distance function (9) is done by the Nelder-Mead or the conjugate-gradient method as implemented in SciPy. These methods do not guarantee a global minimum, but only a local one. Because of this, jumps in the bath parameters might occur as a function of an external (control) parameter like the chemical potential µ, and we would expect that this manifests itself as a hysteresis when cycling over µ. We have observed no such hysteresis in the present study. This being said, the CDMFT algorithm summarized above contains an iteration over impurity models that defines a very complex nonlinear system that rather complicates this simple expectation. Failure to converge often manifests itself by oscillations between two or more sets of bath parameters and experience shows that increasing the parameter set does not necessarily alleviate this problem.

Results and discussion
We have probed the different states listed in Table 2 using the above CDMFT setup. In order to reach a solution from scratch, we have used the following staged approach: (i) Owing to the small value of t 13 , a one-layer model was first studied. (ii) An external field of each of types (2) was then applied to the cluster in order to induce a nonzero average pairing forcefully. This external field was then reduced to zero in a few steps, each time starting from the previous solution. (iii) Once a nontrivial solution was found in this way at zero external field, the second-layer was added (with a complex conjugated bath system, e.g., p − i p instead of p + i p). (iv) the solution found was then scanned as a function of chemical potential within the two-layer model. The most delicate step is to find a first solution; scanning over parameters of the model (such as the chemical potential or the interaction) is easier since the solution at a given set of model parameters provides an initial trial solution for the next parameter set. Computing time varies depending on convergence rate, but is typically of the order of 10 minutes per parameter set once the scan is in motion, with code highly optimized for speed; memory needs are relatively modest at 3-4 gigabytes. 1 We found a nonzero solution for p ± i p pairing extending over a wide range of doping. Fig. 3 shows the average p + i p order parameter on a cluster of the first layer, as a function of electron density on the cluster, for a local repulsion U = 2 meV. The order parameter is the ground-state expectation value of operator (2e) restricted to the cluster within the impurity model. Several variants of the CDMFT procedure are illustrated, which we must now explain. The distance function (9) depend on a set of weights W (iω n ) and a fictitious temperature β −1 . The values of β (in meV −1 ) are indicated in the legend of Fig. 3. The grid of Matsubara frequencies then stops at some cutoff value taken to be ω c = 2 meV in this work. The curve labeled β = 50 (blue dots) is obtained by setting all weights to the same value. The other curves (with a Σ label) are obtained by setting the weights proportional to the self-energy |Σ(iω n )| (the norm of the matrix). This is justified if one considers DMFT from the point of view of the Potthoff functional [34,35]. In particular, it gives more importance to very low frequencies in an insulating state, as the self-energy then grows as ω → 0. We expect the superconducting order parameter to be minimum, if not zero, at quarter (n = 0.5) or three-quarter (n = 1.5) filling, as observed in experiments. Indeed, this commensurate filling leads to an insulating state at the magic angle 1.08 • [5] and superconductivity occurs on either side of this filling value. We see that this is not exactly the case in the data sets of Fig. 3, although using a higher β and, to a lesser extent, a self-energy modulated set of weights, greatly helps. We will stick to the value β = 150 and use a self-energy modulated set of weights in what follows. Figure 4 shows the p + i p order parameter as a function of electron density for the full range of solutions obtained, and five values of the one-site repulsion U (in meV). We note that the system is almost (but not exactly) particle-hole symmetric. Superconductivity is strongly suppressed near half-filling (CDMFT ceases to converge to a superconducting solution when |n − 1| 0. perconductivity is partially suppressed at quarter-and three-quarter filling (n = 0.5, 1.5) and this suppression increases with U. Despite a strong suppression of superconductivity at n = 0.5 and n = 1.5, a Mott state is not fully obtained there for the range of U studied. This is likely caused by our neglect of extended interactions. Note the gap in the solutions in the vicinity of n = 0.3 and n = 1.7; the solutions exist for all values of chemical potential µ around these values, but a discontinuity leads to the forbidden regions when plotted as a function of density.
We also found a weaker singlet solution with d + id symmetry, as illustrated on Fig. 5a for U = 2 meV and U = 5 meV. The singlet solution has a smaller order parameter than the triplet solution, especially in the vicinity of n = 0.5 and n = 1.5, where it is strongly suppressed and suffers from a discontinuity (we only show the hole-doped case for clarity). A possible way to discriminate between the triplet and singlet solutions is to compare the energies of each. An optimal way to estimate the energy in CDMFT is to borrow the expression of the Potthoff selfenergy functional from the variational cluster approximation [36,37], as explained in Ref. [38]. The expression of the Potthoff functional is where E 0 is the ground state energy per site of the impurity model (including the chemical potential contribution), and the functional trace Tr represents an integral over frequencies and wave vector. It is an approximation to the grand potential Ω = E −µN of the system at zero temperature, given that the CDMFT is not far from the solution to Potthoff's variational principle [36]. Figure  Fig. 5b shows the Potthoff functional of the two solutions (p + i p and d + id) at the same time as the corresponding order parameters, as a function of chemical potential µ. The grand potential of the triplet is consistently lower than that of the singlet, except for an isolated point near a discontinuity. We have also compared directly the ground state energies E 0 of the corresponding two impurity models, and the same conclusion holds: the singlet d + id solution has a higher energy, a smaller order parameter, and is thus subdominant. We were not able to resolve the different representations of D 3 , as listed on Table 2. In other words, the energy difference between the A 1 , A 2 and E representations is likely too small to have an effect on the CDMFT convergence procedure. This is due to the small value of the inter-layer hopping t 13 . It is however important to assign opposite chiralities to the two layers.
The effective model used was based on the parameters of Ref. [1], appropriate for a twist angle θ = 1.30 • . Would our conclusions change for different, small twist angles, such as the ones found in Ref. [6] (θ = 1.05 • , 1.16 • )? Maybe. But a similar CDMFT of the nearest-neighbor Hubbard model on the graphene lattice has shown triplet pairing to be dominant [27]; so did a RPA study of bi-layer silicene [39], which is likewise based on the graphene lattice.
Let us compare our conclusions with some other works having found superconductivity in effective models for twisted bilyaer graphene. Ref. [9] finds triplet superconductivity as a Kohn-Luttinger instability, but is essentially a weak-coupling analysis, contrary to ours. Ref. [13] finds triplet superconductivity near n = 0.5, but with f symmetry, using a numerical renormalization group approach expected to be valid from weak to moderate coupling. Our strong-coupling calculations could not stabilize f -wave superconductivity. Kennes et al. [8] find d +id superconductivity near n = 1 using a renormalization-group approach followed by an mean-field analysis. Zhang et al. [14] arrive at the same conclusion, using constrained path Monte Carlo, and so do Chen et al [17]. These three works do not contradict ours, since our prediction concerns mostly regions around n = 0.5 and n = 1.5, not n = 1.
A possible improvement to the present study would be to include extended interactions, for example derived from an on-site Coulomb interaction at the AA sites [23,24]. We expect that including such interactions would hinder pairing at quarter filling. This would require adding inter-orbital interactions U 1,2 (U 3,4 ) between orbitals w 1 and w 2 (w 3 and w 4 ). Unfortunately, since orbitals w 1 and w 2 belong to different clusters in our CDMFT setup, this cannot be implemented as is. The effect could be studied within a different quantum cluster approach, such as the variational cluster approximation [27,34,40], which in practice allows larger clusters. Alternately, intercluster interaction terms could be treated at the mean field level, as done, for instance, in Refs. [27,41]. Interactions that do not have a density-density form (and thus not diagonal in the Wannier basis) would, naturally, complicate matters.
A legitimate question is whether other broken symmetries could compete with superconductivity in the phase diagram. We expect charge order to be a serious contender at commensurate filling (in particular n = 0.5 and n = 1.5), provided extended interactions are taken into account. It is possible that the superconducting order that we found would disappear precisely at these fillings, either because of the extended interactions or out of competition with charge order. Likewise, antiferromagnetism is likely to appear at half-filling (n = 1), where superconductivity is suppressed, because of a suppression of the density of states related to Mott physics. Again we leave this question for future work.