Stability of zero energy Dirac touchings in the honeycomb Hofstadter problem

TL;DR

This paper investigates the stability of zero energy Dirac points in the band structure of electrons on a honeycomb lattice under magnetic flux, revealing symmetry conditions that protect these touchings.

Contribution

It demonstrates that Dirac touchings are protected by anti-unitary particle-hole symmetry and lattice symmetries, and analyzes their robustness to symmetry-lowering perturbations.

Findings

Dirac touchings are guaranteed by combined symmetry conditions.

Locations of Dirac points are symmetry-dependent.

Dirac points exhibit stability under certain perturbations.

Abstract

We study the band structure of electrons hopping on a honeycomb lattice with 1/q (q integer) flux quanta through each elementary hexagon. In the nearest neighbor hopping model the two bands that eventually form the n = 0 Landau level have 2q zero energy Dirac touchings. In this work, we study the conditions needed for these Dirac points and their stability to various perturbations. We prove that these touchings and their locations are guaranteed by a combination of an anti-unitary particle-hole symmetry and the lattice symmetries of the honeycomb structure. We also study the stability of the Dirac touchings to one-body perturbations that explicitly lower the symmetry.

Equations

$${\mathbf{a}}_1$$

$${\mathbf{a}}_2$$

$${\mathbf{r}}_{\mu }(\mathbf{n})$$

$$=$$

$$\frac{eB\sqrt{3}{a}^2}{2\hslash }=\frac{2\pi }{q}\equiv \chi$$

$$\mathbf{A}=-By\mathbf{x}$$

$${\varphi }_{\mu \nu }^L(\mathbf{n};\mathbf{\Delta n})$$

$$H_{\mathrm{nn}}$$

$$d_{\mu ,n_1,n_2}=e^{-i\frac{\chi }{4}{n}_2^2+i\frac{\chi }{6}(n_1-n_2)}{c}_{\mu }(n_1,n_2).$$

$${\varphi }_{\mu \nu }^{\mathrm{OG}}(\mathbf{n};\mathbf{\Delta n})$$

$$H_{\mathrm{nn}}$$

$${\mathbb{T}}_{{\mathbf{a}}_1}{c}_{\mu }(n_1,n_2){\mathbb{T}}_{{\mathbf{a}}_1}^{†}$$

$${\mathbb{T}}_{{\mathbf{a}}_2}{c}_{\mu }(n_1,n_2){\mathbb{T}}_{{\mathbf{a}}_2}^{†}$$

$${\mathbb{R}}_{2\pi /3}{c}_{\mu }(n_1,n_2){\mathbb{R}}_{2\pi /3}^{†}$$

$$=e^{i\chi \left(m_1{m}_2+\frac{m_2(m_2-1)}{2}\right)}{c}_{\mu }(m_1,m_2-\mu ),$$

$${\mathbb{R}}_{\pi }{c}_{\mu }(n_1,n_2){\mathbb{R}}_{\pi }^{†}=e^{-i\chi n_1}{c}_{1-\mu }(-n_1,-n_2)$$

$$\mathbb{S}{c}_A(n_1,n_2){\mathbb{S}}^{-1}=$$

$$\mathbb{S}{c}_B(n_1,n_2){\mathbb{S}}^{-1}=$$

$$\mathbb{S}i{\mathbb{S}}^{-1}=$$

$$\mathbf{K}_{n}=\frac{\pi}{q}\bigg{(}2n-q+\frac{1}{3}\bigg{)}{\hat{x}}-\frac{\pi}{q\sqrt{3}}(2n-q-1){\hat{y}}$$

$${\mathbf{k}}_{\gamma }={\mathbf{k}}_R+{\mathbf{G}}_2\frac{q+1}{2q}+\gamma \frac{{\mathbf{G}}_1}{q}$$

$$H(\mathbf{k})=\left(\begin{array}[]{cc}0&M\\ M^{\dagger}&0\end{array}\right)$$

$$f_{n}(E)=\sum\limits_{\alpha=0}^{q-1}\bigg{(}\psi^{(E)}_{A,\alpha,n}c_{A,\alpha,n}+\psi^{(E)}_{B,\alpha,n}c_{B,\alpha,n}\bigg{)}$$

$${\mathbb{R}}_{\frac{2\pi }{3}}{f}_n(E){\mathbb{R}}_{\frac{2\pi }{3}}^{†}=\sum _{n^{\prime }=0}^{q-1}|t_{n{n}^{\prime }}|e^{i{\varphi }_{n{n}^{\prime }}}{f}_{n^{\prime }}(E)$$

$$\langle\mu,\alpha,n|\mathbb{R}_{\frac{2\pi}{3}}|\nu,\beta,n\rangle=\left(\begin{array}[]{cc}R_{A}(n)&0\\ 0&R_{B}(n)\end{array}\right)$$

$$det\bigg{(}R_{A}(n)\big{(}R_{B}(n)\big{)}^{\dagger}\bigg{)}=real$$

$$q\bigg{(}R_{A}(n)\big{(}R_{B}(n)\big{)}^{\dagger}\bigg{)}_{\alpha\beta}=e^{i\frac{\chi}{3}+i\chi\beta}\delta_{\beta,\alpha-1}$$

$$arg\bigg{(}det\bigg{(}R_{A}(n)\big{(}R_{B}(n)\big{)}^{\dagger}\bigg{)}\bigg{)}=\frac{2\pi}{3}$$

$$\mathcal{Q}(H)=\frac{1}{2\pi i}{\int }_0^{2\pi }d\theta {\nabla }_{\theta }ln\det M(\theta )$$

$${\mathbb{T}}_{{\mathbf{a}}_1}{c}_{\mu }(n_1,n_2){\mathbb{T}}_{{\mathbf{a}}_1}^{†}$$

$${\mathbb{T}}_{{\mathbf{a}}_2}{c}_{\mu }(n_1,n_2){\mathbb{T}}_{{\mathbf{a}}_2}^{†}$$

$${\mathbb{R}}_{2\pi /3}{c}_{\mu }(n_1,n_2){\mathbb{R}}_{2\pi /3}^{†}$$

$$=e^{i\chi \left(m_1{m}_2+\frac{m_2(m_2-1)}{2}\right)}{c}_{\mu }(m_1,m_2-\mu ),$$

$${\mathbb{R}}_{\pi }{c}_{\mu }(n_1,n_2){\mathbb{R}}_{\pi }^{†}=e^{-i\chi n_1}{c}_{1-\mu }(-n_1,-n_2)$$

Full text

Honeycomb Lattice, Optimal Gauge (OG), Dirac Points, Chiral Symmetry ($\mathbb{S}$)

I Introduction

The band structure of electronic energy levels is a fascinating consequence of quantum mechanics applied to a solid. Ashcroft and Mermin (1976) The study of the topology of band structures has received tremendous attention in the last decade highlighted by the discovery of a variety of topological insulators and nodal semi-metals. Kane and Mele (2005a, b); Moore and Balents (2007); Fu et al. (2007); Fu and Kane (2007); Roy (2009); Hasan and Kane (2010) An important theme that has emerged is the importance of symmetries in protecting the distinction between insulating states and also the gaplessness of semi-metals. Mañes et al. (2007); Murakami (2007); Wan et al. (2011); Turner and Vishwanath (2013)

The study of the linear band touching in graphene Mañes et al. (2007); Hou and Chen (2015) has played a profound role in the unfolding of these discoveries. It is now well known that any tight-binding model of graphene (our discussion here will ignore both spin-orbit coupling and electron-electron interactions) with time-reversal symmetry and the symmetry of the honeycomb lattice has two independent Dirac touchings in its Brillouin zone at the K and K′ points. These touchings are stable to a number of quadratic perturbations. If the translational symmetry of the Bravais lattice is preserved, the touchings are stable to any perturbation that preserves rotation by $\pi$ around the honeycomb center and time reversal. Such perturbations can cause the Dirac touchings to move in the BZ, but they cannot gap them out. Breaking inversion or time-reversal symmetry individually lead to a trivial or Chern insulator respectively. A periodic perturbation that breaks the translational symmetry of the original Bravais lattice with a wavevector that connects the K and K′ points can also gap the Dirac points out (e.g. a Kekule dimerization). The derivation of these results is reviewed in Appendix A.

Our goal in this paper is to generalize these results to the honeycomb lattice in a magnetic field. This leads to integer quantum Hall states, historically the first examples of topological states of matterThouless et al. (1982), which in turn inspired the construction of the first lattice model of a Chern band. Haldane (1988)

In this work we study the $2q$ Dirac touchings that arise in the central two bands of the nearest-neighbor tight-binding honeycomb lattice when a flux of $1/q$ ($q$ integer) times the flux quantum is introduced into each elementary honeycomb plaquette – the so-called Hofstadter problem. Hofstadter (1976) There has been quite a bit of work on the Hofstadter problem on the honeycomb lattice. For the nearest-neighbor hopping model, previous work has, among other things, studied the spectrum and the eigenstates Rammal, R. (1985); Kohmoto and Sedrakyan (2006); Agazzi et al. (2014), the Diophantine equation and Chern number characterizing gapped statesSato et al. (2008), the crossover from Dirac-like behavior to conventional nonrelativistic behavior Hatsugai et al. (2006); ANDREI BERNEVIG et al. (2006), and the approach to the continuum limit $q\to\infty$. Kohmoto and Sedrakyan (2006); ANDREI BERNEVIG et al. (2006) The existence of $2q$ Dirac band touchings of the central two bands in the nearest-neighbor hopping model was noticed by several authors, and explored thoroughly more recently. Rhim and Park (2012) It was also pointed out that adding a next-nearest neighbor hopping gaps the Dirac points out. Karnaukhov (2018)

Here we extend the discussion in two ways: We first prove explicitly that certain specific symmetries protect the $2q$ Dirac touchings in a family of hopping models with arbitrary range hoppings, and $1/q$ flux. Second, we study the stability of these linear touchings to various one-body perturbations that lower the symmetry.

II Model

Throughout this paper we will be interested in the problem of spinless fermions hopping on the honeycomb lattice in the presence of a uniform magnetic field. We will study the problem in the tight binding limit and assume that each unit cell of the honeycomb lattice encloses a fraction $1/q$ of the flux quantum.

II.1 Gauge

We use the following conventions to define our honeycomb lattice. The two lattice vectors defining the primitive triangular lattice are,

[TABLE]

With these definitions the vectors describing the sites of the honeycomb lattice are,

[TABLE]

where ${\bf n}=(n_{1},n_{2})$ is a pair of integers and $\mu=0,1$ for the A and B sublattice respectively.

Once we introduce a rational magnetic field

[TABLE]

the hoppings acquire phases and we have to enlarge our unit cell in order to obtain two commuting translations, enabling us to apply Bloch’s theorem to compute the band structure. This enlarged unit cell is the magnetic unit cell (MUC). From now on we will use $a=1$. It will be useful for us to start with a continuum gauge, obtain the hopping phases, and then transform to the final gauge. We begin with the standard Landau gauge,

[TABLE]

We introduce the external magnetic field into the hopping model using the Peierls substitution to calculate the phase of the matrix elements. Using this gauge and the standard formula for the Peierls phase between two lattice points described by ${\bf n},\mu$ and ${\bf n+\Delta n},\nu$:

[TABLE]

where $\chi\equiv eBa^{2}\sqrt{3}/(2\hbar)=2\pi/q$ is the flux per unit cell of our system (in units of the flux quantum $\frac{h}{e}$), and ${\bf\Delta n}=(\Delta n_{1},\Delta n_{2})$. From the expression it is clear that $\phi^{L}_{\mu\nu}({\bf n};{\bf\Delta n})$ depends explicitly on $n_{2}$ and is $2\pi$ periodic only after $2q$ steps in the $n_{2}$ direction. Since $n_{1}$ does not appear it is periodic in every step of $n_{1}$. This means that we need to include $2q$ unit cells of the triangular lattice in our magnetic unit cell. We can see this explicitly by constructing the nearest neighbor hopping Hamiltonian in the Landau gauge,

[TABLE]

which clearly repeats itself with a magnetic unit cell consisting of $2q$ triangular unit cells. This is somewhat unsatisfactory since with a flux of $2\pi/q$ in a triangular unit, a gauge should exist in which there are only $q$ triangular units in the magnetic unit cells – the minimum size of unit cell needed to enclose an integer number of flux quanta. This can be resolved by working in the so-called optimal gauge (OG). To achieve this we make the following gauge transformation from the $d$ fermions (Landau gauge) to a set of $c$ fermions (OG),

[TABLE]

Using the transformation we can now compute the Peierls phase between two arbitrary sites on the honeycomb in the OG,

[TABLE]

From this formula, in the OG it is clear that the phases repeat themselves after $q$ steps in the $n_{2}$ direction, and thus Bloch’s theorem can be applied with only $q$ units of the triangular lattice in the magnetic unit cell (which contain 2$q$ lattice sites). We shall choose the magnetic unit cell shown in Fig. 1 in the rest of the paper and refer to it as the MUC. We can see the periodicity of the MUC explicitly by working out the nearest neighbor Hamiltonian in the OG,

[TABLE]

which clearly repeats itself after $q$ steps in the $n_{2}$ direction. The advantage of constructing the OG starting from the Landau gauge is that we now have a definite prescription to compute the Peierls’ phase for an arbitrary hopping matrix element in this gauge, Eq. (10). This allows us to write down hopping models with an arbitrary range of hopping such that all close paths enclose precisely the flux corresponding to a uniform external magnetic field, all the while still maintaining the MUC containing 2$q$ sites.

III Dirac touchings

Working in the OG we have computed the band structure for various ranges of tight binding models. This involves the diagonlization of a $2q\times 2q$ matrix for each ${\bf k}$ in the first Brillouin zone. The unit cell we have chosen and other lattice conventions are shown in Fig. 1.

Fig. 2 shows the electronic structure of the nearest neighbor model with $q=4$. Our focus in this paper is on the finite-$q$ electronic structure of the two central Bloch bands that eventually form the zero energy $n=0$ continuum Landau levels. In particular, as has been noticed in previous work, for the nearest-neighbor model the two bands have $2q$ linear band touchings that form a honeycomb lattice in reciprocal space. Rhim and Park (2012) As $q$ is increased the bandwidth of these bands decreases exponentially ANDREI BERNEVIG et al. (2006) and eventually as $q\rightarrow\infty$ we recover dispersionless Landau levels.

Are these Dirac touchings special to the nearest neighbor model or are they generic to the inclusion of further neighbor hoppings? It is known Karnaukhov (2019) that a next-nearest neighbor hopping gaps out the Dirac points. The formula Eq. (10) in the optimal gauge consistently allows us to include any range of hopping in the presence of a uniform field. We shall prove below that the Dirac points and their location are stable as long as the further neighbor hoppings are bipartite, i.e. they only connect sites on A with sites on B and maintain spatial the symmetries of the honeycomb lattice. If any A-A and B-B hoppings are included they gap out the Dirac fermions even if the honeycomb spatial symmetries are maintained.

As will be crucial for our discussion, the bipartite hopping structure has an extra symmetry that is broken when hopping between same sublattices is included. We note that in our problem both time-reversal symmetry $\mathbb{T}$ and the standard bipartite particle-hole symmetry $\mathbb{C}$ are broken, since physically they both reverse the direction of the external magnetic field. However, the product of the two, the sublattice symmetry $\mathbb{S}$ (which is an anti-unitary many-body particle-hole transformation) commutes with the Hamiltonian when only bipartite hoppings are included. The sublattice symmetry $\mathbb{S}$ is distinct from the “hidden symmetry” Hou (2013, 2014); Hou and Chen (2017, 2015) which exists on certain lattices and is also antiunitary, but in addition involves a translation and a sublattice exchange.

We prove explicitly that with the added constraint of the presence of $\mathbb{S}$ (in addition to all the lattice symmetries of honeycomb graphene lattice) all hopping models in the presence of a uniform magnetic field on the honeycomb lattice have $2q$ Dirac touchings at zero energy at the same locations in the BZ as the nearest neighbor model. We have tested this assertion by numerical diagonalization for a variety of different choices of the range and magnitude of the hoppings.

III.1 Proof of Dirac touching at special points

We will prove, by contradiction, that for generic Hamiltonians preserving the lattice and $\mathbb{S}$, there are necessarily zero-energy states at the 2$q$ special points in the BZ (the same points for which the nearest neighbor model has Dirac touchings).

Let us briefly introduce the action of various symmetries on the fermion operators in the OG (A more detailed discussion is presented in the appendices). Here $(n_{1},n_{2})$ are the two integers that label the location of a Bravais lattice sites of the original unit cell (not the magnetic unit cell), $\mu,\nu=0,1$ label the $A$ and $B$ sublattices. The two translations in the $\mathbf{a}_{1}$ and $\mathbf{a}_{2}$ directions act as follows:

[TABLE]

A rotation by $\frac{2\pi}{3}$ around the $A$ site at $n_{1}=n_{2}=0$, $\mathbb{R}_{\frac{2\pi}{3}}$, acts as follows:

[TABLE]

where $m_{1}=-n_{1}-n_{2}$ and $m_{2}=n_{1}$. A rotation by $\pi$ around the center of a vertical nearest-neighbor $AB$ bond, $\mathbb{R}_{\pi}$, acts as follows:

[TABLE]

Finally, the anti-unitary particle-hole symmetry $\mathbb{S}$ acts as follows:

[TABLE]

As has been noticed in previous work the nearest-neighbor only hopping model has Dirac touchings in the central two bands. We reproduce the locations, labelled by $n=0,\dots,q-1$, from Appendix D here for convenience. The $q$ $K$-type points, for $q$ odd, are

[TABLE]

with the $q$ $\mathbf{K}^{\prime}$-type points being $\mathbf{K}^{\prime}_{n}=-\mathbf{K}_{n}$.

Here are some properties of the $\mathbf{k}$-space Hamiltonians at these points that we will need. The details are in the appendices.

P1: The translation operator $\mathbb{T}_{\mathbf{a}_{2}}$ sends the Hamiltonian at $\mathbf{K}_{n}$ to the Hamiltonian at $\mathbf{K}_{n-1}$ (mod $\mathbf{G}_{1}$). Similarly for the points $\mathbf{K}^{\prime}_{n}$. Since the real-space Hamiltonian commutes with $\mathbb{T}_{\mathbf{a}_{2}}$, the spectrum must be identical at all the $\mathbf{K}_{n}$ points.

P2: The rotation $\mathbb{R}_{\pi}$ takes the set of $\mathbf{K}_{n}$ points to the set of $\mathbf{K}^{\prime}_{n}$ points. Since this is a symmetry of the Hamiltonian, the spectrum at the $\mathbf{K}^{\prime}_{n}$ points is identical to that at the $\mathbf{K}_{n}$ points. Together with P1, this means that it is sufficient to understand the spectrum at a single $\mathbf{K}_{n}$ point.

P3: A rotation by $\frac{2\pi}{3}$ of the destruction operator at an arbitrary point $\mathbf{k}$ in the BZ leads to a superposition of destruction operators at the points

[TABLE]

where $\mathbf{k}_{R}$ is simply $\mathbf{k}$ geometrically rotated by $\frac{2\pi}{3}$, and $\gamma=0\dots,q-1$. $\mathbf{G}_{1}$ and $\mathbf{G}_{2}$ are the reciprocal lattice vectors of the original lattice. The set of points $\mathbf{K}_{n}$ are taken into each other by this transformation, as are the points $\mathbf{K}^{\prime}_{n}$. The operator transformations for the fermion operators can be found in Appendix C. Note that $\mathbb{R}_{\frac{2\pi}{3}}$ preserves the sublattice index.

P4: The chiral symmetry $\mathbb{S}$ means that at any point $\mathbf{k}$ in the BZ the Hamiltonian can be written in block form

[TABLE]

where the q $A$-type sublattice sites have been listed first, and the q $B$-type sublattices have been listed next.

Now we are ready for the proof by contradication. Let us assume that there are no zero-energy states at a particular $\mathbf{K}_{n}$ point. Let us further assume that there are no degneracies in the spectrum, so there are $2q$ nondegenerate states.

Eq. (18) implies two facts. Firstly, any eigenstate of energy $E\neq 0$ is necessarily a superposition of $A$ and $B$ sublattices $[\psi_{A},\psi_{B}]^{T}$, with non-zero amplitudes on both. Secondly, for every eigenstate with energy $E\neq 0$, there is another eigenstate $[\psi_{A},-\psi_{B}]^{T}$ with energy $-E$. The orthogonality of these two eigenstates implies that each $E\neq 0$ eigenstate has equal probabilities on the $A$-type and $B$-type sublattices.

Let us consider an eigenstate of $H(\mathbf{K}_{n})$ at a particular $n$. For notational convenience we will drop the $\mathbf{K}$ in what follows, and refer to objects at $\mathbf{K}_{n}$ simply by the subscript $n$, e.g., $c_{A\alpha}(\mathbf{K}_{n})\equiv c_{A,\alpha,n}$. By assumption the eigenstate we consider has $E\neq 0$. We can write the destruction operator for this eigenstate as

[TABLE]

Now we apply the $\mathbb{R}_{\frac{2\pi}{3}}$ to this equation. Since $[\mathbb{R}_{\frac{2\pi}{3}},H]=0$, the result will be a superposition of operators corresponding to eigenstates at the same energy $E$ at all the $\mathbf{K}$-type points.

[TABLE]

Now focus on the $n=n^{\prime}$ term on the RHS. From P3 we know that $\mathbb{R}_{\frac{2\pi}{3}}$ does not mix the $A$ and $B$ sublattices. Thus, the restriction of $\mathbb{R}_{\frac{2\pi}{3}}$ to $n=n^{\prime}$ is a block diagonal $2q\times 2q$ matrix.

[TABLE]

where each of $R_{A}(n)$ and $R_{B}(n)$ are $q\times q$ matrices.

Applying this to Eq. (19) we see that $\psi_{A,\alpha,n}$ must be an eigenstate of $R_{A}(n)$, and $\psi_{B,\alpha,n}$ must be an eigenstate of $R_{B}(n)$, with the same eigenvalue. Note that if an eigenstate of $H(n)$ had zero energy, it need not have nonzero amplitudes in both $A$ and $B$ sublattices, and so could evade this conclusion.

By assumption, all the eigenstates have nonzero energy. Thus, all the eigenvalues of $R_{A}(n)$ and $R_{B}(n)$ must be identical. This leads to the conclusion that

[TABLE]

From the explicit forms of $R_{A}(n)$ and $R_{B}(n)$ obtained by restricting the expressions in Appendix C one easily obtains

[TABLE]

From this, one sees that

[TABLE]

This contradicts our conclusion of Eq. (22). Thus, at least some of the states at $\mathbf{K}_{n}$ must have zero energy. From the fact that the chiral symmetry implies that energies must occur in pairs of $\pm E$, an even number of states must have zero energy at any $\mathbf{K}_{n}$.

This shows that there must be band touching at the $\mathbf{K}_{n}$ and $\mathbf{K}^{\prime}_{n}$ points. Carrying out $\mathbf{k}\cdot\mathbf{p}$ perturbation theory around a $\mathbf{K}_{n}$ point, there is no symmetry reason for the first derivative to vanish, and thus the touchings will generically be linear. This completes our proof.

IV Stability to Small Perturbations

In the previous section we have shown that with the graphene lattice symmetry and $\mathbb{S}$, we have $2q$ Dirac nodes at the specific locations: $\{\mathbf{K}_{0},\cdots,\mathbf{K}_{q-1},\mathbf{K}^{\prime}_{0},\cdots,\mathbf{K}^{\prime}_{q-1}\}$ at zero energy. We now study the stability of these Dirac touchings to quadratic perturbations. Before turning to specific perturbations, we address this question in more general topological terms. Schnyder et al. (2008); Kitaev (2009) We know that time-reversal symmetry ($\mathbb{T}$) and particle-hole symmetry ($\mathbb{C}$) are both individually absent, but the composite of the two, the anti-unitary particle hole ($\mathbb{S}$), is present. A band insulator with the symmetry $\mathbb{S}$ would be in class AIII. Our model, with all the symmetries intact, has Dirac touchings and is thus not a band insulator. A band insulator in which $\mathbb{S}$ is broken (say by the introduction of same sub-lattice hopping or a sub-lattice energy difference) would be in class A. It is now understood that the stability of Dirac touchings can be explained by the classification of band insulators in one lower dimension. Ryu and Hatsugai (2002) The argument relies on considering the topological classification of the band insulator Bloch wavefunction on a ring surrounding the Dirac point. Turner and Vishwanath (2013) In one dimension band insulators in class AIII have a $\mathbb{Z}$ classification while in class A they have only trivial band structures. This integer winding number can be computed for each Dirac node by using a simple prescription. Referring to Eq. (18), one computes the winding number of the phase of the determinant of $M$ on a contour in the BZ around the band touching. The expression for the one dimensional winding number isSchnyder et al. (2008); Ryu et al. (2010); Fulga et al. (2012),

[TABLE]

Computing the winding number for the Dirac touchings we find that all $\mathbf{K}$ point have -1 and $\mathbf{K}^{\prime}$ have +1. That all of the $\mathbf{K}$ points have equal winding numbers follows from $\mathbb{T}_{\mathbf{a}_{2}}$ symmetry and that $\mathbf{K}^{\prime}$ have the opposite winding follows from the action of the $\mathbb{R}_{\pi}$ symmetry operation.

From these general topological considerations, we reach the following conclusions for the stability to perturbations:

Generically if $\mathbb{S}$ is broken the Dirac touchings get gapped (at least the argument above does not guarantee stability – below we study a few examples numerically to verify this). The resulting insulator will be in class A, with an integer Chern number.

What perturbation can open a gap if we preserve $\mathbb{S}$? If the perturbation fits in the MUC (so that the BZ is unchanged) the Dirac touchings are stable. We note that locations in the BZ may move if the perturbations reduce the symmetry operations from those present in an undistorted honeycomb lattice. Generally, if we preserve $\mathbb{S}$, small perturbations can open up the gap only if they are at a wavevector that connects Dirac points with opposite winding number. Then in the new smaller BZ (corresponding to the enlarged unit cell), opposite winding number Dirac points lie on top of each other. Encircling such double touchings will give no winding number, invalidating the topological argument for their protection.

We now consider specific lattice examples in which we can study how the Dirac equation gets gapped.

IV.1 ${\mathbb{S}}$-breaking

We first restrict our discussion to perturbations that preserve the MUC. To gap the Dirac touchings requires us to break the $\mathbb{S}$ symmetry. The simplest way to do this is to perturb with a staggered diagonal energy term in the Hamiltonian that has the same magnitude but opposite signs on the each of the sub-lattices. $H=V_{s}\sum_{n_{1},n_{2}}\left(d^{\dagger}_{A,n_{1},n_{2}}d_{A,n_{1},n_{2}}-d^{\dagger}_{B,n_{1},n_{2}}d_{B,n_{1},n_{2}}\right)$. Although this fits in the MUC it breaks some of the lattice symmetry, e.g. $\mathbb{R}_{\pi}$. A second way to break $\mathbb{S}$ symmetry is to include a same sublattice hopping with fixed range for all sites, e.g. a second neighbor hopping, $t_{2}$. We include it here with the correct phase from Eq. (10) corresponding to having a background uniform $B$-field. This perturbation has the feature of preserving every symmetry in the nearest neighbor hopping model except for $\mathbb{S}$. From the arguments made earlier, both perturbations are expected to open up a gap in the Dirac equation, leaving behind a 2-D insulating band structure in class A, characterized by an integer Chern number.

Since the two middle bands become the $n=0$ Landau level of graphene, we expect them to have a combined Chern number of -2. How this -2 is distributed between the two bands obtained after the gap opening perturbation is added depends on the details. The sign of the mass that gaps out a particular Dirac point also determines the transfer of Chern density between the two bands. Perturbations that preserve translations can only realize the total Chern numbers $C=0,-q,q$ because the $\mathbb{T}_{\mathbf{a}_{2}}$ symmetry forces the form of the Bloch Hamiltonian at all the $\mathbf{K}_{n}$ points to be the same, and also all the $\mathbf{K}^{\prime}_{n}$ points to be the same. The perturbation $V_{s}$ results in a trivial insulator, while the other two values of $C$ are realized by the $t_{2}$ perturbation. We have checked all the above assertions by computing the integer invariant numerically, i.e. by integrating the Berry curvature over the Brillouin zone. This is shown and discussed in Fig. 3.

Above we have studied two examples of perturbations that create different Chern numbers, 0 and $\pm q$. The Chern numbers of the bands produced can be any of the intermediate values $1,2\cdots q-1$ as well. This requires a perturbation that breaks $\mathbb{T}_{\mathbf{a}_{2}}$ though it may preserve the MUC.

IV.2 ${\mathbb{S}}$-preserving

To gap out the Dirac nodes with $\mathbb{S}$ preserved, the perturbation must break translational invariance with a momentum that connects Dirac touchings of opposite winding number. The simplest way to achieve this is to include as a perturbation a periodic modulation of the magnitude of the first neighbor hopping, with period corresponding to the Q-vector connecting the Dirac nodes. Since there are $q$ nodes with positive winding and $q$ with negative winding, there appear to be $q^{2}$ different possibilities. However, only $q$ different Q-vectors fit within the magnetic Brillouin zone leading to $q$ different reduced lattice periodicities. Once the nodes are gapped out and in the presence of ${\mathbb{S}}$ we end up with a band insulator in class AIII. Since these are all trivial insulators they are expected to be smoothly connected to each other without a gap closing.

V Conclusions

In conclusion we have studied the stability of the Dirac touchings in the $n=0$ Landau level in the Hofstadter limit when the external magnetic field is very strong, with $\frac{1}{q}$ quanta of flux going through each hexagon.

We started by deriving a formula in the optimal gauge for an arbitary range hopping, so that our problem always fits with $q$ unit cells of the honeycomb. Next, we have shown that the Dirac touchings require the sublattice symmetry ${\mathbb{S}}$ for their protection. Indeed we have proven that every tight binding model with ${\mathbb{S}}$, the correct flux and with the entire lattice symmetry of graphene intact will have $2q$ Dirac touchings at the same location as the nearest neighbor model.

Next, we considered perturbations to the $2q$ Dirac touchings. We showed from general topological stability arguments as well as specific hopping models that perturbatively breaking $\mathbb{S}$ or including a periodic potential that connects Dirac nodes with opposite winding number can gap the Dirac nodes out. All other perturbations preserve the Dirac nodes (though their location in the BZ may move). Of course if these perturbations are made large enough, some finite value of the perturbation may cause the gapping out of the linear touchings.

Introducing electron-electron interactions is known to produce a rich set of symmetry-broken phases in graphene for weak fields. Herbut (2006); Alicea and Fisher (2006); Kharitonov (2012); Feshami and Fertig (2016) There have been a few investigations into interaction effects in the Hofstadter regime, Jung and MacDonald (2009); Mishra et al. (2014); Buividovich et al. (2016); Mishra et al. (2016, 2017) but a full picture remains to be developed.

GM thanks David Mross for a vauable discussion. We acknowledge partial financial support from NSF DMR-1611161 (AD and RKK) and from NSF DMR-1306897 (AD and GM). We thank the hospitality of the Aspen Center for Physics (NSF grant no. 1607611) where this work was finalized for publication. GM also thanks the Gordon and Betty Moore Foundation for sabbatical support at MIT, and the Lady Davis Foundation for sabbatical support at the Technion.

Appendix A Graphene without $B$-field

In discussing the Dirac touchings in graphene with $\mathbb{T}$, Mañes et al. (2007) it is useful to think about the problem in two steps as we have done in our manuscript for the case when $\mathbb{T}$ is broken by an external magnetic field.

First, it is possible using symmetries to prove that for any hopping model that preserves the symmetry of the honeycomb lattice and time reversal, there are two independent Dirac nodes at K and K′. The symmetry argument does not rule out the existence of other additional Dirac nodes in the BZ that may coexist with the two Dirac nodes mandated by symmetry. Proof: Write the graphene Hamiltonian as $h=d_{x}({\bf k})\sigma_{x}+d_{y}({\bf k})\sigma_{y}+d_{z}({\bf k})\sigma_{z}$. A combination of ${\mathbb{R}_{\pi}}$ and $\mathbb{T}$ establish that $d_{z}=0$. Finally requiring $R_{2\pi/3}$ in addition forces $(d_{x}-id_{y})_{({\bf K+q)}}\approx v_{F}(q_{x}-iq_{y})$ and $(d_{x}-id_{y})_{({\bf K^{\prime}+q)}}\approx v_{F}(-q_{x}-iq_{y})$ at leading order. Using $\mathbf{\tau}$ as the valley Pauli matrix we obtain the low energy Hamiltonian as $h\approx\tau_{z}\sigma_{x}q_{x}+\sigma_{y}q_{y}$.

Now we turn to the perturbative stability of the Dirac touchings at K and K′, when the symmetry is lowered by either breaking ${\mathbb{T}}$ or some of the honeycomb lattice symmetries. To gap out the Dirac fermion we need perturbations that generate mass terms that anticommute with both $\tau_{z}\sigma_{x}$ and $\sigma_{y}$. If the translational invariance of the triangular Bravais lattice, ${\mathbb{R}}_{\pi}$ and $\mathbb{T}$ are present the touchings are stable to all perturbations, the only changes to the unperturbed $h$ is a movement of the Dirac touchings. Four mass terms can be added: $\tau_{z}\sigma_{z}$ and $\sigma_{z}$ break $\mathbb{T}$ and $\mathbb{R}_{\pi}$ leading to the Chern insulator and the trivial band insulator respectively. Haldane (1988) $\tau_{x}\sigma_{x}$ and $\tau_{y}\sigma_{x}$ are the last two, they break translational symmetry of the graphene lattice, which could arise e.g. from Kekule dimerization.

Our goal in this paper is to carry out this same two step program for the problem of the honeycomb lattice with $1/q$ flux. The significant increase in complication arises from the fact that the matrices that describe this problem are now $2q$ dimensional.

Appendix B Symmetry Operations in OG

In this section we study the various symmetries that are present in Eq. (11). We will study these operations by asking how the symmetry operations act on the lattice creation and destruction operators in the OG. While the explicit form of the transformations are gauge dependent, these explicit forms exist in every gauge.

Lattice symmetries include translations in the $\mathbf{a}_{1}$ and $\mathbf{a}_{2}$ directions, rotations by $2\pi/3$ about a site and rotation by $\pi$ about the center of a vertical bond. This set of four operations generates all the spatial symmetry operations present in the Hofstatder problem. We note here that mirror symmetries (and generally all improper rotations) which are present for the honeycomb lattice structure, are broken by the presence of the Peierls phases, since they reverse the direction of the magnetic fluxes. As noted above, because of the presence of the Peierls phases the lattice operations must be augmented by a gauge transformation from the naive operations one writes down in the absence of a magnetic field. In the OG they are,

[TABLE]

[TABLE]

where $m_{1}=-n_{1}-n_{2}$ and $m_{2}=n_{1}$, and

[TABLE]

It is easy to verify that each of the above four operations commutes with the Hamiltonian in the OG, Eq. (11).

Finally a very important symmetry for our purposes present in Eq. (11) is an anti-unitary version of the particle-hole symmetry,

[TABLE]

which is easily seen to commute with $H_{\rm nn}$. We note that the conventional time reversal operation $\mathbb{T}$ and the conventional unitary particle-hole $\mathbb{C}$ each reverse the direction of the magnetic field and are hence are absent as symmetries in the present problem. We can then understand that since $\mathbb{S}=\mathbb{TC}$ reverses the magnetic field direction twice it appears as a symmetry of our problem.

Appendix C Symmetries in $\mathbf{k}$ Space

We choose the following periodicity conditions on our fermion operators in the Brillouin zone.

[TABLE]

Real Space translations,

[TABLE]

where, $[\alpha+1]=(\alpha+1)mod~{}q$ and $\alpha\in[0,q-1]$.

Rotation by $\pi$ about a bond center,

[TABLE]

$\frac{2\pi}{3}$ rotation about an $A$ lattice point mixes multiple $\mathbf{k}$ points,

[TABLE]

where $\mathbf{k}_{R}$ is $\mathbf{k}$ rotated by $\frac{2\pi}{3}$.

For $q$ an odd number, we have:

[TABLE]

For $q$ an even number, we have:

[TABLE]

The expressions for the Dirac points, and the restriction of the rotation operators particular Dirac points, naturally fall into two classes, those for $q$ odd, and those for $q$ even.

Appendix D Dirac points and $\mathbb{R}_{\frac{2\pi}{3}}$ for q odd

The Dirac points are,

[TABLE]

where, $n\in[0,q-1]$ and $\mathbf{K}^{\prime}_{n}=-\mathbf{K}_{n}$.

[TABLE]

[TABLE]

where,

[TABLE]

Appendix E Dirac points and $\mathbb{R}_{\frac{2\pi}{3}}$ for q even

The location of the Dirac points are

[TABLE]

where, $n\in[0,q-1]$ and $\mathbf{K}^{\prime}_{n}=-\mathbf{K}_{n}$. Effect of all other operation remain the same as in the odd $q$ case except for the $2\pi/3$ rotations.

[TABLE]

with,

[TABLE]

Appendix F $R_{A}(n)$ and $R_{B}(n)$

The matrix that rotates the wavefunction into itself (multiplying with $\sqrt{q}$ makes it unitary), for $q$ odd and $n=\frac{q-1}{2}$:

[TABLE]

Similarly, for $q$ even and $n=\frac{q}{2}$:

[TABLE]

Appendix G Rhim-Park Wavefunction

We can solve the wave function for the nearest-neighbor Hamiltonian using the methods used by Rhim and Park Rhim and Park (2012). Using Bloch’s theorem we can say that the wave function must have the form

[TABLE]

For the zero energy eigenstate from the Hamiltonian we can write,

[TABLE]

Thus using recursion relation we can write,

[TABLE]

where $k_{1}=\mathbf{k}\cdot\mathbf{a}_{1}$ and $k_{2}=\mathbf{k}\cdot\mathbf{a}_{2}$. Now from the self-consistency for Bloch functions gives the condition on the $\mathbf{k}$ momentum values at which the zero-energy states exist. The conditions are,

[TABLE]

The solutions forms a honeycomb lattice in the momentum space. They consist of two sets,

[TABLE]

and

[TABLE]

Where $j_{1},j_{2}$ can be any integer.

Appendix H Symmetry action in the low energy space

Using the Rhim-Park wave function we can derive the action of the symmetry operations in the low energy space.

Let us start with translations.

[TABLE]

[TABLE]

Where,

[TABLE]

Now, consider the $\pi$ rotation about center of a vertical bond $\mathbb{R}_{\pi}$.

[TABLE]

Now for the $\frac{2\pi}{3}$ rotation $\mathbb{R}_{\frac{2\pi}{3}}$. This is the most complicated of all because it maps a particular Dirac point to a linear combination of all Dirac points with the same winding number. The results we quote below are empirical in the sense that we have not been able to prove them: rather, we fitted the action of $\mathbb{R}_{\frac{2\pi}{3}}$ on the Rhim-Park wavefunctions to an analytic form, and checked them for many values of $q$.

[TABLE]

where,

[TABLE]

Finally the action of the chiral symmetry on the low energy subspace.

[TABLE]