Superfluid weight and Berezinskii-Kosterlitz-Thouless transition temperature of twisted bilayer graphene

Powered by TCPDF (www.tcpdf.org) This material is protected by copyright and other intellectual property rights, and duplication or sale of all or part of any of the repository collections is not permitted, except that material may be duplicated by you for your research use or educational purposes in electronic or print form. You must obtain permission for any other use. Electronic or print copies may not be offered, whether for sale or otherwise to anyone who is not an authorised user. Julku, A.; Peltonen, T. J.; Liang, L.; Heikkila, T. T.; Torma, P.

(Received 17 June 2019; revised manuscript received 20 November 2019; accepted 3 January 2020; published 24 February 2020) We study superconductivity of twisted bilayer graphene with local and nonlocal attractive interactions. We obtain the superfluid weight and Berezinskii-Kosterlitz-Thouless (BKT) transition temperature for microscopic tight-binding and low-energy continuum models. We predict qualitative differences between local and nonlocal interaction schemes which could be distinguished experimentally. In the flat-band limit where the pair potential exceeds the band width we show that the superfluid weight and BKT temperature are determined by multiband processes and quantum geometry of the band. DOI: 10.1103/PhysRevB. 101.060505 Recent experimental discoveries of superconductivity in bilayer graphene twisted close to a "magic angle" θ * [1-3] call for a reconsideration of traditional theories of superconductivity [4,5], in particular because the superconductivity occurs in a regime where the noninteracting electronic states form an asymptotically flat (dispersionless) band [6][7][8][9][10][11][12][13][14][15][16][17]. As the system is two-dimensional, the transition to superconductivity is bound to occur at the Berezinskii-Kosterlitz-Thouless (BKT) temperature T BKT [18][19][20] which can be determined from k B T BKT = π 8 √ det[D s (T BKT )] [21,22]. Here D s is the superfluid weight that yields the size of the supercurrent for a given phase gradient of the order parameter. In conventional theory of superconductivity [23], D s is proportional to the group velocity of electronic bands around the Fermi level. Thus D s = 0 for a flat band, and superconductivity in twisted bilayer graphene (TBG) appears puzzling. One might argue it to be due to the bands not being perfectly flat; however, we show here that a more likely explanation goes beyond the conventional theory. Here we calculate T BKT for TBG as a function of the superconducting order parameter and filling. We use two models of TBG including both the flat and a number of dispersive bands and show that superconductivity in the flat-band regime has essentially a quantum geometric origin.
Recently, it was found that D s has, in addition to the conventional contribution proportional to group velocity, a geometric contribution arising from multiband processes [24][25][26][27][28]. In a flat-band limit the geometric contribution dominates and is bounded from below by the band Berry curvature [27] and Chern number [24] [see also the discussion in the Supplemental Material (SM) [29] and Ref. [30] therein]. Here we show * tero.t.heikkila@jyu.fi † paivi.torma@aalto.fi that the geometric contribution dominates D s and T BKT in the flat-band regime of TBG. Importantly, we show that including only the few flat bands is not sufficient but one needs also a number of dispersive bands to correctly predict the geometric contribution. Therefore, approximate models of TBG such as those with only flat bands, as used for deriving upper [31] and lower [32] bounds of the superfluid weight and in many other works [33][34][35][36][37][38][39][40][41][42][43][44][45], may not be suited for quantitative predictions of TBG superconductivity. Moreover, we predict that, in the flat-band regime, local (s-wave) and nonlocal interactions yield distinct behavior, namely, an anisotropic superfluid weight in the latter case. We propose a four-terminal radio frequency spectroscopy experiment that can detect the possible anisotropy and thus distinguish between the two pairing mechanisms.
Theoretical models. In the renormalized moiré lattice method (RM), we deploy the Fermi-Hubbard Hamiltonian as [29] iασ c jβσ is the kinetic term, N is the total particle number operator, and H int is the effective FIG. 1. (a) The moiré superlattice of TBG depicted with a twist angle θ and the choice of the x and y axes. (b), (c) Single-particle energy band structures of the RM and DP methods, respectively, plotted within the moiré Brillouin zone along the path connecting the high symmetry points shown in the inset of (b). In the DP model (c) the bands coming from the valley K (K ) are drawn as solid (dashed) lines.
attractive interaction described below. Here c iασ annihilates a fermion in the αth lattice site of the ith moiré superlattice unit cell with spin σ ∈ {↑, ↓}, μ is the chemical potential, and the hopping t iα jβ includes both the intra-and interlayer terms.
To reduce the number of lattice sites M within a moiré unit cell (around 12 000 for twist angle θ ∼ 1 • ), we apply a rescaling trick [56,63] under which the Fermi velocity of a monolayer graphene and the moiré periodicity remain invariant but θ becomes larger and thus reduces M. In our computations we use the rescaling such that M = 676 and the rescaled angle is θ = 4.41 • [29] which reproduces the four narrow bands of the bandwidth of 10 meV found experimentally with θ ∼ 1 • [see Fig. 1 In the Dirac point continuum method (DP) we employ the low-energy [11,29,51] Dirac point approximation for the two graphene layers as where c σ ρ,ls (k) is the annihilation operator for spin σ , valley ρ ∈ {+, −}, layer l, sublattice s, and wave vector is the Fourier component [82] of a Slater-Koster [83] parametrized interlayer potential (times an exponential factor), K = R(θ )K − K is the difference vector from the graphene K point to its rotated counterpart, and v F is the graphene Fermi velocity. The k sum is over the the moiré Brillouin zone and the G, G sums are over the (truncated) reciprocal superlattice.
We then write the total Hamiltonian as H = where N is the total particle number operator. To describe the superconducting state with a local pairing interaction λ we use H int = λ ls dr rψ † ↑ρ,ls (r)ψ † ↓ρ,ls (r)ψ ↓ρ,ls (r)ψ ↑ρ,ls (r), which is treated in the mean-field level [29]. Hereρ is the opposite valley of ρ and ψ σ ρ,ls (r) is the continuum electron field operator.
To determine the superfluid weight D s , we first solve order parameters from the BCS gap equations [29]. In Figs. 2(a) and 2(b) we show the spatial profiles for the local and RVB interactions computed with the RM method at ν ≈ −2. Here J is chosen such that the maximum value of the order parameter is max | | ≈ 3.4 meV. From Figs. 2(c) and 2(d) we see max | | depending almost linearly on the interaction constant, which is typical for generic flat-band systems [4,5,24,25,51]. From the obtained order parameter values one can compute D s . For easier comparison between the RM and DP models, below we use max | | as a "parameter." For clarity, note that the used interaction strengths are much larger than the actual energy gap max | |, and they are of the order of the energy separation to the dispersive bands [see Figs. 2(c) and 2(d)].
To obtain D s we use linear response theory. In the meanfield level [84,85] the zero-frequency, long-wavelength limit of the current-current response function In Ref. [27] this was computed for a generic multiorbital lattice geometry with the local interaction. The details on how D s μν is obtained for our different models are discussed in the SM [29]. In Fig. 2(e) we present D s as a function of max | | at ν = −2 for both local (obtained with RM and DP) and RVB interactions (only RM). Figures 2(a), 2(b), and 2(e) depict a striking distinction between the local and RVB pairing schemes related to the pairing symmetry and the resulting form of D s . The local pairing, yielding an s-wave symmetry, conserves the underlying C 3 symmetry of the TBG lattice [ Fig. 2(a)] and D s is isotropic [29], i.e., Inset of (e) shows the total density of states (DOS) for RVB (blue curve) and local interaction (red) at max | | ≈ 3.4 meV computed with RM. The dashed curve is the DOS for local interaction obtained with DP at max | | ≈ 3.5 meV. From the DOS we see the nematic phase being gapless, while the s-wave state is gapped. The RM results are evaluated at T ≈ 0.1 K, whereas the DP results at T = 0.
interaction breaks the C 3 -rotational symmetry and yields a nematic pairing pattern in real space [ Fig. 2(b)] which leads to an anisotropic response, i.e., D s xx = D s yy and D s xy = D s yx = 0. The s wave is gapped, whereas the nematic phase has nodal points in the moiré Brillouin zone [see also the inset of Fig. 2(e)]. The anisotropic D s results in an anisotropic kinetic inductance of TBG, and it can in principle be accessed via radio frequency impedance spectroscopy [86] in a Hall-like four-probe setup.
BKT-transition temperature. By computing D s , one can determine T BKT . In Fig. 3(a) we show T BKT as a function of max | |. We can distinguish two qualitatively different regimes: In the weak-coupling limit the RVB and local interactions yield similar T BKT , whereas for stronger interactions T BKT depends on the pairing model. Moreover, around max | | 2 meV the behavior of the T BKT curves is almost linear, in accordance with previous studies [24,25] where D s of a flat band with the local interaction was shown to depend linearly on the pairing strength. In our case the narrow bands are not exactly flat but slightly dispersive and thus their flat-band characteristics manifest only when the interaction strength is sufficiently large [51]. Because of this, we call the regime with max | | 2 meV as the flat-band limit. In this regime the DP and RM results are in agreement, whereas for weak interactions the results differ due to different band structures.
The difference of the two interaction schemes is further highlighted in Fig. 3(b) which presents the ratio k B T BKT / max | (T = 0)|. At the flat-band limit this ratio approaches a constant whose value depends on the pairing potential. In experiments one can measure T BKT and in principle also deduce [from the local density of states (LDOS)] and thus the ratio of these two quantities can be used to characterize the SC pairing observed in experiments.
In Figs. 3(c) and 3(d) we present T BKT as a function of ν. The weak-coupling regime shows a dome-shaped structure of T BKT which reaches its maxima near the half-fillings of the hole-and electron-doped regimes, similar to experiments [3]. In the RM model the hole-doped region is much stronger due to higher density of states at negative energies [see Fig. 1(b)], while the DP model exhibits approximate electron-hole symmetry. The strong asymmetry of the RM model is due to the applied rescaling approximation which amplifies the finite but small asymmetry of the unscaled model [29]. In the flat-band limit, the shape of the one-particle dispersions are, except for the pronounced particle-hole asymmetry of the RM model, completely dissolved. Geometric contribution. One can decompose D s to conventional, D s conv , and geometric, D s geom , parts, so that D s = D s conv + D s geom [24,27,29]. The conventional term depends on the inverse of the effective mass of the Bloch bands and is thus a single-band contribution, whereas D s geom is a multiband effect depending on the overlap of the Bloch states and their momentum derivatives of the form ∂ k n|m , where |m are the single-particle states of the mth Bloch band and n = m [27], i.e., D s geom = 0 for a single-band system. For a strictly flat band, D s conv = 0 so its superconductivity is purely a multiband process characterized by a finite D s geom . This raises an intriguing question related to the TBG system: How much do the interband terms between dispersive and narrow bands affect D s via D s geom ? We study this question in Fig. 4 for RVB [ Fig. 4(a)] and local [ Fig. 4(b)] pairing by presenting the total D s and its components. We further show results obtained by taking into account [29] either only the four flat bands or the eight lowest (four flat, four dispersive) bands labeled as D s 4 and D s 8 , respectively. In both pairing cases the contribution coming from the four flat bands only is relatively small for larger interactions. The contribution of eight bands is larger due to a larger D s geom , which is caused by the interband terms between dispersive and flat bands as the terms between the dispersive bands only are negligible. The slight dispersion of the narrow bands results in a finite D s conv . Note that D s conv ≈ D s 4,conv , i.e., D s − D s 4,conv gives the total D s geom . From Fig. 4 we see that, at max | | ∼ 1 . . . 2 meV, i.e., when the system enters the flat-band regime, D s geom surpasses D s conv for local pairing and becomes significant in the RVB case. An important implication of Fig. 4 is the importance of the dispersive bands when computing D s and the insufficiency of four-band models even if the pairing occurs predominantly within the four narrow bands. Also for a noninteracting system, higher bands have been argued to be necessary, but for different (symmetry) reasons [89].
Discussion. Our work shows that TBG is characterized by two distinct superconducting regimes. When is much smaller than the flat-band bandwidth, the superfluid weight D s and the BKT transition temperature T BKT are well described by conventional theory of superconductivity. On the other hand, in this weak-coupling regime the results are somewhat different for the RM and DP models. This is consistent with the low-energy dispersion in TBG being very sensitive to the details of the model used [90]. In the flat-band regime where is larger than the width of the significant density of states in the flat bands, a major contribution to the superfluid weight D s originates from the geometric properties of the bands. The geometric contribution D s geom is proportional to the quantum metric [24] whose importance in physics has been recently emerging [24][25][26][27][28][91][92][93][94][95][96][97][98][99][100][101][102][103][104]. Moreover, in the flat-band regime, both D s and T BKT depend sensitively on the pairing mechanism, but not strongly on the employed microscopic model. In particular, for a non-local RVB interaction D s becomes anisotropic, which could be seen in four-terminal radio frequency spectroscopy experiments to reveal information about the pairing mechanism more directly than the LDOS measurements of [105][106][107].
Within both of our models, at θ * the crossover between the two regimes takes place for = 1 . . . 2 meV, implying T BKT ≈ 1.5 . . . 3 K. This is also the ballpark of the experimentally accessed critical temperatures [2,3]. Thus the geometric contribution of the superfluid weight and the dependence on the pairing mechanism should be relevant for current experiments. An interesting future direction of research is to include other interaction channels than pairing and explore the insulating states observed in TBG [2, 3,108]. Based on our results, one can anticipate that quantum geometry and multiband processes are important in superconductivity and correlated states of other twisted multilayer materials [42,[109][110][111][112][113][114][115][116][117][118][119][120][121][122].
Note added. Recently, a related work [123] appeared at the arXiv preprint server.