Valley-momentum locking in a graphene superlattice with Y-shaped Kekul\'e bond texture

Recent experiments by Guti\'errez $\textit{et al.}$ [Nature Phys. $\textbf{12}$, 950 (2016)] on a graphene-copper superlattice have revealed an unusual Kekul\'e bond texture in the honeycomb lattice --- a Y-shaped modulation of weak and strong bonds with a wave vector connecting two Dirac points. We show that this so-called"Kek-Y"texture produces two species of massless Dirac fermions, with valley isospin locked parallel or antiparallel to the direction of motion. In a magnetic field $B$ the valley degeneracy of the $B$-dependent Landau levels is removed by the valley-momentum locking --- but a $B$-independent and valley-degenerate zero-mode remains.


I. INTRODUCTION
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][2][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 p · σ ≡ p x σ x + p y σ 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 p · τ 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 ( √ 3 × √ 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 Fig. 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].
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][9][10][11][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][14][15]), we find that the Y-shaped Kekulé bond ordering by itself does not impose a mass on the Dirac fermions [16]. Instead, the valley degeneracy is broken by the helicity operator p · τ , 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.

II. TIGHT-BINDING MODEL
A. Real-space formulation A monolayer of carbon atoms has the tight-binding Hamiltonian describing the hopping with amplitude t r, between an atom at site r = na 1 + ma 2 (n, m ∈ Z) on the A sublattice (annihilation operator a r ) and each of its three nearest neighbors at r + s on the B sublattice (annihilation operator b r+s ). The lattice vectors are defined by s 1 = 1 2 ( All lengths are measured in units of the unperturbed C-C bond length a 0 ≡ 1.
For the uniform lattice, with t r, ≡ t 0 , the band structure is given by [17] There is a conical singularity at the Dirac points For later use we note the identities The bond-density wave that describes the Kek-O and Kek-Y textures has the form t r, /t 0 = 1 + 2 Re ∆e i(pK++qK−)·s +iG·r (4a) The Kekulé wave vector couples the Dirac points. The coupling amplitude ∆ = ∆ 0 e iφ may be complex, but the hopping amplitudes t r, are real in order to preserve time-reversal symmetry. As illustrated in Fig. 1, the index 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 [18].

B. Transformation to momentum space
The Kek-O and Kek-Y superlattices have the same hexagonal Brillouin zone, with reciprocal lattice vectors K ± -smaller by a factor 1/ √ 3 and rotated over 30 • with respect to the original Brillouin zone of graphene (see Fig. 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 tightbinding Hamilonian (1), The momentum k still varies over the original Brillouin zone. In order to restrict it to the superlattice Brillouin zone we collect the annihilation operators at k and k ± G in the column vector c k = (a k , a k−G , a k+G , b k , b k−G , b k+G ) and write the Hamiltonian in a 6 × 6 matrix form: where we used Eq. (3).

A. Gapless spectrum
The low-energy spectrum is governed by the four modes u k = (a k−G , a k+G , b k−G , b k+G ), which for small k lie near the Dirac points at ±G. (We identify the K valley with +G and the K valley with −G.) Projection onto this subspace reduces the six-band Hamiltonian (8) to an effective four-band Hamiltonian, Corrections to the low-energy spectrum from virtual transitions to the higher bands are of order ∆ 2 0 . We will include these corrections later, but for now assume ∆ 0 1 and neglect them. The k-dependence of ε n may be linearized near k = 0,

4-component Dirac equation has the form
The spinor Ψ K contains the wave amplitudes on the A and B sublattices in valley K and similarly Ψ K for valley K , but note the different ordering of the components [19]. We have defined the momentum operator p = −i ∂/∂r, with p · σ = p x σ x + p y σ y . The Pauli matrices σ x , σ y , σ z , with σ 0 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], The Kek-Y texture, instead, has a gapless spectrum, consisting of a pair of linearly dispersing modes with different velocities v 0 (1 ± ∆ 0 ). The two qualitatively different dispersions are contrasted in Fig. 2.

B. 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 ∆ = ∆ 0 . (Complex∆ and ν = −1 are equivalent upon a unitary transformation.) The Dirac Hamiltonian (11) can be written in the compact form with the help of a second set of Pauli matrices τ x , τ y , τ z and unit matrix τ 0 acting on the valley degree of freedom. The two velocities are defined by An eigenstate of the current operator with eigenvalue v σ ± v τ 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 p, σ, and τ , and hence leaves H unaffected [20]: The valley-momentum locking does break the sublattice symmetry, since H no longer anticommutes with σ z , but another chiral symmetry involving both sublattice and valley degrees of freedom remains:

C. Landau level quantization
A perpendicular magnetic field B in the z-direction (vector potential A in the x-y plane), breaks the time-reversal symmetry (16) via the substitution p → −i ∂/∂r + eA(r) ≡ Π. 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 Π · σ and Π · τ do not commute for A = 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 App. A). We find Landau levels at energies E + n , E − n , −E + n , −E − n , n = 0, 1, 2, . . ., given by  The absence of a splitting of the zeroth-Landau level can be understood as a topological protection in the context of an index theorem [22][23][24][25], which requires that , we see that both S + and S − have a rank-two null space [26], spanned by the spinors ψ

IV. EFFECT OF VIRTUAL TRANSITIONS TO HIGHER BANDS
So far we have assumed ∆ 0 1, and one might ask how robust our findings are to finite-∆ 0 corrections, involving virtual transitions from the ε ±1 bands near E = 0 to the ε 0 band near E = 3t 0 . We have been able to include these to all orders in ∆ 0 (see App. 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 ∆ = ∆ 0 the renormalization is given by For complex ∆ = ∆ 0 e iφ 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.

V. PSEUDOSPIN-VALLEY COUPLING
In zero magnetic field the low-energy Hamiltonian (14) does not couple the pseudospin σ and valley τ degrees of freedom. A σ ⊗τ 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 Fig. 1. We consider this effect for the ν = 1 Kek-Y texture with a real∆ = ∆ 0 .
The ionic potential acts on one-third of the A sublattice sites, labeled r Y . (For ν = −1 it would act on one-third of the B sublattice sites.) Fourier transformation of the on-site contribution µ Y rY a † rY a rY to the tight-binding Hamiltonian (1) gives with the help of the lattice sum the momentum-space Hamiltonian The E 1 block is still given by Eq. (8). The additional M Y -block breaks the chiral symmetry. Projection onto the subspace spanned by u k = (a k−G , a k+G , b k−G , b k+G ) gives the effective Hamiltonian The corresponding Dirac Hamiltonian has the form (11) with an additional σ ⊗ τ coupling, The energy spectrum, has two bands that cross linearly in p at E = 0, while the other two bands have a quadratic p-dependence. (See Fig. 4.)

The three bands E
(1) − 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 σ z ⊗ τ z term [17].

VI. DISCUSSION
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][9][10][11] (Kek-O superlattice), and contrary to expectations from its experimental realization [6,7].
The gapless low-energy Hamiltonian H = v σ p · σ + v τ p · τ is the sum of two helicity operators, with the momentum p coupled independently to both the sublattice pseudospin σ and the valley isospin τ . This valleymomentum locking is distinct from the coupling of the valley to a pseudo-magnetic field that has been explored as an enabler for valleytronics [29], 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 [30].
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 Fig. 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 neigbor. 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 valleydegenerate zero-mode. The graphene analogues in optics and acoustics [31] could also provide an interesting platform for a Kek-Y superlattice with a much stronger amplitude modulation than can be realized with electrons. [27] In a Kek-O superlattice the Landau levels are given by E 2 n = (3t0∆0) 2 + 2n eBv 2 0 , n = 0, 1, 2, . . ., with a twofold valley degeneracy for n ≥ 1 and a nondegenerate zeroth Landau level at ±3t0∆0.
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 Π = p + eA, in the form in terms of the raising and lowering operators 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 ω n = 2n eB, n = 0, 1, 2, . . ., in view of the commutator [Π − , Π + ] = 2 eB. So the strategy is to express the secular equation det(E − H) = 0 in a form that involves only the mixed products Π + Π − , and no Π 2 + or Π 2 − . This is achieved by means of a unitary transformation, as follows.
We define the unitary matrix and reduce the determinant of a 4 × 4 matrix to that of a 2 × 2 matrix: The matrix product RR † is not of the desired form, but R † R is, involving only Π + Π − = Ω and Π − Π + = Ω + ω 1 . Hence the determinant is readily evaluated for E = −µ, where we have abbreviatedv = v 2 σ + v 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 In the second equation we introduced the energy scale E B = v/l m , with l m = /eB the magnetic length. The B-independent level E − 0 = µ becomes a zero-mode in the limit µ → 0.
As a check on the calculation, we note that for µ = 0, v τ = 0 we recover the valley-degenerate Landau level spectrum of graphene [17], Another special case of interest is µ = 0, v σ = v τ ≡ v 0 , when the two modes of Dirac fermions have velocities v σ ± v τ equal to 0 and 2v 0 . From Eq. (A8) we find the Landau level spectrum The mode with zero velocity remains B-independent, while the mode with velocity 2v 0 produces a sequence of Landau levels with a 1/2 offset in the n-dependence.
Appendix B: Calculation of the low-energy Hamiltonian to all orders in the Kek-Y bond modulation We seek to reduce the six-band Hamiltonian (8) to an effective 4 × 4 Hamiltonian that describes the low-energy spectrum near k = 0. For ∆ 0 1 we can simply project onto the 2 × 2 lower-right subblock of E ν , which for the |ν| = 1 Kek-Y bond modulation vanishes linearly in k. This subblock is coupled to the ε 0 band near E = 3t 0 by matrix elements of order ∆ 0 , so virtual transitions to this higher band contribute to the low-energy spectrum in order ∆ 2 0 . We will now show how to include these effects to all order in ∆ 0 .
One complication when we go beyond the small-∆ 0 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 don't need to distinguish between the phase of∆ = e 2πi(p+q)/3 ∆ and the phase of ∆. The choice between ν = ±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 with D 0 = 1 + 2∆ 2 0 and evaluate The matrix elements that couple the lower-right 2 × 2 subblock ofẼ 1 to ε 0 are now of order k, so the effect on the low-energy spectrum is of order k 2 and can be neglected -to all orders in ∆ 0 .

(B5)
Finally, we arrive at the effective Hamiltonian (14), with renormalized velocities: To third order in ∆ 0 we have For real ∆, when φ = 0 and ρ ± is real, Eq. (B7) simplifies to The velocities of the two Dirac modes are then given by More generally, for complex ∆ = ∆ 0 e iφ both v 1 and v 2 become φ-dependent to second order in ∆ 0 , see Fig. 5.
Note that the asymmetry in ±∆ 0 vanishes for φ = π/6. For this phase the superlattice has three different bond strengths (see Fig. 6) that are symmetrically arranged around the unperturbed value t 0 .