Valley-momentum locking in a graphene superlattice with Y-shaped Kekulé bond texture

The coupling of orbital and spin degrees of freedom is a promising new direction in nano-electronics, referred to as 'spin-orbitronics', that aims at non-magnetic control of information carried by charge-neutral spin currents [1–3]. Graphene offers a rich platform for this research [4, 5], because the conduction electrons have three distinct spin quantum numbers: in addition to the spin magnetic moment s = ±1/2, there is the sublattice pseudospin σ = A, B and the valley isospin τ = K, K'. While the coupling of the electron spin s to its momentum p is a relativistic effect, and very weak in graphene, the coupling of σ to p is so strong that one has a pseudospin-momentum locking: the pseudospin points in the direction of motion, as a result of the helicity operator ${\boldsymbol{p}}\cdot {\boldsymbol{\sigma }}\equiv {p}_{x}{\sigma }_{x}+{p}_{y}{\sigma }_{y}$ in the Dirac Hamiltonian of graphene.

The purpose of this paper is to propose a way to obtain a similar handle on the valley isospin, by adding a term ${\boldsymbol{p}}\cdot {\boldsymbol{\tau }}$ to the Dirac Hamiltonian, which commutes with the pseudospin helicity and locks the valley to the direction of motion. We find that this valley-momentum locking should appear in a superlattice that has recently been realized experimentally by Gutiérrez et al [6, 7]: a superlattice of graphene grown epitaxially onto Cu(111), with the copper atoms in registry with the carbon atoms. One of six carbon atoms in each superlattice unit cell ($\sqrt{3}\times \sqrt{3}$ larger than the original graphene unit cell) have no copper atoms below them and acquire a shorter nearest-neighbor bond. The resulting Y-shaped periodic alternation of weak and strong bonds (see figure 1) is called a Kekulé-Y (Kek-Y) ordering, with reference to the Kekulé dimerization in a benzene ring (called Kek-O in this context) [7].

Zoom In Zoom Out Reset image size

The Kek-O and KeK-Y superlattices have the same Brillouin zone, with the K and K' valleys of graphene folded on top of each other. The Kek-O ordering couples the valleys by opening a gap in the Dirac cone [8–12], and it was assumed by Gutiérrez et al that the same applies to the Kek-Y ordering [6, 7]. While it is certainly possible that the graphene layer in the experiment is gapped by the epitaxial substrate (for example, by a sublattice-symmetry breaking ionic potential [13–15]), we find that the Y-shaped Kekulé bond ordering by itself does not impose a mass on the Dirac fermions⁴ . Instead, the valley degeneracy is broken by the helicity operator ${\boldsymbol{p}}\cdot {\boldsymbol{\tau }}$, which preserves the gapless Dirac point while locking the valley degree of freedom to the momentum. In a magnetic field the valley-momentum locking splits all Landau levels except for the zeroth Landau level, which remains pinned to zero energy.

2.1. Real-space formulation

A monolayer of carbon atoms has the tight-binding Hamiltonian

$$\begin{eqnarray}&&H=-\displaystyle \sum _{{\boldsymbol{r}}}\displaystyle \sum _{{\ell }=1}^{3}{t}_{{\boldsymbol{r}},{\ell }}\,{a}_{{\boldsymbol{r}}}^{\dagger }{b}_{{\boldsymbol{r}}+{{\boldsymbol{s}}}_{{\ell }}}^{}+{\rm{H}}.{\rm{c}}.,\end{eqnarray} \tag{ 1 }$$

describing the hopping with amplitude ${t}_{{\boldsymbol{r}},{\ell }}$ between an atom at site ${\boldsymbol{r}}=n{{\boldsymbol{a}}}_{1}+m{{\boldsymbol{a}}}_{2}$ ($n,m\in {\mathbb{Z}}$) on the A sublattice (annihilation operator ${a}_{{\boldsymbol{r}}}$) and each of its three nearest neighbors at ${\boldsymbol{r}}+{{\boldsymbol{s}}}_{{\ell }}$ on the B sublattice (annihilation operator ${b}_{{\boldsymbol{r}}+{{\boldsymbol{s}}}_{{\ell }}}$). The lattice vectors are defined by ${{\boldsymbol{s}}}_{1}=\tfrac{1}{2}(\sqrt{3},-1)$, ${{\boldsymbol{s}}}_{2}=-\tfrac{1}{2}(\sqrt{3},1)$, ${{\boldsymbol{s}}}_{3}=(0,1)$, ${{\boldsymbol{a}}}_{1}={{\boldsymbol{s}}}_{3}-{{\boldsymbol{s}}}_{1}$, ${{\boldsymbol{a}}}_{2}={{\boldsymbol{s}}}_{3}-{{\boldsymbol{s}}}_{2}$. All lengths are measured in units of the unperturbed C–C bond length a₀ ≡ 1.

For the uniform lattice, with ${t}_{{\boldsymbol{r}},{\ell }}\equiv {t}_{0}$, the band structure is given by [16]

$$\begin{eqnarray}&&E({\boldsymbol{k}})=\pm | \varepsilon ({\boldsymbol{k}})| ,\,\,\varepsilon ({\boldsymbol{k}})={t}_{0}\displaystyle \sum _{{\ell }=1}^{3}{{\rm{e}}}^{{\rm{i}}{\boldsymbol{k}}\cdot {{\boldsymbol{s}}}_{{\ell }}}.\end{eqnarray} \tag{ 2 }$$

There is a conical singularity at the Dirac points ${{\boldsymbol{K}}}_{\pm }=\tfrac{2}{9}\pi \sqrt{3}(\pm 1,\sqrt{3})$, where $E({{\boldsymbol{K}}}_{\pm })=0$. For later use we note the identities

$$\begin{eqnarray}&&\varepsilon ({\boldsymbol{k}})=\varepsilon ({\boldsymbol{k}}+3{{\boldsymbol{K}}}_{\pm })={{\rm{e}}}^{2\pi {\rm{i}}/3}\varepsilon ({\boldsymbol{k}}+{{\boldsymbol{K}}}_{+}+{{\boldsymbol{K}}}_{-}).\end{eqnarray} \tag{ 3 }$$

The bond-density wave that describes the Kek-O and Kek-Y textures has the form

$$\begin{eqnarray}{t}_{{\boldsymbol{r}},{\ell }}/{t}_{0} & = & 1+2\,\mathrm{Re}\,[{\rm{\Delta }}{{\rm{e}}}^{{\rm{i}}(p{{\boldsymbol{K}}}_{+}+q{{\boldsymbol{K}}}_{-})\cdot {{\boldsymbol{s}}}_{{\ell }}+{\rm{i}}{\boldsymbol{G}}\cdot {\boldsymbol{r}}}]\\ & = & 1+2{{\rm{\Delta }}}_{0}\cos \left[\phi +\tfrac{2}{3}\pi (m-n+{N}_{{\ell }})\right],\end{eqnarray} \tag{ 4a }$$

$$\begin{eqnarray}&&{N}_{1}=-q,\,\,{N}_{2}=-p,\,\,{N}_{3}\,=\,p+q,\,\,p,q\in {{\mathbb{Z}}}_{3}.\end{eqnarray} \tag{ 4b }$$

The Kekulé wave vector

$$\begin{eqnarray}&&{\boldsymbol{G}}\equiv {{\boldsymbol{K}}}_{+}-{{\boldsymbol{K}}}_{-}=\tfrac{4}{9}\pi \sqrt{3}(1,0)\end{eqnarray} \tag{ 5 }$$

couples the Dirac points. The coupling amplitude ${\rm{\Delta }}={{\rm{\Delta }}}_{0}{{\rm{e}}}^{{\rm{i}}\phi }$ may be complex, but the hopping amplitudes ${t}_{{\boldsymbol{r}},{\ell }}$ are real in order to preserve time-reversal symmetry. (We note that our definition of Δ differs by a factor 3 from that of [8].)

As illustrated in figure 1, the index

$$\begin{eqnarray}&&\nu =1+q-p\mathrm{mod}3\end{eqnarray} \tag{ 6 }$$

distinguishes the Kek-O texture (ν = 0) from the Kek-Y texture (ν =  ±1). Each Kekulé superlattice has a 2π/3 rotational symmetry, reduced from the 2π/6 symmetry of the graphene lattice. The two ν = ±1 Kek-Y textures are each others mirror image⁵ .

2.2. Transformation to momentum space

The Kek-O and Kek-Y superlattices have the same hexagonal Brillouin zone, with reciprocal lattice vectors ${{\boldsymbol{K}}}_{\pm }$—smaller by a factor $1/\sqrt{3}$ and rotated over 30° with respect to the original Brillouin zone of graphene (see figure 1). The Dirac points of unperturbed graphene are folded from the corner to the center of the Brillouin zone and coupled by the bond-density wave.

To study the coupling we Fourier transform the tight-binding Hamiltonian (1),

$$\begin{eqnarray}H({\boldsymbol{k}}) & = & -\varepsilon ({\boldsymbol{k}}){a}_{{\boldsymbol{k}}}^{\dagger }{b}_{{\boldsymbol{k}}}^{}-{\rm{\Delta }}\varepsilon ({\boldsymbol{k}}+p{{\boldsymbol{K}}}_{+}+q{{\boldsymbol{K}}}_{-}){a}_{{\boldsymbol{k}}+{\boldsymbol{G}}}^{\dagger }{b}_{{\boldsymbol{k}}}^{}\\ & & -{{\rm{\Delta }}}^{* }\varepsilon ({\boldsymbol{k}}-p{{\boldsymbol{K}}}_{+}-q{{\boldsymbol{K}}}_{-}){a}_{{\boldsymbol{k}}-{\boldsymbol{G}}}^{\dagger }{b}_{{\boldsymbol{k}}}^{}+{\rm{H}}.{\rm{c}}.\end{eqnarray} \tag{ 7 }$$

The momentum ${\boldsymbol{k}}$ still varies over the original Brillouin zone. In order to restrict it to the superlattice Brillouin zone we collect the annihilation operators at ${\boldsymbol{k}}$ and ${\boldsymbol{k}}\pm {\boldsymbol{G}}$ in the column vector ${c}_{{\boldsymbol{k}}}=({a}_{{\boldsymbol{k}}},{a}_{{\boldsymbol{k}}-{\boldsymbol{G}}},{a}_{{\boldsymbol{k}}+{\boldsymbol{G}}},{b}_{{\boldsymbol{k}}},{b}_{{\boldsymbol{k}}-{\boldsymbol{G}}},{b}_{{\boldsymbol{k}}+{\boldsymbol{G}}})$ and write the Hamiltonian in a 6 × 6 matrix form:

$$\begin{eqnarray}H({\boldsymbol{k}})=-{c}_{{\boldsymbol{k}}}^{\dagger }\left(\begin{array}{cc}0 & {{ \mathcal E }}_{\nu }({\boldsymbol{k}})\\ {{ \mathcal E }}_{\nu }^{\dagger }({\boldsymbol{k}}) & 0\end{array}\right){c}_{{\boldsymbol{k}}},\end{eqnarray} \tag{ 8a }$$

$$\begin{eqnarray}{{ \mathcal E }}_{\nu }=\left(\begin{array}{ccc}{\varepsilon }_{0} & \tilde{{\rm{\Delta }}}{\varepsilon }_{\nu +1} & {\tilde{{\rm{\Delta }}}}^{* }{\varepsilon }_{-\nu -1}\\ {\tilde{{\rm{\Delta }}}}^{* }{\varepsilon }_{1-\nu } & {\varepsilon }_{-1} & \tilde{{\rm{\Delta }}}{\varepsilon }_{\nu }\\ \tilde{{\rm{\Delta }}}{\varepsilon }_{\nu -1} & {\tilde{{\rm{\Delta }}}}^{* }{\varepsilon }_{-\nu } & {\varepsilon }_{1}\end{array}\right),\end{eqnarray} \tag{ 8b }$$

$$\begin{eqnarray}&&\tilde{{\rm{\Delta }}}={{\rm{e}}}^{2\pi {\rm{i}}(p+q)/3}{\rm{\Delta }},\,\,{\varepsilon }_{n}=\varepsilon ({\boldsymbol{k}}+n{\boldsymbol{G}}),\end{eqnarray} \tag{ 8c }$$

where we used equation (3).

3.1. Gapless spectrum

The low-energy spectrum is governed by the four modes ${u}_{{\boldsymbol{k}}}=({a}_{{\boldsymbol{k}}-{\boldsymbol{G}}},{a}_{{\boldsymbol{k}}+{\boldsymbol{G}}},{b}_{{\boldsymbol{k}}-{\boldsymbol{G}}},{b}_{{\boldsymbol{k}}+{\boldsymbol{G}}})$, which for small ${\boldsymbol{k}}$ lie near the Dirac points at $\pm {\boldsymbol{G}}$. (We identify the K valley with $+{\boldsymbol{G}}$ and the K' valley with $-{\boldsymbol{G}}$.) Projection onto this subspace reduces the six-band Hamiltonian (8) to an effective four-band Hamiltonian,

$$\begin{eqnarray}{H}_{\mathrm{eff}}=-{u}_{{\boldsymbol{k}}}^{\dagger }\left(\begin{array}{cc}0 & {h}_{\nu }\\ {h}_{\nu }^{\dagger } & 0\end{array}\right){u}_{{\boldsymbol{k}}},\,\,{h}_{\nu }=\left(\begin{array}{cc}{\varepsilon }_{-1} & \tilde{{\rm{\Delta }}}{\varepsilon }_{\nu }\\ {\tilde{{\rm{\Delta }}}}^{* }{\varepsilon }_{-\nu } & {\varepsilon }_{1}\end{array}\right).\end{eqnarray} \tag{ 9 }$$

Corrections to the low-energy spectrum from virtual transitions to the higher bands are of order ${{\rm{\Delta }}}_{0}^{2}$. We will include these corrections later, but for now assume Δ₀ ≪ 1 and neglect them.

The ${\boldsymbol{k}}$-dependence of ${\varepsilon }_{n}$ may be linearized near ${\boldsymbol{k}}=0$,

$$\begin{eqnarray}&&{\varepsilon }_{0}=3{t}_{0},\,\,{\varepsilon }_{\pm 1}={\hslash }{v}_{0}(\mp {k}_{x}+{{\rm{i}}{k}}_{y})+\mathrm{order}\,({k}^{2}),\end{eqnarray} \tag{ 10 }$$

with Fermi velocity ${v}_{0}=\tfrac{3}{2}{t}_{0}{a}_{0}/{\hslash }$. The corresponding 4-component Dirac equation has the form

$$\begin{eqnarray}&&{ \mathcal H }\left(\begin{array}{c}{{\rm{\Psi }}}_{K^{\prime} }\\ {{\rm{\Psi }}}_{K}\end{array}\right)=E\left(\begin{array}{c}{{\rm{\Psi }}}_{K^{\prime} }\\ {{\rm{\Psi }}}_{K}\end{array}\right),\,\,{ \mathcal H }=\left(\begin{array}{cc}{v}_{0}{\boldsymbol{p}}\cdot {\boldsymbol{\sigma }} & \tilde{{\rm{\Delta }}}{Q}_{\nu }\\ {\tilde{{\rm{\Delta }}}}^{* }{Q}_{\nu }^{\dagger } & {v}_{0}{\boldsymbol{p}}\cdot {\boldsymbol{\sigma }}\end{array}\right),\end{eqnarray} \tag{ 11a }$$

$$\begin{eqnarray}&&{{\rm{\Psi }}}_{K^{\prime} }=\left(\begin{array}{c}-{\psi }_{B,K^{\prime} }\\ {\psi }_{A,K^{\prime} }\end{array}\right),\,\,{{\rm{\Psi }}}_{K}=\left(\begin{array}{c}{\psi }_{A,K}\\ {\psi }_{B,K}\end{array}\right),\end{eqnarray} \tag{ 11b }$$

$$\begin{eqnarray}{Q}_{\nu }=\left(\begin{array}{cc}{\varepsilon }_{-\nu }^{* } & 0\\ 0 & -{\varepsilon }_{\nu }\end{array}\right)=\left\{\begin{array}{lc}3{t}_{0}{\sigma }_{z} & \mathrm{if}\,\,\nu =0,\\ {v}_{0}(\nu {p}_{x}-{{\rm{i}}{p}}_{y}){\sigma }_{0} & \mathrm{if}\,\,| \nu | =1.\end{array}\right.\end{eqnarray} \tag{ 11c }$$

The spinor Ψ_K contains the wave amplitudes on the A and B sublattices in valley K and similarly ${{\rm{\Psi }}}_{K^{\prime} }$ for valley K', but note the different ordering of the components [17]⁶ . We have defined the momentum operator ${\boldsymbol{p}}=-{\rm{i}}{\hslash }\partial /\partial {\boldsymbol{r}}$, with ${\boldsymbol{p}}\cdot {\boldsymbol{\sigma }}={p}_{x}{\sigma }_{x}+{p}_{y}{\sigma }_{y}$. The Pauli matrices ${\sigma }_{x},{\sigma }_{y},{\sigma }_{z}$, with σ₀ the unit matrix, act on the sublattice degree of freedom.

For the Kek-O texture we recover the gapped spectrum of Kekulé dimerized graphene [8],

$$\begin{eqnarray}&&{E}^{2}={v}_{0}^{2}| {\boldsymbol{p}}{| }^{2}+{(3{t}_{0}{{\rm{\Delta }}}_{0})}^{2}\,\,\mathrm{for}\,\,\nu =0.\end{eqnarray} \tag{ 12 }$$

The Kek-Y texture, instead, has a gapless spectrum,

$$\begin{eqnarray}&&{E}_{\pm }^{2}={v}_{0}^{2}{(1\pm {{\rm{\Delta }}}_{0})}^{2}| {\boldsymbol{p}}{| }^{2},\,\,\mathrm{for}\,\,| \nu | =1,\end{eqnarray} \tag{ 13 }$$

consisting of a pair of linearly dispersing modes with different velocities ${v}_{0}(1\pm {{\rm{\Delta }}}_{0})$. The two qualitatively different dispersions are contrasted in figure 2.

Zoom In Zoom Out Reset image size

3.2. Valley-momentum locking

The two gapless modes in the Kek-Y superlattice are helical, with both the sublattice pseudospin and the valley isospin locked to the direction of motion. To see this, we consider the ν = 1 Kek-Y texture with a real $\tilde{{\rm{\Delta }}}={{\rm{\Delta }}}_{0}$. (Complex $\tilde{{\rm{\Delta }}}$ and ν = − 1 are equivalent upon a unitary transformation.) The Dirac Hamiltonian (11) can be written in the compact form

$$\begin{eqnarray}&&{ \mathcal H }={v}_{\sigma }({\boldsymbol{p}}\cdot {\boldsymbol{\sigma }})\otimes {\tau }_{0}+{v}_{\tau }\,{\sigma }_{0}\otimes ({\boldsymbol{p}}\cdot {\boldsymbol{\tau }}),\end{eqnarray} \tag{ 14 }$$

with the help of a second set of Pauli matrices τ_x, τ_y, τ_z and unit matrix τ₀ acting on the valley degree of freedom⁷ . The two velocities are defined by v_σ = v₀ and v_τ = v₀Δ₀.

An eigenstate of the current operator

$$\begin{eqnarray}&&{j}_{\alpha }=\partial { \mathcal H }/\partial {p}_{\alpha }={v}_{\sigma }\,{\sigma }_{\alpha }\otimes {\tau }_{0}+{v}_{\tau }\,{\sigma }_{0}\otimes {\tau }_{\alpha }\end{eqnarray} \tag{ 15 }$$

with eigenvalue ${v}_{\sigma }\pm {v}_{\tau }$ is an eigenstate of σ_α with eigenvalue +1 and an eigenstate of τ_α with eigenvalue ±1. (The two Pauli matrices act on different degrees of freedom, so they commute and can be diagonalized independently.) This valley-momentum locking does not violate time-reversal symmetry, since the time-reversal operation in the superlattice inverts all three vectors ${\boldsymbol{p}}$, ${\boldsymbol{\sigma }}$, and ${\boldsymbol{\tau }}$, and hence leaves ${ \mathcal H }$ unaffected⁸ :

$$\begin{eqnarray}&&({\sigma }_{y}\otimes {\tau }_{y}){{ \mathcal H }}^{* }({\sigma }_{y}\otimes {\tau }_{y})={ \mathcal H }.\end{eqnarray} \tag{ 16 }$$

The valley-momentum locking does break the sublattice symmetry, since ${ \mathcal H }$ no longer anticommutes with σ_z, but another chiral symmetry involving both sublattice and valley degrees of freedom remains:

$$\begin{eqnarray}&&({\sigma }_{z}\otimes {\tau }_{z}){ \mathcal H }=-{ \mathcal H }({\sigma }_{z}\otimes {\tau }_{z}).\end{eqnarray} \tag{ 17 }$$

3.3. Landau level quantization

A perpendicular magnetic field B in the z-direction (vector potential ${\boldsymbol{A}}$ in the x–y plane), breaks the time-reversal symmetry (16) via the substitution ${\boldsymbol{p}}\mapsto -{\rm{i}}{\hslash }\partial /\partial {\boldsymbol{r}}+e{\boldsymbol{A}}({\boldsymbol{r}})\equiv {\boldsymbol{\Pi }}$. The chiral symmetry (17) is preserved, so the Landau levels are still symmetrically arranged around E = 0, as in unperturbed graphene. Because the two helicity operators ${\boldsymbol{\Pi }}\cdot {\boldsymbol{\sigma }}$ and ${\boldsymbol{\Pi }}\cdot {\boldsymbol{\tau }}$ do not commute for ${\boldsymbol{A}}\ne 0$, they can no longer be diagonalized independently. In particular, this means the Landau level spectrum is not simply a superposition of two spectra of Dirac fermions with different velocities.

It is still possible to calculate the spectrum analytically (see appendix A). We find Landau levels at energies ${E}_{n}^{+},{E}_{n}^{-},-{E}_{n}^{+},-{E}_{n}^{-}$, $n=0,1,2,\,\ldots$, given by

$$\begin{eqnarray}&&{E}_{n}^{\pm }={E}_{B}{[2n+1\pm \sqrt{1+n(n+1){(4{v}_{\sigma }{v}_{\tau })}^{2}{\bar{v}}^{-4}}]}^{1/2},\end{eqnarray} \tag{ 18 }$$

with the definitions $\bar{v}=\sqrt{{v}_{\sigma }^{2}+{v}_{\tau }^{2}}$ and ${E}_{B}=\bar{v}\sqrt{{\hslash }{eB}}$.

In unperturbed graphene all Landau levels have a twofold valley degeneracy⁹ : ${E}_{n}^{+}={E}_{n+1}^{-}$ for v_τ = 0. This includes the zeroth Landau level: ${E}_{0}^{-}=0=-{E}_{0}^{-}$. A nonzero v_τ breaks the valley degeneracy of all Landau levels at $E\ne 0$, but a valley-degenerate zero-mode ${E}_{0}^{-}=0$ remains, see figure 3.

Zoom In Zoom Out Reset image size

The absence of a splitting of the zeroth Landau level can be understood as a topological protection in the context of an index theorem [18–21], which requires that either ${{\rm{\Pi }}}_{+}\equiv {{\rm{\Pi }}}_{x}+{\rm{i}}{{\rm{\Pi }}}_{y}$ or ${{\rm{\Pi }}}_{-}\equiv {{\rm{\Pi }}}_{x}-{\rm{i}}{{\rm{\Pi }}}_{y}$ has a zero-mode. If we decompose ${ \mathcal H }={{\rm{\Pi }}}_{+}{S}_{-}+{{\rm{\Pi }}}_{-}{S}_{+}$, with ${S}_{\pm }={v}_{\sigma }({\sigma }_{x}\pm {\rm{i}}{\sigma }_{y})+{v}_{\tau }({\tau }_{x}\pm {\rm{i}}{\tau }_{y})$, we see that both S₊ and S₋ have a rank-two null space¹⁰ , spanned by the spinors ${\psi }_{\pm }^{(1)}$ and ${\psi }_{\pm }^{(2)}$. So if ${{\rm{\Pi }}}_{\pm }{f}_{\pm }=0$, a twofold degenerate zero-mode of ${ \mathcal H }$ is formed by the states ${f}_{\pm }{\psi }_{\mp }^{(1)}$ and ${f}_{\pm }{\psi }_{\mp }^{(2)}$.

All of this is distinctive for the Kek-Y bond order: for the Kek-O texture it is the other way around—the Landau levels have a twofold valley degeneracy except for the nondegenerate Landau level at the edge of the band gap¹¹ .

So far we have assumed ${{\rm{\Delta }}}_{0}\ll 1$, and one might ask how robust our findings are to finite-Δ₀ corrections, involving virtual transitions from the ${\varepsilon }_{\pm 1}$ bands near E = 0 to the ${\varepsilon }_{0}$ band near E = 3t₀. We have been able to include these to all orders in Δ₀ (see appendix B), and find that the entire effect is a renormalization of the velocities v_σ and v_τ in the Hamiltonian (14), which retains its form as a sum of two helicity operators. For real Δ = Δ₀ the renormalization is given by v_σ = v₀ρ₊, v_τ = v₀ρ₋ with

$$\begin{eqnarray}&&{\rho }_{\pm }=\tfrac{1}{2}(1-{{\rm{\Delta }}}_{0})\left(\displaystyle \frac{1+2{{\rm{\Delta }}}_{0}}{\sqrt{1+2{{\rm{\Delta }}}_{0}^{2}}}\pm 1\right).\end{eqnarray} \tag{ 19 }$$

For complex ${\rm{\Delta }}={{\rm{\Delta }}}_{0}{{\rm{e}}}^{{\rm{i}}\phi }$ the nonlinear renormalization introduces a dependence on the phase ϕ modulo 2π/3.

What this renormalization shows is that, as expected for a topological protection, the robustness of the zeroth Landau level to the Kek-Y texture is not limited to perturbation theory—also strong modulations of the bond strength cannot split it away from E = 0.

In zero magnetic field the low-energy Hamiltonian (14) does not couple the pseudospin σ and valley τ degrees of freedom. A ${\boldsymbol{\sigma }}\otimes {\boldsymbol{\tau }}$ coupling is introduced in the Kek-Y superlattice by an ionic potential μ_Y on the carbon atoms that line up with the carbon vacancies—the atoms located at each center of a red Y in figure 1. We consider this effect for the ν = 1 Kek-Y texture with a real $\tilde{{\rm{\Delta }}}={{\rm{\Delta }}}_{0}$.

The ionic potential acts on one-third of the A sublattice sites, labeled ${{\boldsymbol{r}}}_{{\rm{Y}}}$. (For ν = −1 it would act on one-third of the B sublattice sites.) Fourier transformation of the on-site contribution ${\mu }_{{\rm{Y}}}{\sum }_{{{\boldsymbol{r}}}_{{\rm{Y}}}}{a}_{{{\boldsymbol{r}}}_{{\rm{Y}}}}^{\dagger }{a}_{{{\boldsymbol{r}}}_{{\rm{Y}}}}$ to the tight-binding Hamiltonian (1) gives with the help of the lattice sum

$$\begin{eqnarray}&&\displaystyle \sum _{{{\boldsymbol{r}}}_{{\rm{Y}}}}{{\rm{e}}}^{{\rm{i}}{\boldsymbol{k}}\cdot {{\boldsymbol{r}}}_{{\rm{Y}}}}\propto \delta ({\boldsymbol{k}})+\delta ({\boldsymbol{k}}-{\boldsymbol{G}})+\delta ({\boldsymbol{k}}+{\boldsymbol{G}})\end{eqnarray} \tag{ 20 }$$

the momentum-space Hamiltonian

$$\begin{eqnarray}H({\boldsymbol{k}})=-{c}_{{\boldsymbol{k}}}^{\dagger }\left(\begin{array}{cc}{M}_{{\rm{Y}}} & {{ \mathcal E }}_{1}({\boldsymbol{k}})\\ {{ \mathcal E }}_{1}^{\dagger }({\boldsymbol{k}}) & 0\end{array}\right){c}_{{\boldsymbol{k}}},\end{eqnarray} \tag{ 21a }$$

$$\begin{eqnarray}{M}_{{\rm{Y}}}=-{\mu }_{{\rm{Y}}}\left(\begin{array}{ccc}1 & 1 & 1\\ 1 & 1 & 1\\ 1 & 1 & 1\end{array}\right).\end{eqnarray} \tag{ 21b }$$

The ${{ \mathcal E }}_{1}$ block is still given by equation (8). The additional ${M}_{{\rm{Y}}}$-block breaks the chiral symmetry.

Projection onto the subspace spanned by ${u}_{{\boldsymbol{k}}}=({a}_{{\boldsymbol{k}}-{\boldsymbol{G}}},{a}_{{\boldsymbol{k}}+{\boldsymbol{G}}},{b}_{{\boldsymbol{k}}-{\boldsymbol{G}}},{b}_{{\boldsymbol{k}}+{\boldsymbol{G}}})$ gives the effective Hamiltonian

$$\begin{eqnarray}{H}_{\mathrm{eff}}=-{u}_{{\boldsymbol{k}}}^{\dagger }\left(\begin{array}{cc}{m}_{{\rm{Y}}} & {h}_{1}\\ {h}_{1}^{\dagger } & 0\end{array}\right){u}_{{\boldsymbol{k}}},\,\,{m}_{{\rm{Y}}}=-{\mu }_{{\rm{Y}}}\left(\begin{array}{cc}1 & 1\\ 1 & 1\end{array}\right).\end{eqnarray} \tag{ 22 }$$

The corresponding Dirac Hamiltonian has the form (11) with an additional ${\boldsymbol{\sigma }}\otimes {\boldsymbol{\tau }}$ coupling,

$$\begin{eqnarray}{ \mathcal H } & = & {v}_{\sigma }({\boldsymbol{p}}\cdot {\boldsymbol{\sigma }})\otimes {\tau }_{0}+{v}_{\tau }\,{\sigma }_{0}\otimes ({\boldsymbol{p}}\cdot {\boldsymbol{\tau }})+\tfrac{1}{2}{\mu }_{{\rm{Y}}}\\ & & +\tfrac{1}{2}{\mu }_{{\rm{Y}}}({\sigma }_{x}\otimes {\tau }_{x}+{\sigma }_{y}\otimes {\tau }_{y}-{\sigma }_{z}\otimes {\tau }_{z}).\end{eqnarray} \tag{ 23 }$$

The energy spectrum,

$$\begin{eqnarray}{E}_{\pm }^{(1)} & = & \pm ({v}_{\sigma }-{v}_{\tau })| {\boldsymbol{p}}| ,\\ {E}_{\pm }^{(2)} & = & {\mu }_{{\rm{Y}}}\pm \sqrt{{({v}_{\sigma }+{v}_{\tau })}^{2}| {\boldsymbol{p}}{| }^{2}+{\mu }_{{\rm{Y}}}^{2}},\end{eqnarray} \tag{ 24 }$$

has two bands that cross linearly in p at E = 0, while the other two bands have a quadratic p-dependence (see figure 4). The pseudospin and valley isospin orientation for the two bands is illustrated in figure 5.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The three bands ${E}_{+}^{(1)}$, ${E}_{-}^{(1)}$, ${E}_{-}^{(2)}$ that intersect at p = 0 are reminiscent of a spin-one Dirac one. Such a dispersion is a known feature of a potential modulation that involves only one-third of the atoms on one sublattice [14, 15]. The spectrum remains gapless even though the chiral symmetry is broken. This is in contrast to the usual staggered potential between A and B sublattices, which opens a gap via a ${\sigma }_{z}\otimes {\tau }_{z}$ term [16].

In summary, we have shown that the Y-shaped Kekulé bond texture (Kek-Y superlattice) in graphene preserves the massless character of the Dirac fermions. This is fundamentally different from the gapped band structure resulting from the original Kekulé dimerization [8–11] (Kek-O superlattice), and contrary to expectations from its experimental realization [6, 7].

The gapless low-energy Hamiltonian ${ \mathcal H }={v}_{\sigma }{\boldsymbol{p}}\cdot {\boldsymbol{\sigma }}+{v}_{\tau }{\boldsymbol{p}}\cdot {\boldsymbol{\tau }}$ is the sum of two helicity operators, with the momentum ${\boldsymbol{p}}$ coupled independently to both the sublattice pseudospin ${\boldsymbol{\sigma }}$ and the valley isospin ${\boldsymbol{\tau }}$. This valley-momentum locking is distinct from the coupling of the valley to a pseudo-magnetic field that has been explored as an enabler for valleytronics [23], and offers a way for a momentum-controlled valley precession. The broken valley degeneracy would also remove a major obstacle for spin qubits in graphene [24].

A key experimental test of our theoretical predictions would be a confirmation that the Kek-Y superlattice has a gapless spectrum, in stark contrast to the gapped Kek-O spectrum. In the experiment by Gutiérrez et al on a graphene/Cu heterostructure the Kek-Y superlattice is formed by copper vacancies that are in registry with one out of six carbon atoms [6, 7]. These introduce the Y-shaped hopping modulations shown in figure 1, but in addition will modify the ionic potential felt by the carbon atom at the center of the Y. Unlike the usual staggered potential between A and B sublattices, this potential modulation in an enlarged unit cell does not open a gap [14, 15]. We have also checked that the Dirac cone remains gapless if we include hoppings beyond nearest neighbor. All of this gives confidence that the gapless spectrum will survive in a realistic situation.

Further research in other directions could involve the Landau level spectrum, to search for the unique feature of a broken valley degeneracy coexisting with a valley-degenerate zero-mode. The graphene analogs in optics and acoustics [25] could also provide an interesting platform for a Kek-Y superlattice with a much stronger amplitude modulation than can be realized with electrons.

We have benefited from discussions with A Akhmerov, V Cheianov, J Hutasoit, P Silvestrov, and D Varjas. This research was supported by the Netherlands Organization for Scientific Research (NWO/OCW) and an ERC Synergy Grant.

We calculate the spectrum in a perpendicular magnetic field of a graphene sheet with a Kekulé-Y bond texture. We start by rewriting the Hamiltonian (14), with ${\boldsymbol{\Pi }}={\boldsymbol{p}}+e{\boldsymbol{A}}$, in the form

$$\begin{eqnarray}&&{ \mathcal H }=\tfrac{1}{2}{{\rm{\Pi }}}_{+}{S}_{-}+\tfrac{1}{2}{{\rm{\Pi }}}_{-}{S}_{+}+\mu {\sigma }_{z}\otimes {\tau }_{z},\end{eqnarray} \tag{ A1 }$$

in terms of the raising and lowering operators

$$\begin{eqnarray}{{\rm{\Pi }}}_{\pm } & = & {{\rm{\Pi }}}_{x}\pm {\rm{i}}{{\rm{\Pi }}}_{y},\,\,{\sigma }_{\pm }={\sigma }_{x}\pm {\rm{i}}{\sigma }_{y},\,\,{\tau }_{\pm }={\tau }_{x}\pm {\rm{i}}{\tau }_{y},\\ {S}_{\pm } & = & {v}_{\sigma }\,{\sigma }_{\pm }\otimes {\tau }_{0}+{v}_{\tau }\,{\sigma }_{0}\otimes {\tau }_{\pm }.\end{eqnarray} \tag{ A2 }$$

The chiral-symmetry breaking term μσ_z ⨂ τ_z that we have added will serve a purpose later on.

We know that the Hermitian operator Ω = Π₊Π₋ has eigenvalues ${\omega }_{n}=2n{\hslash }{eB}$, n = 0, 1, 2, ..., in view of the commutator $[{{\rm{\Pi }}}_{-},{{\rm{\Pi }}}_{+}]=2{\hslash }{eB}$. So the strategy is to express the secular equation $\det (E-{ \mathcal H })=0$ in a form that involves only the mixed products ${{\rm{\Pi }}}_{+}{{\rm{\Pi }}}_{-}$, and no ${{\rm{\Pi }}}_{+}^{2}$ or ${{\rm{\Pi }}}_{-}^{2}$. This is achieved by means of a unitary transformation, as follows.

We define the unitary matrix

$$\begin{eqnarray}&&U=\exp \left[\tfrac{1}{4}{\rm{i}}\pi ({\sigma }_{0}+{\sigma }_{z})\otimes {\tau }_{y}\right]\end{eqnarray} \tag{ A3 }$$

and reduce the determinant of a 4 × 4 matrix to that of a 2 × 2 matrix:

$$\begin{eqnarray}\det ({ \mathcal H }-E) & = & \det {U}^{\dagger }({ \mathcal H }-E)U\\ & = & \det \left(\begin{array}{cc}-E+\mu & {R}^{\dagger }\\ R & -E-\mu \end{array}\right)\\ & = & \left\{\begin{array}{cc}\det ({E}^{2}-{\mu }^{2}-{{RR}}^{\dagger }) & \mathrm{if}\,\,E\ne \mu ,\\ \det ({E}^{2}-{\mu }^{2}-{R}^{\dagger }R) & \mathrm{if}\,\,E\ne -\mu ,\end{array}\right.\end{eqnarray} \tag{ A4 }$$

$$\begin{eqnarray}\mathrm{with}\,\,R=\left(\begin{array}{cc}-{v}_{\tau }{{\rm{\Pi }}}_{-} & {v}_{\sigma }{{\rm{\Pi }}}_{-}\\ -{v}_{\sigma }{{\rm{\Pi }}}_{+} & {v}_{\tau }{{\rm{\Pi }}}_{+}\end{array}\right).\end{eqnarray} \tag{ A5 }$$

The matrix product RR^† is not of the desired form, but R^†R is,

$$\begin{eqnarray}{R}^{\dagger }R=\left(\begin{array}{cc}{v}_{\sigma }^{2}{{\rm{\Pi }}}_{-}{{\rm{\Pi }}}_{+}+{v}_{\tau }^{2}{{\rm{\Pi }}}_{+}{{\rm{\Pi }}}_{-} & -{v}_{\sigma }{v}_{\tau }({{\rm{\Pi }}}_{-}{{\rm{\Pi }}}_{+}+{{\rm{\Pi }}}_{+}{{\rm{\Pi }}}_{-})\\ -{v}_{\sigma }{v}_{\tau }({{\rm{\Pi }}}_{-}{{\rm{\Pi }}}_{+}+{{\rm{\Pi }}}_{+}{{\rm{\Pi }}}_{-}) & {v}_{\sigma }^{2}{{\rm{\Pi }}}_{+}{{\rm{\Pi }}}_{-}+{v}_{\tau }^{2}{{\rm{\Pi }}}_{-}{{\rm{\Pi }}}_{+}\end{array}\right),\end{eqnarray} \tag{ A6 }$$

involving only ${{\rm{\Pi }}}_{+}{{\rm{\Pi }}}_{-}={\rm{\Omega }}$ and ${{\rm{\Pi }}}_{-}{{\rm{\Pi }}}_{+}={\rm{\Omega }}+{\omega }_{1}$. Hence the determinant is readily evaluated for $E\ne -\mu$,

$$\begin{eqnarray}\det ({ \mathcal H }-E)=\det ({E}^{2}-{\mu }^{2}-{R}^{\dagger }R)=\displaystyle \prod _{n=0}^{\infty }\det \left(\begin{array}{cc}{E}^{2}-{\mu }^{2}-{\bar{v}}^{2}{\omega }_{n}-{v}_{\sigma }^{2}{\omega }_{1} & {v}_{\sigma }{v}_{\tau }(2{\omega }_{n}+{\omega }_{1})\\ {v}_{\sigma }{v}_{\tau }(2{\omega }_{n}+{\omega }_{1}) & {E}^{2}-{\mu }^{2}-{\bar{v}}^{2}{\omega }_{n}-{v}_{\tau }^{2}{\omega }_{1}\end{array}\right),\end{eqnarray} \tag{ A7 }$$

where we have abbreviated $\bar{v}=\sqrt{{v}_{\sigma }^{2}+{v}_{\tau }^{2}}$.

Equating the determinant to zero and solving for E we find four sets of energy eigenvalues ${E}_{n}^{+},{E}_{n}^{-},-{E}_{n}^{+},-{E}_{n}^{-}$, given by

$$\begin{eqnarray}&&{({E}_{n}^{\pm })}^{2}-{\mu }^{2}=\left({\omega }_{n}+\tfrac{1}{2}{\omega }_{1}\right){\bar{v}}^{2}\pm \tfrac{1}{2}\sqrt{{\omega }_{1}^{2}{\bar{v}}^{4}+{(4{v}_{\sigma }{v}_{\tau })}^{2}{\omega }_{n}{\omega }_{n+1}}={E}_{B}^{2}[2n+1\pm \sqrt{1+n(n+1){(4{v}_{\sigma }{v}_{\tau })}^{2}{\bar{v}}^{-4}}].\end{eqnarray} \tag{ A8 }$$

In the second equation we introduced the energy scale ${E}_{B}={\hslash }\bar{v}/{l}_{m}$, with ${l}_{m}=\sqrt{{\hslash }/{eB}}$ the magnetic length. The B-independent level ${E}_{0}^{-}=\mu$ becomes a zero-mode in the limit $\mu \to 0$.

As a check on the calculation, we note that for μ = 0, v_τ = 0 we recover the valley-degenerate Landau level spectrum of graphene [16],

$$\begin{eqnarray}&&{E}_{n}^{-}=({\hslash }{v}_{\sigma }/{l}_{m})\sqrt{2n},\,\,{E}_{n}^{+}={E}_{n+1}^{-}.\end{eqnarray} \tag{ A9 }$$

Another special case of interest is μ = 0, ${v}_{\sigma }={v}_{\tau }\equiv {v}_{0}$, when the two modes of Dirac fermions have velocities ${v}_{\sigma }\pm {v}_{\tau }$ equal to 0 and 2v₀. From equation (A8) we find the Landau level spectrum

$$\begin{eqnarray}&&{E}_{n}^{-}=0,\,\,{E}_{n}^{+}=2({\hslash }{v}_{0}/{l}_{m})\sqrt{2n+1}.\end{eqnarray} \tag{ A10 }$$

The mode with zero velocity remains B-independent, while the mode with velocity 2v₀ produces a sequence of Landau levels with a 1/2 offset in the n-dependence.

We seek to reduce the six-band Hamiltonian (8) to an effective 4 × 4 Hamiltonian that describes the low-energy spectrum near ${\boldsymbol{k}}=0$. For ${{\rm{\Delta }}}_{0}\ll 1$ we can simply project onto the 2 × 2 lower-right subblock of ${{ \mathcal E }}_{\nu }$, which for the $| \nu | =1$ Kek-Y bond modulation vanishes linearly in ${\boldsymbol{k}}$. This subblock is coupled to the ${\varepsilon }_{0}$ band near $E=3{t}_{0}$ by matrix elements of order Δ₀, so virtual transitions to this higher band contribute to the low-energy spectrum in order ${{\rm{\Delta }}}_{0}^{2}$. We will now show how to include these effects to all order in Δ₀.

One complication when we go beyond the small-Δ₀ regime is that the phase ϕ of the modulation amplitude can no longer be removed by a unitary transformation. As we will see, the low-energy Hamiltonian depends on ϕ modulo 2π/3—so we do not need to distinguish between the phase of $\tilde{{\rm{\Delta }}}={{\rm{e}}}^{2\pi {\rm{i}}(p+q)/3}{\rm{\Delta }}$ and the phase of Δ. The choice between $\nu =\pm 1$ still does not matter, the two Kek-Y modulations being related by a mirror symmetry. For definiteness we take ν = +1.

We define the unitary matrix

$$\begin{eqnarray}V=\left(\begin{array}{cc}{\rm{\Phi }} & 0\\ 0 & {\rm{\Phi }}\end{array}\right)\left(\begin{array}{cc}{ \mathcal V } & 0\\ 0 & {\mathbb{1}}\end{array}\right),\,\,{\rm{\Phi }}=\left(\begin{array}{ccc}1 & 0 & 0\\ 0 & {{\rm{e}}}^{-{\rm{i}}\phi } & 0\\ 0 & 0 & {{\rm{e}}}^{{\rm{i}}\phi }\end{array}\right),\end{eqnarray} \tag{ B1a }$$

$$\begin{eqnarray}{ \mathcal V }=\displaystyle \frac{1}{2{D}_{0}}\left(\begin{array}{ccc}2 & -2{{\rm{\Delta }}}_{0} & -2{{\rm{\Delta }}}_{0}\\ 2{{\rm{\Delta }}}_{0} & 1+{D}_{0} & 1-{D}_{0}\\ 2{{\rm{\Delta }}}_{0} & 1-{D}_{0} & 1+{D}_{0}\end{array}\right),\end{eqnarray} \tag{ B1b }$$

with ${D}_{0}=\sqrt{1+2{{\rm{\Delta }}}_{0}^{2}}$ and evaluate

$$\begin{eqnarray}{V}^{\dagger }\left(\begin{array}{cc}0 & {{ \mathcal E }}_{1}\\ {{ \mathcal E }}_{1}^{\dagger } & 0\end{array}\right)V=\left(\begin{array}{cc}0 & {\tilde{{ \mathcal E }}}_{1}\\ {\tilde{{ \mathcal E }}}_{1}^{\dagger } & 0\end{array}\right),\end{eqnarray} \tag{ B2a }$$

$$\begin{eqnarray}{\tilde{{ \mathcal E }}}_{1}={{ \mathcal V }}^{\dagger }{{ \mathcal E }}_{1}=\left(\begin{array}{ccc}{D}_{0}{\varepsilon }_{0} & {\rho }_{0}^{* }{\varepsilon }_{-1} & {\rho }_{0}{\varepsilon }_{1}\\ 0 & {\rho }_{+}{\varepsilon }_{-1} & {\rho }_{-}^{* }{\varepsilon }_{1}\\ 0 & {\rho }_{-}{\varepsilon }_{-1} & {\rho }_{+}^{* }{\varepsilon }_{1}\end{array}\right),\end{eqnarray} \tag{ B2b }$$

$$\begin{eqnarray}&&{\rho }_{\pm }=\displaystyle \frac{1}{2{D}_{0}}[1-2{{\rm{\Delta }}}_{0}^{2}\pm {D}_{0}+{{\rm{e}}}^{-3{\rm{i}}\phi }{{\rm{\Delta }}}_{0}(1\mp {D}_{0})],\end{eqnarray} \tag{ B2c }$$

$$\begin{eqnarray}&&{\rho }_{0}=\displaystyle \frac{{{\rm{\Delta }}}_{0}}{{D}_{0}}(2+{{\rm{e}}}^{3{\rm{i}}\phi }{{\rm{\Delta }}}_{0}).\end{eqnarray} \tag{ B2d }$$

The matrix elements that couple the lower-right 2 × 2 subblock of ${\tilde{{ \mathcal E }}}_{1}$ to ${\varepsilon }_{0}$ are now of order k, so the effect on the low-energy spectrum is of order k² and can be neglected—to all orders in Δ₀.

The resulting effective low-energy Hamiltonian has the 4 × 4 form (9), with h₁ replaced by

$$\begin{eqnarray}{h}_{1}=\left(\begin{array}{cc}{\rho }_{+}{\varepsilon }_{-1} & {\rho }_{-}^{* }{\varepsilon }_{1}\\ {\rho }_{-}{\varepsilon }_{-1} & {\rho }_{+}^{* }{\varepsilon }_{1}\end{array}\right).\end{eqnarray} \tag{ B3 }$$

The phases of ${\rho }_{\pm }=| {\rho }_{\pm }| {{\rm{e}}}^{{\rm{i}}{\theta }_{\pm }}$ can be eliminated by one more unitary transformation, with the 4 × 4 diagonal matrix

$$\begin{eqnarray}&&{\rm{\Theta }}=\mathrm{diag}({{\rm{e}}}^{{\rm{i}}{\theta }_{-}},{{\rm{e}}}^{{\rm{i}}{\theta }_{+}},{{\rm{e}}}^{{\rm{i}}{\theta }_{+}+{\rm{i}}{\theta }_{-}},1),\end{eqnarray} \tag{ B4 }$$

which results in

$$\begin{eqnarray}{{\rm{\Theta }}}^{\dagger }\left(\begin{array}{cc}0 & {h}_{1}\\ {\tilde{h}}_{1}^{\dagger } & 0\end{array}\right){\rm{\Theta }}=\left(\begin{array}{cc}0 & {\tilde{h}}_{1}\\ {\tilde{h}}_{1}^{\dagger } & 0\end{array}\right),\,\,{\tilde{h}}_{1}=\left(\begin{array}{cc}| {\rho }_{+}| {\varepsilon }_{-1} & | {\rho }_{-}| {\varepsilon }_{1}\\ | {\rho }_{-}| {\varepsilon }_{-1} & | {\rho }_{+}| {\varepsilon }_{1}\end{array}\right).\end{eqnarray} \tag{ B5 }$$

Finally, we arrive at the effective Hamiltonian (14), with renormalized velocities:

$$\begin{eqnarray}&&{ \mathcal H }={v}_{\sigma }({\boldsymbol{p}}\cdot {\boldsymbol{\sigma }})\otimes {\tau }_{0}+{v}_{\tau }\,{\sigma }_{0}\otimes ({\boldsymbol{p}}\cdot {\boldsymbol{\tau }}),\,\,{v}_{\sigma }=| {\rho }_{+}| {v}_{0},\,\,{v}_{\tau }=| {\rho }_{-}| {v}_{0},\end{eqnarray} \tag{ B6 }$$

$$\begin{eqnarray}&&| {\rho }_{\pm }{| }^{2}=\displaystyle \frac{1}{2{D}_{0}^{2}}(1+3{{\rm{\Delta }}}_{0}^{4}\pm {D}_{0}(1-3{{\rm{\Delta }}}_{0}^{2})+2{{\rm{\Delta }}}_{0}^{3}(\pm {D}_{0}-2)\cos 3\phi ).\end{eqnarray} \tag{ B7 }$$

To third order in Δ₀ we have

$$\begin{eqnarray}&&{v}_{\sigma }/{v}_{0}=1-\tfrac{3}{2}{{\rm{\Delta }}}_{0}^{2}-\tfrac{1}{2}{{\rm{\Delta }}}_{0}^{3}\cos 3\phi ,\,\,{v}_{\tau }/{v}_{0}={{\rm{\Delta }}}_{0}-\tfrac{3}{2}{{\rm{\Delta }}}_{0}^{2}\cos 3\phi +\tfrac{1}{16}{{\rm{\Delta }}}_{0}^{3}(1-9\cos 6\phi )+{ \mathcal O }({{\rm{\Delta }}}_{0}^{4}).\end{eqnarray} \tag{ B8 }$$

For real Δ, when ϕ = 0 and ρ_± is real, equation (B7) simplifies to

$$\begin{eqnarray}&&{\rho }_{\pm }=\tfrac{1}{2}(1-{{\rm{\Delta }}}_{0})\left(\displaystyle \frac{1+2{{\rm{\Delta }}}_{0}}{\sqrt{1+2{{\rm{\Delta }}}_{0}^{2}}}\pm 1\right).\end{eqnarray} \tag{ B9 }$$

The velocities of the two Dirac modes are then given by

$$\begin{eqnarray}{v}_{1} & = & {v}_{\sigma }+{v}_{\tau }={v}_{0}\displaystyle \frac{(1-{{\rm{\Delta }}}_{0})(1+2{{\rm{\Delta }}}_{0})}{\sqrt{1+2{{\rm{\Delta }}}_{0}^{2}}},\\ {v}_{2} & = & {v}_{\sigma }-{v}_{\tau }={v}_{0}(1-{{\rm{\Delta }}}_{0}).\end{eqnarray} \tag{ B10 }$$

More generally, for complex ${\rm{\Delta }}={{\rm{\Delta }}}_{0}{{\rm{e}}}^{{\rm{i}}\phi }$ both v₁ and v₂ become ϕ-dependent to second order in Δ₀, see figure B1.

Zoom In Zoom Out Reset image size

Note that the asymmetry in $\pm {{\rm{\Delta }}}_{0}$ vanishes for ϕ = π/6. For this phase the superlattice has three different bond strengths (see figure B2) that are symmetrically arranged around the unperturbed value t₀.

Zoom In Zoom Out Reset image size