Quasiperiodic pairing in graphene quasicrystals

Introduction.— Quasicrystals [1] (QCs) stand as an intermediary between disordered and periodic structures, giving rise to a unique blend of physics derived from both realms [2, 3, 4]. Consequently, they exhibit captivating properties such as quantum criticality [5], confined states [6, 7, 8, 9, 10, 11], self-similarity/fractals [12, 13, 14, 15], and higher-dimensionality [16, 17, 18]. The absence of translational symmetry in QC, despite their long-range patterns, opens up many questions regarding the emergence of conventional crystalline phenomena, such as topology [19, 20, 21, 22, 23, 24, 25, 26, 27], correlated states [28, 29, 30, 31, 32, 33], and magnetism [34, 35, 36, 5, 37, 38, 38, 39]. Furthermore, despite their quasiperiodicity, QCs exhibit a trace of energy dispersion in their spectral function, suggesting the presence of partially Bloch-like wavefunctions [40, 41, 42, 43, 44], which gives rise to exotic phenomena including quantum oscillations [40] as in periodic metallic systems.

Interestingly, in line with the existence of many intriguing crystalline phenomena in QCs, superconductivity has recently been discovered also in QC materials [45, 46]. Theoretically, superconductivity in QC can arise either intrinsically [47, 48, 49, 50] or through the proximity effect [51, 52]. Despite the absence of well-defined Bloch wavefunctions, modified versions of crystalline superconductivity theories [53] were proposed to explain superconductivity in QCs [48]. However, the quasiperiodicity effect can strongly affect superconducting properties in QC [54, 42, 55, 56, 57, 58, 59, 60, 61]. For instance, superconductivity in QC can potentially give a “fractal superconductor” due to the underlying self-similar patterns [45]. Interestingly, despite the strong dependence of the order parameter on local environments [30], the superconducting gap remains uniform throughout the system with a well-defined critical temperature [62]. However, specific heat jumps at the superconducting transition and Bogoliubov peaks may be suppressed due to the underlying quasiperiodic wavefunction [54].

Recent advances in the study of twisted multilayer moiré materials have garnered new attention to extrinsic QCs [63]. Unlike intrinsic QCs with inherent quasiperiodicity, extrinsic QCs originate from stacking two-dimensional periodic lattices with relatively incommensurate basis vectors. In certain cases, this stacking leads to atomic QC structures with a quasiperiodic tiling of clusters of atoms. In particular, twisting two graphene layers with a $30^{\circ}$ angle leads to twisted bilayer graphene quasicrystal (TBGQC), which manifests 12-fold dodecagonal Stampfli tiling [64, 65]. The TBGQC can commonly appear alongside more energetically favorable AB stacked bilayer graphene [66]. Furthermore, TBGQC was experimentally confirmed by observing low-energy Dirac cones and their mirrored images due to Umklapp interlayer scattering. It was suggested that these Dirac cones may give rise to relativistic Dirac physics within QC realms [67, 68, 69]. However, it was later realized that these Dirac cone replicas do not contribute to the electrical transport [70], and effectively low energy electronic properties are described by the original graphene’s Dirac cones [71]. Moreover, huge interlayer coupling is needed to achieve quasicrystal physics at the Fermi energy [72, 73]. To avoid such unrealistic conditions, the Fermi energy must be shifted away from the Dirac point to the higher energy window, where new van Hove singularities (vHS) with quasicrystalline wavefunctions appear [63, 74]. Therefore, TBGQC hosts both periodic energy ranges (PER) and quasiperiodic energy ranges (QER) simultaneously, and quasiperiodicity can be toggled on or off by external stimuli such as pressure, and electrostatic gating [75].

In this Letter, motivated by the recent discoveries of superconductivity in extrinsic QCs [76], we study the superconducting instability of TBGQC, in a broad energy regime including both the PER and QER. We note that superconductivity has been observed in both twisted or non-twisted multi-layer graphene, particularly when the electronic dispersion is flattened [77, 78, 79, 80, 81, 82]. In the case of TBGQC, the tunability of its quasiperiodicity using Fermi energy manipulation offers a unique opportunity to study the quasiperiodicity effect on QC superconductors. It has been shown that TBGQC may host high-angular-momentum topological superconductivity for PER [83]. However, the significance and effects of quasiperiodicity on superconducting instability in QER still remain elusive. We found that for Fermi energy in PER, the pairing instability can be understood as a simple superposition of two monolayer graphene superconductors. On the other hand, in QER, the quasiperiodicity not only induces the quasiperiodic modulation of pairing amplitudes, but also enhances the superconductivity itself, especially for onsite spin-singlet pairing. Our work paves a way to understand the fundamental properties of the emergent superconductivity in various systems with inherent quasiperiodicity such as TBGQC and moiré quasicrystals [76].

Model and method.— To elucidate the irrational structure of TBGQC, let us consider a general TBG. The basis lattice vectors of the first layer are given by $\vec{a}_{1,2}=a(\tfrac{3}{2},\pm\tfrac{\sqrt{3}}{2})$, and the corresponding second layer vectors $\vec{a}^{\prime}_{1,2}$ are obtained by $\theta$ rotation, where $a$ is the shortest carbon-carbon distance. The TBG can form a commensurate structure when $\cos\theta=\frac{3p^{2}+3pq+\frac{q^{2}}{2}}{3p^{2}+3pq+q^{2}}$ holds where the integer $p$ and $q$ satisfy $p\vec{a}_{2}+(p+q)\vec{a}_{1}=(p+q)\vec{a}^{\prime}_{2}+p\vec{a}^{\prime}_{1}$  [84]. Therefore, the twisted layers make a super-lattice with larger basis vectors $\vec{A}_{1}=p\vec{a}_{2}+(p+q)\vec{a}_{1}$, $\vec{A}_{2}=(2p+q)\vec{a}_{2}-(p+q)\vec{a}_{1}$, and its super-cell containing $N_{\text{cell}}=4(3p^{2}+3pq+q^{2})$ atoms. The TBGQC is obtained by choosing $\theta=30^{\circ}$, at which $p$ and $q$ satisfy $\tau_{D}\equiv\tfrac{q}{p}=\sqrt{3}$. The irrational number $\tau_{D}$ represents the quasi-periodicity of TBGQC similar to the golden ratio in Penrose QC [17, 8].

Analyzing QC structures can be simplified by studying their corresponding periodic approximants [85]. To find different approximants of TBGQC  [86], one can approximate $\tau_{D}$ by truncating the continued fraction of $\sqrt{3}$ and obtain two coprime integers $p_{g}$ and $q_{g}$ satisfying $\tau_{D}(g)=\tfrac{q_{g}}{p_{g}}$ where $g$ is an index for approximant generations, and by definition satisfies $\lim_{g\rightarrow\infty}\tau_{D}(g)=\sqrt{3}$. For instance, selecting values like $\tau_{D}(1)=\tfrac{1}{1},\tau_{D}(2)=\tfrac{2}{1},\tau_{D}(3)=\tfrac{5}{3}$, and $\tau_{D}(4)=\tfrac{7}{4}$ results in approximants containing $N_{\text{cell}}=28,52,388,$ and $724$ atoms, respectively. In this work, our primary focus is g=4 ($\tau_{D}(4){=}\tfrac{7}{4}$), or equivalently $\theta=30.1583^{\circ}$. Fig. 1(a) shows the corresponding approximant supercell, which is large enough to capture the quasiperiodicity while ensuring computational feasibility.

To describe the energy spectrum of TBGQC, we use a tight-binding Hamiltonian that describes hopping between carbon $p_{z}$ orbitals [87] (see Supplemental Material (SM) at ¹¹1see Supplemental Material at **** which also includes citation [111, 112, 113, 106, 114, 115, 116, 117]. ). The translational symmetry of the TBGQC approximant allows us to diagonalize the Hamiltonian in the supercell momentum space, which gives the energy eigenvalues $\epsilon_{\vec{k},\lambda}$ ($\lambda=1,\dots,N_{\text{cell}}$) and the energy eigenvectors $\ket{\psi_{\vec{k},\lambda}}$. In Fig. 1(b), the average local density of states (LDOS) $\rho(\omega)\equiv\tfrac{1}{N_{\text{cell}}}\sum_{i=1}^{N_{\text{cell}}}\rho_{i}(\omega)$ is plotted with (purple plot) and without (grey plot) interlayer coupling, respectively. Here, $\rho_{i}(\omega)=\sum_{\vec{k},\lambda}\left|\innerproduct{i}{\psi_{\vec{k},\lambda}}\right|^{2}\delta(\epsilon_{\vec{k},\lambda}-\omega)$ is the LDOS for the $i$-th site in a supercell. We note that the interlayer coupling generates singular peaks in LDOS, which arise from new vHS in TBGQC, dubbed the quasicrystal vHS (QCvHS), in the energy range of $\omega\in[-2.5,-1]\text{eV}$ [63]. To examine the influence of quasiperiodicity on the wavefunction, we show [green plot in Fig. 1(b)] the standard deviation of the LDOS, $\sigma_{\rho}(\omega)\equiv\frac{1}{N_{\text{cell}}}\sqrt{\sum_{i}(\rho_{i}(\omega)-\rho(\omega))^{2}}$. $\sigma_{\rho}(\omega)$ is prominent near the newly generated QCvHS, indicating the importance of the underlying QC structure for these energy ranges. In a previous study, the inverse participation ratio (IPR) enhancement was observed in a similar energy range, indicating the non-uniformity of the corresponding wavefunctions [63]. In the subsequent discussion, we refer to the energy range $\omega\in[-2.5,-1]\text{eV}$ with substantial LDOS fluctuations as QER [indicated in Fig. 1(b)].

Since the TBGQC approximant contains many atoms, the resulting band structure is complicated and contains many folded energy bands of original graphene. However, one can unfold the band structure by calculating spectral functions [87], which captures how original graphene band structures are affected by interlayer coupling. Explicitly, the spectral function is given by $A^{l}_{\alpha}(\vec{k},\omega)=\sum_{\lambda,\alpha}\left|\innerproduct{\phi_{\vec{k},l,\alpha}}{\psi_{\vec{k},\lambda}}\right|^{2}\delta(\epsilon_{\vec{k},\lambda}-\omega)$ where $\ket{\phi_{\vec{k},l,\alpha}}$ indicate the original graphene wavefunctions with the layer index $l=\pm 1$ and the sublattice index $\alpha=A,B$. In Fig. 1(c), we plot the unfolded band structure for the TBGQC approximant, which around zero energy closely resembles the Dirac cone of the original graphene while significant modification appears in QER.

In Fig. 2(a), we plot the LDOS $\rho_{i}(\omega)$ at different energies $\omega$, which exhibits quasiperiodicity within QER, reflecting the inhomogeneous local environment distribution of lattice sites. Contrary to each graphene layer where the lattice sites belonging to a given sublattice (A or B) share the same local environment, the local environments exhibit large fluctuations in TBGQC. To clarify the relation between the local environment and the LDOS inhomogeneity, we introduce an “environment-based classification space” (EBCS) as depicted in Fig. 1(d). To obtain an EBCS, we first project all lattice sites of the top layer belonging to a specific sublattice onto the bottom layer. Then, we record the position of each projected lattice site relative to the hexagonal unit cell of the bottom layer enclosing the projected site. By marking the relative positions of all the projected sites within a single hexagon, we obtain an EBCS. Because EBCS in TBGQC contains information about local environments, it resembles the perpendicular space [89, 90] in intrinsic QCs which can classify vertices based on their local environments [91, 30]. In Fig. 2(b), we plot the LDOS distribution in EBCS, which shows smoother variation compared to the original real-space depiction in Fig. 2(a). This implies that in QER, vertices with similar local environments yield similar LDOS, thereby reflecting underlying quasiperiodicity, which also alters emergent superconducting instability in TBGQC, as illustrated below.

Pairing instability– To study the superconducting instability of TBGQC, we assume the following phenomenological effective attraction

where $U_{ij}^{\eta}=U_{ji}^{\eta}<0$ is the attraction between electrons at the $i$-th and $j$-th vertices. $\eta=+1,-1$ indicate the spin-singlet and spin-triplet pairing channels, respectively. We investigate the pairing instability using the linearized gap equation approach, which is a powerful tool for studying competing superconducting orders near the critical temperature, and has been extensively applied to various superconductors [92, 93, 94, 95, 96] including the quasicrystalline superconductivity [49].

As shown in SM [88], we develop a real-space linearized gap equation suitable for TBGQC approximant given by $M_{L,L^{\prime}}\Delta_{L^{\prime}}^{\eta}=\zeta\Delta_{L}^{\eta}$, where $M$ is the pairing matrix. Here, $L\equiv(i,j)$ represents a pair of vertices that denote the electron positions forming a Cooper pair. To construct the pairing matrix, we consider the pairs of vertices with nonzero attraction $U_{ij}^{\eta}$. The order parameter is given by a vector $\Delta^{\eta}_{L=(i,j)}=U_{ij}^{\eta}\big{<}c_{j\downarrow}c_{i\uparrow}+\eta c_{i\downarrow}c_{j\uparrow}\big{>}$ that contain local pairing information. The pairing eigenvalues ($\zeta$) and pairing eigenvectors ($\Delta$) of $M$ give the possible superconducting pairing instabilities of TBGQC. As the $\zeta$ increases, one can expect a larger critical temperature for the corresponding pairing channel. Thus, the dominant pairing instability appears in the channel with the largest $\zeta_{\text{max}}$. In deriving the linearized gap equation, we do not assume any particular symmetry or distribution of $\Delta$. This allows us to capture all possible superconducting states including the inhomogeneous pairings such as FFLO [97] or fractal superconducting states [45]. We find that the FFLO pairing instabilities always exhibit small pairing eigenvalues in our calculations, thus they are neglected in the following discussion.

Let’s first consider the onsite attraction represented by $U_{i,j}=-\delta_{i,j}$, which exclusively results in spin-singlet pairings [98, 99]. As higher DOS leads to higher superconducting transition temperatures in BCS theory [100], the LDOS enhancement due to quasiperiodicity may augment the quasiperiodic superconducting instabilities in TBGQC [101, 102, 55]. To confirm it, we consider the largest eigenvalues of the pairing matrix for decoupled and coupled graphene layers in Fig. 3(a) and (b), respectively. One can observe that superconductivity is indeed enhanced in QER when two layers are coupled. To demonstrate the relation between the LDOS enhancement and enhanced superconductivity, we compare the largest eigenvalues of $M$ and $\overline{M}$ (denoted by $\zeta_{\text{max}}$, $\overline{\zeta}_{\text{max}}$, respectively). Here $\overline{M}_{l,l^{\prime}}=\tfrac{1}{N_{\text{cell}}}\sum_{L\in l;L^{\prime}\in l^{\prime}}M_{L,L^{\prime}}$ is a $2\times 2$ pairing matrix obtained by averaging the matrix elements of $M$ with identical layer indices ($l,l^{\prime}=\pm$) [88]. Therefore, $\overline{M}$ describes uniform superconductivity within each layer, thus benefiting only from DOS but not LDOS. In Fig. 3(c), we plot $(\zeta_{\text{max}}-\overline{\zeta}_{\text{max}})/\zeta_{\text{max}}$ which shows relative enhancement of the dominant pairing eigenvalue. Notably, the pairing eigenvalue exhibits a significant enhancement within QER. In contrast, no such enhancement is observed in other energy ranges.

To systematically examine the dominant pairing profile, we compare $\Delta$ with trial pairing eigenvectors with specific pairing symmetry defined in the layer decoupled limit. We compute the pairing overlap amplitude (POA) of $\Delta$ with the trial pairing eigenvectors. If POA approaches 1 (0), the trial pairing eigenvector can (cannot) fully describe the given pairing eigenvector. Generally, POA smaller than 1 indicates the presence of non-uniform and aperiodic pairing. We denote the characteristic of the given pairing eigenvector using the symmetry of the trial pairing eigenvector with the largest POA. In Fig. 3(a) and (b), we characterize the pairing instabilities with red and blue colors indicating the s-wave and s^∗-wave symmetries [obtained in the layer decoupled limit, respectively]. Here and in the following, the pairing symmetry with (without) ^∗ symbol indicates that the pairing eigenvector has the same (opposite) signs between layers. As shown in Fig. 3(a), both s-wave and s^∗-wave symmetries are degenerate in the layer decoupled limit. On the other hand, as demonstrated in Fig. 3(b), when two graphene layers are coupled, $\Delta$ with s-wave characteristic dominates over s^∗-wave pairing. As shown in Fig. 3(f), the POA of the dominant $\Delta$ deviates significantly away from the s-wave pairing in QER while in other energy ranges, POA is one, indicating uniform s-wave pairing in both layers. We note that increasing approximant generation toward quasiperiodic structure further enhances the quasiperiodic pairing instabilities [88].

In Fig. 3(d) and (e), we depict the distribution of the dominant pairing eigenvector within a supercell and an EBCS, respectively. The dominant $\Delta$ does not change its sign across the entire system, which, in line with Fig. 3(b), can be considered as an inhomogeneous s-wave pairing. Comparing Fig. 2(a,b) with Fig. 3(d,e), it becomes apparent that $\Delta$ mirrors the underlying LDOS distribution. Similar to LDOS distribution, $\Delta$ also exhibits quasiperiodicity within QER, while it is uniform at other energy ranges (see e.g. Fig. 3 (d,e) for $\omega=0.5\text{eV}$ ). The order parameter and LDOS are distributed around their average values, remaining nonzero at all vertices, which resembles the characteristics observed in weakly disordered superconductors [103] but is distinct from highly disordered superconductors with superconducting islands [104]. The absence of superconducting islands despite the underlying quasiperiodicity indicates that superconductivity remains coherent and is not suppressed by phase fluctuations [55].

Furthermore, we consider the effect of longer-range [nearest-neighbor (NN) and next-nearest-neighbor (NNN)] attractive interactions in spin-singlet (SS) and spin-triplet (ST) channels [105] (see further detail in [88]). For monolayer graphene, due to the underlying $D_{6h}$ symmetry, p-wave ($p_{x}$, and $p_{y}$) and d-wave ($d_{xy}$, and $d_{x^{2}-y^{2}}$) pairings are doubly degenerate, while f-wave and s-wave pairings are not [106]. In TBGQC, due to the $D_{6d}$ symmetry [63], f-wave pairing also becomes doubly degenerate ($f_{y^{3}-3yx^{2}}$, and $f_{x^{3}-3xy^{2}}$), similar to p-wave ($p_{x}$, and $p_{y}$), p*-wave ($p^{*}_{x}$, and $p^{*}_{y}$), d-wave ($d_{xy}$, and $d_{x^{2}-y^{2}}$) and d*-wave ($d^{*}_{xy}$, and $d^{*}_{x^{2}-y^{2}}$) pairings.

In the NN-SS attraction, mostly d-wave and d^∗-wave pairing instabilities are competing. However, within QER d-wave pairing becomes dominant pairing instability, but its POA becomes smaller than 1 [Fig. 3(g)]. Especially, highly quasiperiodic pairing with small POA and mostly s-wave characteristic emerges at around $\mu=-1.96$ eV [see red dots in Fig. 3(g)]. In the case of the NNN-SS attraction, mostly s-wave pairing is dominant, although close to the main vHS $\mu\lessapprox-2.2$ eV, d-wave becomes dominant. As expected, within QER [see Fig. 3(h)], POA is less than one, indicating underlying quasiperiodic pairing.

For the NN-ST attraction, the dominant instability is p-wave or p^∗-wave pairing. However, near zero energy, four-fold degenerate Kekule ordering [107, 108] dominates [88]. Like other channels, POA of dominant pairings (p-wave or p*-wave) drops away from 1 in QER [Fig. 3 (i)]. Finally, in the NNN-ST attraction, f-wave pairings dominate and show no significant aperiodicity even within QER [109, 110], which may arise from maximizing f-wave pairing eigenvalue in QER as shown in SM using perturbation theory (see S3.D [88]).

Discussion.— To conclude, we have demonstrated the appearance of quasiperiodic pairing instabilities with s-wave characteristics in QER for different spin-singlet attraction channels, especially near QCvHS. This suggests that quasiperiodic s-wave pairing instability is generally expected to be dominant in TBGQC. Moreover, our theory can also be applied to the superconductivity in the recently discovered moiré graphene quasicrystal constructed by stacking three layers with small rotation. Interestingly, moiré graphene quasicrystal shows both features of moiré energy range and QER in which superconductivity emerges [76]. In the spectral function of the moiré quasicrystal, akin to TBGQC [63], many weakly dispersing states (partial flat bands) are observed within QER[76]. We speculate that the large LDOS fluctuations induced by quasiperiodicity, and the existence of partial flat bands, may be a crucial driving force behind the emergent superconductivity in moiré quasicrystals.