Multi-flavor Dirac fermions in Kekul\'{e}-distorted graphene bilayers

TL;DR

This paper predicts a new class of Kekul'e-textured graphene bilayers with multiple valleys and tunable Dirac quasiparticles, opening pathways for strongly correlated electronic states.

Contribution

It introduces a novel approach to engineer graphene's band structure via layer stacking, resulting in multiple valleys and tunable Dirac quasiparticles.

Findings

Up to six valleys in the band structure.

Room-temperature Dirac quasiparticles with tunable masses.

Fermi velocities half of pristine graphene, indicating strong coupling.

Abstract

Graphene's electronic structure can be fundamentally altered when a substrate- or adatom-induced Kekul\'e superlattice couples the valley and isospin degrees of freedom. Here, we show that the band structure of Kekul\'e-textured graphene can be re-engineered through layer stacking. We predict a family of Kekul\'e graphene bilayers that exhibit band structures with up to six valleys, and room-temperature Dirac quasiparticles whose masses can be tuned electrostatically. Fermi velocities half as large as in pristine graphene put this system in the strongly coupled regime, where correlated ground states can be expected.

Equations

$$H_0=\sum _{\mathbf{R}}\sum _{j=1}^3\left[t_{\mathrm{B}}(\mathbf{R},{\mathbf{\delta }}_j)a_{\mathbf{R}}^{†}{b}_{\mathbf{R}+{\mathbf{\delta }}_j}+t_{\mathrm{T}}(\mathbf{R}-{\mathbf{\delta }}_j,{\mathbf{\delta }}_j)c_{\mathbf{R}-{\mathbf{\delta }}_j}^{†}{d}_{\mathbf{R}}+\text{H.c.}\right],$$

$$\begin{array}{rl}{t}_{\mathrm{B}}(\mathbf{r},{\mathbf{\delta }}_j)=& -{\gamma }_0[1+2{\Delta }_0\cos (\mathbf{G}\cdot \mathbf{r})],\\ t_{\mathrm{T}}(\mathbf{r},{\mathbf{\delta }}_j)=& -{\gamma }_0[1+2{\Delta }_0\cos (\mathbf{G}\cdot [\mathbf{r}+{\mathbf{\delta }}_j])],\end{array}$$

$$\begin{split}H_{1}=\sum_{\mathbf{R}}\Big{[}&\sum_{j,m=1}^{3}\tfrac{1-\delta_{j,m}}{2}t_{3}(\mathbf{R}-\bm{\delta}_{j},\bm{\delta}_{j}+\bm{\delta}_{m})c_{\mathbf{R}-\bm{\delta}_{j}}^{\dagger}b_{\mathbf{R}+\bm{\delta}_{m}}\\ &+\gamma_{1}a_{\mathbf{R}}^{\dagger}d_{\mathbf{R}}+\text{H.c.}\Big{]}\end{split}$$

$$t_3(\mathbf{r},\mathbf{\delta })=-{\gamma }_3\left[1-2{\Delta }_3\cos ([2{\mathbf{K}}_{+}+{\mathbf{K}}_{-}]\cdot \mathbf{\delta }+\mathbf{G}\cdot \mathbf{r})\right],$$

$${\Delta }_3=\frac{2{\beta }_{\perp }}{{\beta }_{\parallel }}{\left[1+\frac{3h^2}{a^2}\right]}^{-1}{\Delta }_0.$$

$$H=\sum _{\mathbf{k}\in \mathrm{KBZ}}\left(\begin{array}{cc}{\psi }_{\mathbf{k}}^{†},& {\chi }_{\mathbf{k}}^{†}\end{array}\right)\left(\begin{array}{cc}{h}_{\psi }(\mathbf{k})& t(\mathbf{k})\\ t^{†}(\mathbf{k})& h_{\chi }(\mathbf{k})\end{array}\right)\left(\begin{array}{c}{\psi }_{\mathbf{k}}\\ {\chi }_{\mathbf{k}}\end{array}\right),$$

$$\begin{array}{rl}& h_{\psi }(\mathbf{k})=-{\gamma }_3\left(\begin{array}{cccc}0& s_{-1}^{\ast }& 0& -{\Delta }_3{s}_0^{\ast }\\ s_{-1}& 0& -{\Delta }_3{s}_0& 0\\ 0& -{\Delta }_3{s}_0^{\ast }& 0& s_1^{\ast }\\ -{\Delta }_3{s}_0& 0& s_1& 0\end{array}\right),h_{\chi }(\mathbf{k})={\gamma }_1\left(\begin{array}{cccc}0& 1& 0& 0\\ 1& 0& 0& 0\\ 0& 0& 0& 1\\ 0& 0& 1& 0\end{array}\right),t(\mathbf{k})=-{\gamma }_0\left(\begin{array}{cccc}0& {\Delta }_0{s}_{-1}& 0& s_{-1}\\ {\Delta }_0{s}_{-1}^{\ast }& 0& s_{-1}^{\ast }& 0\\ 0& s_1& 0& {\Delta }_0{s}_1\\ s_1^{\ast }& 0& {\Delta }_0{s}_1^{\ast }& 0\end{array}\right).\end{array}$$

$$\begin{array}{rl}{H}_{\mathrm{eff}}(\mathbf{k})=& -v({\sigma }_1{\tau }_3{k}_x+{\sigma }_2{\tau }_0{k}_y)+(3{\gamma }_3{\Delta }_3+2w^2{\Delta }_0^2{k}^2){\sigma }_1{\tau }_1\\ & -w^2(1+{\Delta }_0^2)\mathrm{Re}\{(k_x+i{k}_y{)}^2\}{\sigma }_1{\tau }_0\\ & +w^2(1+{\Delta }_0^2)\mathrm{Im}\{(k_x+i{k}_y{)}^2\}{\sigma }_2{\tau }_3.\end{array}$$

$$H_{\mathrm{eff}}(\eta {\mathbf{k}}_{\mathrm{D}}+\mathbf{q})\approx \hslash \left(\eta V_x{q}_x{\sigma }_1+V_y{q}_y{\sigma }_2\right),$$

$$\Delta V_{\mathrm{eff}}=\frac{e\Delta V}{2}\left[1-\frac{\left(1-\frac{{\gamma }_4^2}{{\gamma }_0^2}\right)\left(1+{\Delta }_0^2\right)}{{\gamma }_1+\frac{e^2\Delta V^2}{4{\gamma }_1}}{w}^2{k}^2\right]{\sigma }_3{\tau }_0,$$

Full text

Superlattices naturally emerge in multilayered two-dimensional (2D) crystals, either by differences in lattice constants or layer orientation. Heterostructures formed by stacking nearly commensurate van der Waals crystalsGeim and Grigorieva (2013) display moiré patterns, with approximate periodicities tens of times larger than the materials’ lattice constants. This long-range periodicity has enabled the observation of Hofstadter’s butterfly in graphene on hexagonal boron nitride (hBN)Hunt et al. (2013); Dean et al. (2013), and more recently of exciton minibandsAlexeev et al. (2019); Tran et al. (2019); Xia (2019); Jin et al. (2019) and moiré-umklapp optical signaturesAlexeev et al. (2019) in hetero-bilayers of transition-metal dichalcogenides. Moiré patterns also form in twisted graphene bilayers, where the conduction and valence bands flatten at so-called “magic angles” Bistritzer and MacDonald (2011), leading to strongly correlated phenomena such as unconventional superconductivity and many-body insulating behaviorCao et al. (2018a, b).

A different type of superlattice, caused by Kekulé distorsions, has also been theorizedChamon (2000); Cheianov et al. (2009) and recently reported in graphene monolayers on copper substrates at room temperatureGutiérrez et al. (2016). In the latter case, so-called ghost copper adatoms distributed periodically beneath the graphene plane produce a bond density wave that triples the graphene unit cell over nanometre-sized regions, where charge carriers comprise two different species of valley-pseudospin-locked Dirac fermionsGamayun et al. (2018); Andrade et al. (2019). One may envision Kekulé-distorted graphene as a new building block for multilayered 2D materials, with new electronic properties that can be engineered by controlling the stacking type. However, whether the substrate effect causing the Kekulé distortion will propagate to higher layers is, to our knowledge, still an open question.

In this Letter, we show theoretically that Kekulé-distorted graphene bilayers support up to six flavors of valley-degenerate Dirac quasiparticles, with electrostatically tunable relativistic masses. While these low-energy Dirac bands are directly related to the band topology transition below energies of order $1\,{\rm meV}$ in bilayer grapheneMcCann and Fal’ko (2006), we predict band widths of order $20\,{\rm meV}$, preserving multi-valley physics almost at room temperature, and for carrier dopings of up to $10^{12}{\rm cm}^{-2}$.

We propose a system where two graphene layers with Kekulé distortions are stacked with $AB$ (Bernal) configuration, for instance, following the procedure depicted in Fig. 1. Such a system may be produced by combining the fabrication of Kekulé-textured graphene on copper samplesGutiérrez et al. (2016) with state-of-the-art manual transfer techniques used to fabricate van der Waals heterostructuresCao et al. (2018a, b); Alexeev et al. (2019); Tran et al. (2019); Xia (2019); Jin et al. (2019). Without relative rotation of the layers, and by considering the possible Kekulé-center stacking types, we predict a family of $5$ patterns that may be catalogued based on their approximate 2D symmetries by the wallpaper groups $p6m$, $cmm$, $p6$, and $cm$ (Supplement). This allows for a rich playground to explore quantum phases.

Among these, we consider the low-energy bands of bilayer superlattices with $p6m$ wallpaper symmetry, which exhibit a multi-valley band structure with six flavors of massless Dirac fermions. The corresponding states are delocalized across the bilayer, such that a mass gap can be opened and tuned by out-of-plane electric fields, by contrast to Dirac fermions in monolayer graphene. The Dirac quasiparticles in the $p6m$ Kekulé bilayer are six times slower than those in single graphene layers, and thus much more susceptible to many-body phenomena driven by Coulomb interactions, such as spontaneous symmetry breakingLemonik et al. (2010); Semenoff (2012) and many-body insulating phasesVelasco Jr et al. (2012).

Model. The envisioned system may be fabricated from a single layer of epitaxial graphene grown on a copper substrate, cut and folded such that two regions with Y-shaped Kekulé distortions (Kek-Y)Gutiérrez et al. (2016) become stacked, forming interlayer bonds. One of the possible outcomes, which will be the focus of this Letter, is sketched in Fig. 1, where the stacking is $AB$ (Bernal) type, and the Kek-Y distortion centers of both layers coincide. Due to the folding, the two Kek-Y distorsions are oriented with a relative angle of $180^{\circ}$, yielding a structure whose in-plane symmetry is described by the $p6m$ wallpaper group. In three dimensions, this configuration is centrosymmetric, resulting in a spin-degenerate band structure in the absence of external magetic fields. Thus, neglecting spin-orbit interactions induced by the substrate, a spinless Hamiltonian will be considered in the following.

We model the graphene monolayers as

[TABLE]

where the annihilation operators $a_{\mathbf{r}}$ and $b_{\mathbf{r}}$ ($c_{\mathbf{r}}$ and $d_{\mathbf{r}}$) correspond to the $A$ and $B$ sublattices of the bottom (top) layer, respectively; $\bm{\delta}_{j}$ are the nearest-neighbor vectors $\bm{\delta}_{1}=\tfrac{a}{\sqrt{3}}(\tfrac{\sqrt{3}}{2},\,\tfrac{1}{2})$, $\bm{\delta}_{2}=\tfrac{a}{\sqrt{3}}(-\tfrac{\sqrt{3}}{2},\,\tfrac{1}{2})$ and $\bm{\delta}_{3}=\tfrac{a}{\sqrt{3}}(1,\,0)$, with $a=2.461\mbox{\AA}$ the lattice constant; and $\mathbf{R}=m\bm{\mathbf{a}}_{1}+n\bm{\mathbf{a}}_{2}$ ($\bm{\mathbf{a}}_{1}=\bm{\delta}_{3}-\bm{\delta}_{1}$, $\bm{\mathbf{a}}_{2}=\bm{\delta}_{3}-\bm{\delta}_{2}$; $m$ and $n$ integers) are the bottom-layer $A$ sublattice sites. In writing Eq. (1), we ignore intralayer sublattice asymmetries; the role of other on-site energy asymmetries in breaking the $p6m$ symmetry will be discussed belowMucha-Kruczyński et al. (2010).

In Bernal stacking, the top-layer lattice can be obtained from the bottom one, by an in-plane displacement of $-\bm{\delta}_{3}$, and an out-of-plane shift of magnitudeBaskin and Meyer (1955) $h=3.35\mbox{\AA}$. The Kek-Y deformations of the two layers are accounted for by the modulated hopping termsGamayun et al. (2018) (see Supplement)

[TABLE]

where $\gamma_{0}=3.16\,{\rm eV}$ is the intralayer hopping energy, $\Delta_{0}$ is the bond strength modulation due to the Kek-Y distorsion, $\mathbf{G}=\mathbf{K}_{+}-\mathbf{K}_{-}$ is a Kekulé Bragg vector, and $\mathbf{K}_{\pm}=\tfrac{4\pi}{3a}(\pm\tfrac{1}{2},\tfrac{\sqrt{3}}{2})$ are inequivalent graphene Brillouin zone (GBZ) corners (Fig. 2a).

Based on the Slonczewski-Weiss-McClure model for graphiteMcClure (1957); Slonczewski and Weiss (1958), we describe interlayer bonding with the Hamiltonian111Here, we ignore the interlayer intra-sublattice hopping term $\gamma_{4}$, which merely introduces particle-hole asymmetryMcCann and Fal’ko (2006). We show in the Supplement that all our results remain valid when that term is included.

[TABLE]

Here, $\gamma_{1}=0.381\,{\rm eV}$ describes interlayer hopping between Kek-Y centers in the two layers, which is constant throughout the superlattice. By contrast, terms of the form (Supplement)

[TABLE]

with $\gamma_{3}=0.38\,{\rm eV}$, produce the interlayer bond modulation shown in Fig. 2b. Projected onto the plane, the pattern resembles an O-shaped Kekulé distortion (Kek-O, Gamayun et al. (2018); Chamon (2000)), geometrically analogous to what occurs naturally in benzene rings. Although Kek-O patterns have not been observed in graphene, but only in artificial graphene-like systemsGomes et al. (2012), here, a Kek-O superlattice arises naturally for the $\gamma_{3}$ bonds as the separation between top-layer $A$- and bottom-layer $B$-sublattice atoms are modified by the in-plane Kekulé distortions. This is described by $\Delta_{3}$ in Eq. (4). Considering in- and out-of-plane Grüneisen parameters $\beta_{\parallel}$ and $\beta_{\perp}$, respectively, whereSun et al. (2018); Amorim et al. (2016); Naumis et al. (2017) $\beta_{\perp}\approx 1.35\beta_{\parallel}\approx 2.26$, we find

[TABLE]

Fourier-transforming the full Hamiltonian $H=H_{0}+H_{1}$, and folding it into the Kekulé Brillouin zone (KBZ, Fig. 2a), we obtain

[TABLE]

where $\psi_{\mathbf{k}}=(c_{\mathbf{k}-\mathbf{G}},\,b_{\mathbf{k}-\mathbf{G}},\,c_{\mathbf{k}+\mathbf{G}},\,b_{\mathbf{k}+\mathbf{G}})^{T}$ and $\chi_{\mathbf{k}}=(a_{\mathbf{k}+\mathbf{G}},\,d_{\mathbf{k}+\mathbf{G}},\,a_{\mathbf{k}-\mathbf{G}},\,d_{\mathbf{k}-\mathbf{G}})^{T}$ are four-spinors; $\alpha_{\mathbf{k}}=N^{-1/2}\sum_{\mathbf{r}}e^{i\mathbf{k}\cdot\mathbf{r}}\alpha_{\mathbf{r}}$ ($\alpha=a,\,b,\,c,\,d$), with $N$ the number of unit cells; and

[TABLE]

We have defined the functions $s_{\mu}(\mathbf{k})\equiv\sum_{j=1}^{3}e^{i(\mathbf{k}+\mu\mathbf{G})\cdot\bm{\delta}_{j}}$ ($\mu=-1,\,0,\,1$), and for simplicity, ignored all states near the GBZ $\Gamma$ point, which are $\sim 10\,{\rm eV}$ higher than the low-lying states, thus reducing the Hamiltonian’s dimensions from $12\times 12$ to $8\times 8$. The $\chi$ subspace consists of interlayer dimers at energies $\pm\gamma_{1}$, which we may project out to obtain a $4\times 4$ low-energy modelMcCann and Fal’ko (2006). As explained in the Supplement, we do so up to second order in perturbation theory by performing a standard Löwdin partitioningLöwdin (1951), and obtain the effective Bloch Hamiltonian about the $\Gamma$ point

[TABLE]

We have defined $v=a\sqrt{3}\gamma_{3}/2$ and $w^{2}=3a^{2}\gamma_{0}^{2}/4\gamma_{1}$, and kept terms proportional to $\gamma_{3}$ up to first order in momentum. $\sigma_{m}$ and $\tau_{n}$ are Pauli matrices acting upon the sublattice and valley subspaces, respectively. Eq. (8) is a good approximation to the full model (6), as long as the Kekulé distortion is small ($\Delta_{0}\lesssim 0.2$; see Fig. 2c and Supplement).

Figs. 2c–e show the main findings of this Letter: considering a small in-plane bond modulation of $\Delta_{0}\sim 0.1$, the low-energy dispersion of the $p6m$ Kek-Y graphene bilayer consists of a hexad of low-eccentricity elliptical Dirac bands. These bands can be divided into two sets, given by the Dirac points appearing at $\eta\mathbf{k}_{\rm D}=\eta k_{\rm D}\hat{\bm{\mathbf{x}}}$, with valley index $\eta=\pm 1$, and their corresponding symmetry partners obtained by $C_{3}$ rotations. We provide an analytic expression for $k_{\rm D}$, estimated from the effective model (8), in the Supplement. Here, we only point out that $k_{\rm D}$ is much smaller than the KBZ size ($k_{\rm D}\ll|\mathbf{G}|/\sqrt{3}$).

The origin of the multi-valley structure can be traced to the Lifshitz transition occurring at either $K$ point in the band structure of Bernal-stacked bilayer grapheneMcCann et al. (2007), where four Dirac cones appear below energies of order $1\,{\rm meV}$ (not shown; see, e.g., Ref. Varlet et al., 2015). Folding these cones into the KBZ produces the hexagonal node structure, with a doubly-degenerate cone exactly at the KBZ center, shown in gray in the top inset of Fig. 2c. It is clear, however, that the internode distance does not correspond to $k_{\rm D}$, which is instead determined by the Kekulé distortion strength, as described in the Supplement. More importantly, Fig. 2c shows that the Dirac cone bandwidth in the $p6m$ Kekulé bilayer is about 30 times larger than the Lifsitz transition energy. This means that multi-valley physics should survive at room temperature, and for carrier densities as large as $7\times 10^{11}{\rm cm}^{-2}$, or possibly higher ifGutiérrez et al. (2016) $\Delta_{0}>0.08$. By contrast, ultra-clean suspended or encapsulated samples are required to observe the Lifshitz transition in bilayer grapheneMartin et al. (2010); Varlet et al. (2014). An important conclusion may be drawn by this connection; namely, that massless Dirac fermions appear only for Bernal-stacked Kekulé bilayers, where the Lifshitz transition is present, and not for $AA$ stacking, which we also classify in the Supplement.

Linearizing $H_{\rm eff}(\mathbf{k})$ about $\eta\mathbf{k}_{\rm D}$, and projecting out the higher-energy parabolic bands up to first order in momentum, we obtain

[TABLE]

with expressions for the Fermi velocites $V_{x},\,V_{y}$ given in the Supplement. For $\Delta_{0}=0.08$, we find $V_{x}=1.93\times 10^{-3}c$ and $V_{y}=1.23\times 10^{-3}c$, where $c$ is the speed of light in vacuum, corresponding to elliptical Dirac bands of eccentricity $0.77$. Similarly to the case of graphene, Eq. (9) gives opposite chiralities $\eta=\pm 1$ for states about the Dirac points located at $\eta C_{3}^{n}\mathbf{k}_{\rm D}$ ($n=0,1,2$). This is verified numerically by the opposite Berry curvatures of adjacent valleys, evaluated using the full model (6), and shown in Fig. 2d.

By contrast to graphene, however, the bilayer nature of the present system allows opening a relativistic mass gap when an out-of-plane electric field is applied (Fig. 2e), producing a potential drop $\Delta V$ across the structure. This lowers the symmetry from $p6m$ to $p3$ by distinguishing between second-nearest interlayer neighbors, whose 2D projections would have otherwise been connected by $C_{6}$ rotations. For $|e\Delta V|\lesssim 50\,{\rm meV}$, the parameters of the effective model (8) become slightly renormalized, and a new term,

[TABLE]

is introduced, where $e$ is the (positive) electronic charge, and $v_{0,4}=a\gamma_{0,4}\sqrt{3}/2$. Similar effects (gap opening and symmetry lowering) are produced by intralayer sublattice asymmetries. By contrast, both the $p6m$ symmetry and the gapless Dirac spectrum are preserved by an on-site energy asymmetry between dimerizing- and non-dimerizing sites, caused by different molecular environmentsMucha-Kruczyński et al. (2010).

So far, we have focused on the most symmetric configuration that can result from stacking the two Kekulé-distorted graphene layers. A different arrangement, described by the space group $cmm$, is shown in Fig. 3a, and its corresponding low-energy dispersion in Fig. 3b. By performing numerical tight-binding calculations, we found that this system exhibits 4 inequivalent Dirac cones distributed preferentially along the $\hat{\bm{\mathbf{y}}}$ axis, mirroring the formation of chains of distorted interlayer bonds in its structure (see Fig. 3a). Other possible patterns are shown in the Supplement, including a chiral structure described by the $p6$ wallpaper group, for which we anticipate interesting optical activity. An in-depth exploration of the symmetry and electronic properties of this large family of materials will be published elsewhereRuiz-Tijerina et al. .

Multi-valley systems host vastly complex phenomena, both from the single- and many-particle perspectives. Here, we have proposed a family of Bernal-stacked bilayer graphene Kekulé superlattices that display Dirac bands centered at multiple valleys near the mini Brillouin zone center. In the most symmetric case, where the Kekulé distortion centers of the two layers coincide, we have predicted a system of multi-valley Dirac fermions with electrostatically tunable mass gaps, with great potential for many-body and quantum Hall physics at room temperature. The small momentum-space separation ($\sim 10^{-2}\,\mbox{\AA}^{-1}$) between Dirac valleys in all cases studied makes the band structure susceptible to long-range disorder, requiring samples with high crystalline quality to explore their multi-valley spectra. If Coulomb interactions are considered, the Fermi velocities $V_{x,y}$ reported above give an effective fine structure constant of $\alpha_{\rm eff}=e^{2}/\hbar V_{x,y}\approx 4.6$, corresponding to strong couplingSemenoff (2012). This makes the different Kekulé graphene bilayers studied in this Letter good candidates to observe many-body effects, including the anomalous quantum Hall insulatingNandkishore and Levitov (2010) and gapless symmetry-brokenLemonik et al. (2010); Vafek and Yang (2010) phases predicted for bilayer graphene. We believe that systems like the one proposed in this Letter may just be within the reach of current experimental capabilities, drawing from state-of-the art techniques currently used for manual stacking of van der Waals crystals. We anticipate that uncertainties in the exact interlayer alignment, inherent in the type of systems we propose, will lead to moiré pattern formation, adding a further layer of complexity to the already rich phenomenology predicted in this Letter.

D.R.-T. acknowledges funding from UNAM-DGAPA. E.A. thanks CONACyT for a doctoral scholarship. R.C.-B. thanks 20va Convocatoria Interna (UABC). E.A. and G.N. thank UNAM-DGAPA-PAPIIT project IN-102717 for finantial support. F.M. thanks the support of PAPIIT-UNAM, through project IN111317. D.R.-T. thanks M. Asmar for useful discussions during preparation of this paper.