Metal-insulator transition and dominant $d+id$ pairing symmetry in twisted bilayer graphene

TL;DR

This study models twisted bilayer graphene using a two-orbital Hubbard model, revealing a metal-insulator transition, antiferromagnetic order, and dominant $d+id$ superconductivity, highlighting its tunability for correlated phenomena.

Contribution

It provides a detailed theoretical analysis of the electronic phases and pairing symmetry in twisted bilayer graphene, emphasizing the role of electron correlations and interlayer coupling.

Findings

Metal-insulator transition confirmed by temperature-dependent conductivity.

Dominant $d+id$ superconducting pairing over a wide filling range.

Suppression of conductivity, spin correlation, and superconductivity with increased interlayer coupling.

Abstract

Motivated by recent experimental studies that have found signatures of a correlated insulator phase and tuning superconductivity in twisted bilayer graphene, we study the temperature-dependent conductivity, the spin correlation and the superconducting pairing correlation within a two-orbital Hubbard model on an emergent honeycomb lattice. The evaluation of the temperature dependence of the conductivity demonstrates that there is a metal-insulator transition, and the Mott phase at strong coupling is accompanied by antiferromagnetic order. The electronic correlation drives a $d+id$ superconducting pairing to be dominant over a wide filling region. All of the dc conductivity, the spin correlation and the superconductivity are suppressed as the interlayer coupling strength increases, and the critical $U_c$ for the metal-insulator transition is also reduced. Our intensive numerical results reveal that twisted bilayer graphene should be a uniquely tunable platform for exploring strongly correlated phenomena.

Equations

$$H_{tb}$$

$$H_0$$

$$H_1$$

$$H_2$$

$$H_U=U\sum _{i,m}{\mathbf{n}}_{im\uparrow }{\mathbf{n}}_{im\downarrow }+V\sum _i{\mathbf{n}}_{ix}{\mathbf{n}}_{iy},$$

$${\sigma }_{dc}(T)=\frac{{\beta }^2}{\pi }{\Lambda }_{xx}(\text{q}=0,\tau =\frac{\beta }{2})$$

$$S_{AFM}=\frac{1}{N_s}⟨[\sum _{mr}({\hat{S}}_{mar}^z-{\hat{S}}_{mbr}^z){]}^2⟩,$$

$$P_{\alpha }=\frac{1}{N_s}\sum _{l,i,j}{\int }_0^{\beta }d\tau ⟨{\Delta }_{l,\alpha }^{†}(i,\tau ){\Delta }_{l,\alpha }^{}(j,0)⟩,$$

$${\Delta }_{l\alpha }^{†}(i)=\sum _{\mathbf{l}}{f}_{\alpha }^{†}({\delta }_{\mathbf{l}})(a_{li\uparrow }{b}_{li+{\delta }_{\mathbf{l}}\downarrow }-a_{li\downarrow }{b}_{li+{\delta }_{\mathbf{l}}\uparrow }{)}^{†},$$

Full text

I Introduction

In accordance with recent experiments on twisted bilayer graphene (TBG), the arresting phenomena including novel phases, unconventional superconductivity and a Mott-like insulator behavior have been discovered in an excellent two-dimensional system that demonstrates the importance of correlations effects*Yuan *et al. (2018); *Cao et al. (2018). Most recently, a tuning superconductivity is induced by varying the interlayer spacing with hydrostatic pressure, which sparked more intense interest in TBG, as it may be a uniquely tunable platform for exploring correlated statesYankowitz *et al. (2019).

With two layers of graphene twisted at a narrow range of particular magic angle, the band structure of twisted bilayer graphene becomes nearly flat, and as a result, the Fermi velocity decreases to zero in the vicinity of the Fermi energy. Being interpreted as a correlated Mott insulator at half-fillingYuan et al. (2018), this system is doped with a few additional charge carriers that then change the system from the initial insulator to a superconductorCao et al. (2018), and shares a strikingly similar trend as that in doped cupratesBednorz and Muller (1986), heavy-fermion*Steglich et al. (1979), iron-basedKamihara et al. (2008) and organic superconductorsJerome *et al. (1980). Consequently, it has been the subject of intense study since the discovery of high-temperature superconductors Scalapino (2012), which may shed light on several long-standing problems encompassing the understanding of unconventional superconductivity, and may even prove to be a significant step in the search for room-temperature superconductors.

Substantial theoretical effort has gone into this field*Rozhkov *et al. (2016); *Ray *et al. (2019a); Roy and Juričić (2019); *Kennes *et al. (2018); *Liu *et al. (2018); Zhang (2019); Su and Lin (2018); *Tang *et al. (2019); *Huang *et al. (2019); Yuan and Fu (2018a, b); *Ray *et al. (2019b); *Ochi *et al. (2018); *Chen *et al. (2018a, b); *Zhu *et al. (2018); *Koshino *et al. (2018); Kang and Vafek (2018); González and Stauber (2019); *Song *et al. (2019); *Ahn *et al. (2019); *Zou *et al. (2018); *Liu *et al. (2019); *Carter *et al. (2018); Jian and Xu (2018); *Wu *et al. (2019a); Wu (2019); Venderbos and Fernandes (2018); *Gu *et al. (2019); *Wu *et al. (2019b); *Lian *et al. (2018); Choi and Choi (2018); Xie and MacDonald (2020); *Peltonen *et al. (2018); *Mora *et al. (2019); *Chebrolu *et al. (2019); *Da Liao *et al. (2019); *Zhang *et al. (2019); *Classen *et al. (2019); *Codecido *et al. (2019); *Zhang et al. (2020), and many possibilities of the exotic electronic structures further reflect the fact that TBG can be a realistic platform for various kinds of largely unknown physics. However, it is difficult to identify an effective low-energy model to address the strong correlation effects in TBG, as the moiré pattern in TBG with small twist angles requires a very large system size and more than 10, 000 atoms in one unit cell, which makes the usual first-principles electronic structure calculations almost impossible. The origin of the insulating state and the possible pairing symmetry of the observed superconductivity are both highly debatedRay *et al. (2019a); Roy and Juričić (2019); *Kennes *et al. (2018); *Liu *et al. (2018); Zhang (2019); Su and Lin (2018); *Tang *et al. (2019); *Huang et al. (2019) due to the weakness of the proposed effective model and the method to treat strong electronic correlation, which are the major challenge in this active field. For example, some mean filed results argued that part of the experimental findings could well be understood to be a $s$-wave superconducting state resulting from an attractive electron-electron interaction mediated by electron-phonon couplingPeltonen et al. (2018) and some recent works regarded the insulating state at half filling as a non-Mott picture but maybe some Kekulé valence-bond orderXu et al. (2018), chiral spin-density wave orderLiu et al. (2018), a nematic phasesDodaro et al. (2018), or crystalline statesPadhi *et al. (2018).

Recently, a two-orbital Hubbard model, sketched in Fig.1, constructed from Wannier orbitals that extend over the size of supercells, is proposed to capture the electronic structure of narrow minibands and the effect of Coulomb interaction in TBG, and the centers of these Wannier orbitals form an emergent honeycomb latticeYuan and Fu (2018a, b). This model has been verified by explicit numerical calculations*Koshino *et al. (2018); Kang and Vafek (2018) and provides us an opportunity to explore the rich physics in TBG by the unbiased numerical method, which is the most reliable way to establish the phase diagram in the system where the strong correlation effect dominates. In this paper, we are aiming to identify the insulating state at half filling and the dominant pairing symmetry of the superconductivity with the proposed effective model.

The metal-insulator transition is one of the most fundamental and profound physical phenomenon of quantum mechanics. However, our understanding of interaction-driven metal-insulator transitions still remains rather controversial because strongly correlated systems are difficult to solve by using both analytical and numerical methods*Otsuka *et al. (2016). In this paper, with the two-orbital Hubbard model on an emergent honeycomb lattice, we perform a quantum Monte carlo study of the metal-insulator transition and pairing interaction in TBG. The calculations of the current-current correlation function show that repulsion between electrons can significantly reduce the conductivity, and at low temperatures change the system from conducting behavior to insulating behavior. After performing a careful finite size scaling analysis, we demonstrate that the metal-insulator transition is associated with the presence of antiferromagnetic (AFM) long-range order. Our intensive numerical simulations suggest that the state at half filling is a Mott-like insulating state with long range AFM order, and with a finite doping, the dominant pairing symmetry of the superconducting state is $d+id$-wave.

II Model and methods

In the effective model, the tight-binding part $H_{tb}$ mainly contains the intralayer hoppings $H_{0}$, and $H_{1}$ and $H_{2}$ are also introduced into $H_{tb}$ to further break the SU(4) and SU(1) symmetriesYuan and Fu (2018a, b). The detailed $H_{tb}$ can be written as

[TABLE]

where ${\bf c}_{i}=(c_{i,x},c_{i,y})^{\text{T}}$ with $c_{i,x(y)}$ annihilates an electron with $p_{x(y)}$-orbital at site $i$. $t_{1}$ and $t_{2}$ are the hopping amplitude between nearest-neighbor(NN) and fifth-NN sites. The sketch of hoppings is shown in the insert of Fig. 1(b). The chiral basis $c_{\pm}=(c_{x}\pm ic_{y})/\sqrt{2}$ is associated with $p_{x}\pm ip_{y}$ orbitals. $\bm{e}^{\parallel,\perp}_{ij}$ denotes in-plane unit vectors in the direction parallel and perpendicular to the NN bond $\langle ij\rangle$. The on-site Coulomb interaction part is written as

[TABLE]

where $m$ is the $p_{x(y)}$ orbitals, and ${\bf n}_{im\sigma}={\bf c}^{\dagger}_{im\sigma}\cdot{\bf c}_{im\sigma}$. We have mainly used $V$=0 in this paper, except as explicitly noted elsewhere. In this model, $U_{v}(1)$ is kept and $C_{6}$ is completely ignored. Thus, in the Dirac points at charge neutrality are not symmetry-protected robust features*Zou *et al. (2018). The unconventional superconductivity that we focused in this paper is away from the charge neutrality, and therefore, the controversial point of this effective model has no influence on our numerical calculations. Moreover, the electronic filling is controlled by the chemical potential, $H_{\mu}=\mu\sum_{i,m,\sigma}{\bf n}_{im\sigma}$. This model provides a theoretical framework for studying correlated electron phenomena in TBG. In our following simulations, the system we performed is sketched in Fig. 1(a) with periodic boundary conditions, and we take $t_{1}$ as the unit. With this effective model, the hopping parameters dictate bandwidth (at $\Gamma$ point), Dirac velocity (at $K$ point) and valley splitting (along $\Gamma-M$ lines). Such band properties can be calculated as functions of twist angle, and one could then compare them to extract hopping parameters. We can then, choose the parameters that we find at the magic angle. The parameter we used is $t_{1}=1.0,t_{1}^{\prime}=0.1,t_{2}=0.025$, and $t_{2}^{\prime}=0.1$, which is taken from that of Refs.Yuan and Fu (2018a, b). With this set of parameters, one can see that the band structure in Fig. 1(b) is degenerate with finite $t_{1}^{\prime}$ including $\Gamma$ and $K$ points.

Our simulations are mostly performed on lattice of $L$=4, and the total number of lattice sites is $N_{s}$=2$\times$2$\times$3$L^{2}$, in which the first $2$ indicates two orbits, and the second $2$ means two interpenetrating triangular sublattices with hexagonal shape such that it preserves most geometric symmetries of graphene*Ma *et al. (2011, 2014, 2010). The number 3 means that each triangular sublattice is consistent of 3 square lattices with $L^{2}$ sites. The basic strategy of the finite temperature determinant quantum Monte Carlo (DQMC) method is to express the partition function as a high-dimensional integral over a set of random auxiliary fields. The integral is then accomplished by Monte Carlo techniques. In our simulations, 8 000 sweeps were used to equilibrate the system, and an additional 30 000$\sim$ 240 000 sweeps were then made, each generating a measurement. These measurements were divided into ten bins that provide the basis of coarse-grain averages, and errors were estimated based on standard deviations from the average. To assess our results and their accuracy with respect to the infamous sign problem as the particle-hole symmetry is broken, a very careful analysis on the average of sign is shown.

To explore the phase transitions between the metal and insulator behaviors, we compute the $T$-dependent DC conductivity, which is calculated from the wave vector q- and imaginary time $\tau$-dependent current-current correlation functionTrivedi and Randeria (1995) $\Lambda_{xx}(\textbf{q},\tau)$,

[TABLE]

where $\Lambda_{xx}(\textbf{q},\tau)=\left<\hat{j}_{x}(\textbf{q},\tau)\hat{j}_{x}(\textbf{-q},0)\right>$, $\beta=1/T$, $\hat{j}_{x}(\textbf{q},\tau)$ is the $(\textbf{q},\tau)$-dependent current operator in $x$ direction. Eq.3 has been employed for metal-insulator transitions in the Hubbard model in many works, and it has already proved its validityTrivedi and Randeria (1995); *Denteneer *et al. (1999); *Ma *et al. (2018).

To examine how the system evolves with the variation in the magnetic order, we study the AFM spin structure factor

[TABLE]

which indicates the onset of long-range AFM order if $\lim_{N_{s}\rightarrow\infty}(S_{AFM}/N_{s})>$0. Here, ${\hat{S}^{z}_{mar}}$(${\hat{S}^{z}_{mbr}}$) is the $z$ component spin operator on A (B) sublattice of orbit $m$. $S_{AFM}$ for different interactions are calculated on lattices with $L=3,4,5,6$, and are extrapolated to the thermodynamic limit using polynomial functions in $1/\sqrt{N_{s}}$.

III Results and discussion

First, we examine the AFM spin structure factor behaviors versus the inverse temperature $\beta$ for the different lattice size, $L$, and interaction strength, $U$. Here, we use $t_{1}=1.0,t_{1}^{\prime}=0.1,t_{2}=0.025$, and $t_{2}^{\prime}=0.1$. In the following, we fix $t_{1}=1.0,t_{2}=0.025$, and may vary $t^{\prime}_{1}=t^{\prime}_{2}$ to explore the tunable physics. From Fig. 2(a), we can see that the AFM spin structure factor increases as the lattice size increases slightly as $U\leq 3.5$, and the AFM spin structure factor clearly increases at $U=4.0$, which indicates that it has a potential to have a long range order as $U>3.5$. To identify the accurate critical value of AFM long range order from Fig. 2(b), we further extrapolate the finite size results to the thermodynamic limit by using polynomial functions in $1/\sqrt{N_{s}}$, and it can be seen from the figure that the AFM long range order starts to appear around $U\simeq 3.6\sim 3.8$. In Fig. 2(c), we also examine the behavior of the spin correlation as the interlay coupling strength varies, which shows that the spin correlation is suppressed as $t^{\prime}_{1}=t^{\prime}_{2}$ increases, which may be understood from the shown band width in Fig. 1 (b). The increasing $t^{\prime}_{1}=t^{\prime}_{2}$ results in larger band width, and the effective interaction tends to be reduced and thus the spin correlation is suppressed.

For the finite temperature DQMC method, the notorious sign problem prevents exact results for lower temperature, higher interaction, or larger lattice for cases without particle-hole symmetry. To examine the reliability of the present data shown in Fig.2, we show the average of sign in Fig. 2(c) and (d), which is dependent on the different temperature $\beta$ at different interactions $U$ (a) and different lattice sizes (b) with the Monte Carlo parameters of 30 000 times runs. For the present results, our numerical results are reliable as one can see that the average of the corresponding sign is mostly larger than 0.70 for the $U$ from 3.0 to 4.0 with 30 000 times measurements. To obtain the same quality of data as $\left\langle sign\right\rangle\simeq 1.0$, much longer runs are necessary to compensate the fluctuations. Indeed, we can estimate that the runs need to be stretched by a factor on the order of $\left\langle sign\right\rangle^{-2}$*Blankenbecler *et al. (1981); Santos (2003); *Yang *et al. (2016). In our simulations, some of the results are obtained with more than 240 000 times runs, and thus the results for the current parameters are reliable.

Second, we examine the temperature dependence dc conductivity $\sigma_{dc}(T)$ with $L=4,5$ across several interaction strengths, as shown in Fig 3 (a) and (b). For $U<3.6$, $\sigma_{dc}$ diverges as the temperature $T$ decreases to the zero, while the conductivity curve is concave down and closed to zero with decreasing temperature for $U\geq 3.8$. These behaviors of the $\sigma_{dc}$ curve in the low temperature region suggest that there is a metal-insulator transition at $U\simeq 3.6\sim 3.8$. Moreover, the critical Hubbard interaction for AFM long range order is almost the same as that from semimetal to AFM insulator, which reflects that there is no spin liquid phase in magic angle twisted bilayer graphene. To further extract the critical interaction $U_{c}$ of metal-insulator transition and perform finite-size scaling analysis, we plot the dc conductivity $\sigma_{dc}$ as a function of $U$ with three different temperatures in Fig. 3 (c) at $L=3,4,5$. For each $L$, every three curves shown in Fig. 3 (c) intersect at one point, which reflects $U_{c}(L)$. In the range of $U<U_{c}(L)$, the conductivity increases as the temperature decreases, which demonstrates the metallic phase. The opposite situation emerges within the range of $U>U_{c}(L)$. The conductivity $\sigma_{dc}$ values at higher temperature exceed those of lower temperature for the same $U$. The system maintains an insulating phase. The extracted $U_{c}$ for different lattice size have been shown in Fig. 3 (d), and the critical interaction $U_{c}\simeq 3.66$ is almost independent on the lattice size, which ensure the $U_{c}$ for metal-insulator transition.

It is interesting to ask what the fate of the correlated insulating phases is by varying the interlayer coupling $t^{\prime}_{1}$ and $t^{\prime}_{2}$, which may shed light on the relation between the superconducting phases and insulating phases. In Fig.4 (a), the dc conductivity, $\sigma_{dc}$, versus temperature, $T$, computed at various interaction strengths for $t^{\prime}_{1}=t^{\prime}_{2}=0.05$ is shown, which indicates that the metal-insulating transition is in the region of $U_{c}=3.8\sim 4.0$, which is slightly larger than that of $t^{\prime}_{1}=t^{\prime}_{2}=0.10$. The dc conductivity for different interlayer coupling at $U=3.6$ is shown in Fig.4 (b), which clearly shows that the dc conductivity reduces as the interlay coupling increases, and that a smaller $U_{c}$ for metal-insulator transition is required for larger interlay coupling. Even the band width decreases as the $t^{\prime}_{1}=t^{\prime}_{2}$ decreases, while the introduction of $t^{\prime}_{1}$ and $t^{\prime}_{2}$ may open a larger band gap, which favors the insulating phases. The simulations that presented here is from an effective model, and one possible shortcoming here is that, it is not easy to explore the effect of the twist angle or interlayer coupling on the metal-insulator transition. It would be interesting to establish the correspondence with these physical parameters of the effective model, which allows, for instance, to discuss the effect of approaching the magic twist angle. However, due to the limitation of the DQMC method of lattice size or the sign problem, it should be hard to extend that in current paper.

To identify the pairing symmetry of superconducting state in TBG, we studied the effective pairing interaction with different pairing symmetries as a function of the electronic fillings for $t^{\prime}_{1}=t^{\prime}_{2}=0.10$ in Fig.5. Following previous literatures*White *et al. (1989a); *Ma *et al. (2013); *White *et al. (1989b), the effective pairing interaction ${\bf P_{\alpha}}$, is the difference between the pairing susceptibility $P_{\alpha}$ and the bubble contribution$\widetilde{P}_{\alpha}$, and it is defined as ${\bf P_{\alpha}}=P_{\alpha}-\widetilde{P}_{\alpha}$ with

[TABLE]

where $\alpha$ stands for the pairing symmetry and the corresponding order parameter reads

[TABLE]

with $f_{\alpha}(\bf{\delta}_{\bf l})$ being the form factor of the pairing function. To extract the intrinsic pairing interaction in the finite system, one should subtract from $P_{\alpha}$ its uncorrelated single-particle contribution $\widetilde{P}_{\alpha}$, which is achieved by replacing $\langle a_{{li}\downarrow}^{{\dagger}}a_{{lj}\downarrow}b_{i+\delta_{\bf l}\uparrow}^{{\dagger}}b_{j+\delta_{\bf l^{\prime}}\uparrow}\rangle$ in Eq. (6) with $\langle a_{{i}\downarrow}^{{\dagger}}a_{{j}\downarrow}\rangle\langle b_{i+\delta_{\bf l}\uparrow}^{{\dagger}}b_{j+\delta_{\bf l^{\prime}}\uparrow}\rangle$.

In Fig.5(a), it is shown that for the investigated filling region, the pairing with $d+id$ symmetry dominates over pairings with other symmetry, and agrees with our previous results on the original geometry*Huang *et al. (2019). The positive effective pairing interaction indicates that there is indeed the possibility of electronic correlation driven superconductivity. For the exact paring form, we refer readers to Ref.*Huang *et al. (2019). In Fig.5 (b), the filling dependent effective pairing interaction with $d+id$ symmetry is shown for different $t^{\prime}_{1}=t^{\prime}_{2}$. The effective pairing interaction with $d+id$ symmetry, $P_{d+idNN}$ is suppressed as $t^{\prime}_{1}=t^{\prime}_{2}$ increases. In Figs.5(c) and (d), we compare the effective pairing interaction with $d+id$ symmetry on different lattices. One can see that, the behavior that the effective pairing interaction increasing as the temperature decreases is robust, indicating the possibility of superconducting transition in real systems.

In Fig.6(a), we examine the effect of the inter-orbital Coulomb interaction on the ground state pairing correlation, $C_{\alpha}(r=\textbf{R}_{i}-\textbf{R}_{j})=\Sigma_{l}\langle\Delta^{{\dagger}}_{l\alpha}(i)\Delta_{l\alpha}(j)\rangle$, where one can see that the pairing correlation is slightly suppressed by the inter-orbital Coulomb interaction while the leading pairing symmetry does not change. To learn more on the $d+id$ pairing correlation in the thermodynamic limit, we examine the evolution of $C_{d+id}$ with increasing lattice size. In the inset of Fig. 6(b), the average of long-range $d+id$ pairing correlation, $\overline{C}_{d+id}$=$\frac{1}{\sqrt{N^{\prime}}}\sum_{r>4a}C_{d+id}(r)$, where $N^{\prime}$ is the number of electronic pairs with $r>4a$, is plotted as a function of $\frac{1}{\sqrt{N_{s}}}$ for $U$=$3.0$. We observe that $\overline{C}_{d+id}$ has a finite positive value. This result suggests the possible presence of long-range $d+id$ superconducting order in the parameter regime investigated.

In summary, within an effective two-orbit model for TBG, we study the spin correlation, the dc conductivity and the superconducting pairing interaction by using nonbiased DQMC method. At half filling, an antiferromagnetically ordered Mott insulator is proposed beyond a critical $U_{c}=3.6\sim 3.8$. By varying $t^{\prime}_{1}=t^{\prime}_{2}$, we report a tunable metal-insulator transition and superconductivity in TBG, where the dc conductivity is suppressed as $t^{\prime}_{1}=t^{\prime}_{2}$ increases, and a smaller $U_{c}$ for metal-insulator transition is required. With a finite doping, the pairing with $d+id$ symmetry dominates over the other pairing symmetries, and it could be suppressed by increasing the interlay coupling strength close to half filling. Our exact numerical results demonstrate that the TBG holds a very similar interaction driven phase diagram of doped cuprates and other high temperature superconductors, and the TBG opens up exciting opportunities to explore tunable correlated states in two-dimensional moiré superlattice heterostructures.

Acknowledgement* — This work was supported by NSFC (Nos. 11774033 and 11974049) and Beijing Natural Science Foundation (No.1192011). The numerical simulations in this work were performed at HSCC of Beijing Normal University and Tianhe in Beijing Computational Science Research Center. *