Generalized eigenproblem without fermion doubling for Dirac fermions on a lattice

The spatial discretization of the single-cone Dirac Hamiltonian on the surface of a topological insulator or superconductor needs a special"staggered"grid, to avoid the appearance of a spurious second cone in the Brillouin zone. We adapt the Stacey discretization from lattice gauge theory to produce a generalized eigenvalue problem, of the form ${\mathcal H}\psi=E {\mathcal P}\psi$, with Hermitian tight-binding operators ${\mathcal H}$, ${\mathcal P}$, a locally conserved particle current, and preserved chiral and symplectic symmetries. This permits the study of the spectral statistics of Dirac fermions in each of the four symmetry classes A, AII, AIII, and D.


Introduction
Three-dimensional topological insulators are Nature's way of working around the Nielsen-Ninomiya no-go theorem [1], which forbids the existence of a single species of massless Dirac fermions on a lattice. The fermion doubling required by the theorem is present in a topological insulator slab, but the two species of Dirac fermions are spatially separated on opposite surfaces [2,3]. On each surface the two-dimensional (2D) Dirac Hamiltonian emerges as the effective low-energy Hamiltonian, with a single Dirac cone at k = (k x , k y ) = 0.
Since it is computationally expensive to work with a three-dimensional (3D) lattice, one would like to be able to discretize the 2D Dirac Hamiltonian, without introducing a second Dirac cone. We can draw inspiration from lattice gauge theory, where a variety of strategies have been developed to avoid fermion doubling [4,5]. The condensed matter context introduces its own complications, notably the lack of translational invariance and breaking of chiral symmetry by disorder and boundaries.
In Ref. 6 it was shown how the transfer matrix of the Dirac equation in a disorder potential can be discretized without fermion doubling. This allows for efficient calculation of the conductance and other transport properties in an open system [7][8][9]. Here we apply the same approach to the Hamiltonian of a closed system, in order to study the spectral statistics.
The Nielsen-Ninomiya theorem forbids a local discretization of the eigenvalue problem H D ψ = Eψ without fermion doubling and without breaking the chiral symmetry relation (1.2) One way to circumvent the no-go theorem, is to abandon the locality by introducing longrange hoppings in the discretized Dirac Hamiltonian [10]. Here we follow an alternative route, following Stacey [11], which is to work with a generalized eigenvalue problem with local tight-binding operators H and P on both sides of the equation. Going beyond Ref. 11, we transform the operators H and P such that they remain, respectively, Hermitian and positive definite in the absence of translational invariance. This favors a stable and efficient numerical solution, and moreover guarantees that the resulting spectrum is real, not only in the continuum limit but at any grid size.
A key feature of our approach, compared with the more familiar approaches of Wilson fermions [12] and Susskind fermions [13], is that both the chiral symmetry (1.2) is preserved and the symplectic time-reversal symmetry [14] σ y H * D σ y = H D . (1.4) This also implies the conservation of the product of the chiral and symplectic symmetries, which is a particle-hole symmetry, To demonstrate the capabilities of our approach we calculate the spectral statistics of a disordered system and show how the numerics distinguishes broken versus preserved chiral or symplectic symmetry in each of the four symmetry classes of random-matrix theory [15].
The outline of the paper is as follows: In the next section we formulate the generalized eigenproblem, first following Stacey [11] for a translationally invariant system, and then including disorder. The symmetrization that produces a Hermitian H and positive definite P is introduced in Sec. 3. The locality of the discretization scheme is demonstrated by the construction of a locally conserved current in Sec. 4. By applying different types of disorder, in scalar potential, vector potential, or mass, we can access the different symmetry classes and obtain the characteristic spectral statistics for each, as we show in Sec. 5. We conclude in Sec. 6.
2 Construction of the generalized eigenproblem

Staggered discretization
If we discretize the Dirac Hamiltonian (1.1) on a lattice (lattice constant a), the replacement of the momentum k by a −1 sin ka produces a second Dirac cone at the edge of the Brillouin zone (k = π/a). To place our work into context, we summarize methods to remove this spurious low-energy excitation.
If one is willing to abandon the locality of the Hamiltonian, one can eliminate the fermion doubling by a discretization of the spatial derivative that involves all lattice points, df /dx → n (−1) n n −1 f (x − na). The resulting dispersion remains strictly linear in the first Brillouin zone. This discretization scheme goes by the name of slac fermions [10] in the high-energy physics literature. It has recently been implemented in a condensed matter context [16].
An alternative line of approach preserves the locality at the expense of a symmetry breaking. The simplest way is to couple the top and bottom surfaces of the topological insulator slab [17,18]. The coupling adds a momentum dependent mass term µσ z (1−cos ka) which gaps out the second cone, while breaking both chiral symmetry and symplectic symmetry. This is the Wilson fermion regularization of lattice gauge theory [12,19]. The product of chiral and symplectic symmetry is preserved by Wilson fermions, which may be sufficient for some applications [20,21].
It is possible to maintain the chiral symmetry by discretizing the Dirac Hamiltonian on a pair of staggered grids. Much of the lattice gauge theory literature is based on the Susskind discretization [13], which applies a different grid to each of the two components of the spinor wave function ψ. On a 2D lattice it reduces the number of Dirac cones in the Brillouin zone from 4 to 2. Chiral symmetry is preserved, but symplectic symmetry is broken by the Susskind discretization (see App. A).
Hammer, Pötz, and Arnold [22,23] have developed an ingenious single-cone discretization method for the time-dependent Dirac equation. As in the Susskind discretization, different grids are used for each of the spinor components, but these are staggered not only in space but also in time. While this method is well suited for dynamical simulations [24,25], it is not easily adapted to energy-resolved spectral studies.
An altogether different approach, introduced by Stacey [11,26], is to evade the fermiondoubling no-go theorem by the replacement of the conventional eigenvalue problem H D ψ = Eψ by a generalized eigenproblem U ψ = EΦψ. There is now no obstruction to having a local U and Φ and also preserving chiral and symplectic symmetry.
The Stacey discretization of the transfer matrix was implemented in Ref. 6. In what follows we show how to apply it to the Hamiltonian, to solve the time-independent Dirac equation on a 2D lattice. In the next subsection we first summarize the results of Ref. 11 for a translationally invariant system, and then will present the modifications needed to apply the method in the presence of a disorder potential.

Translationally invariant system
We seek to discretize the Dirac equation H D ψ = Eψ on a 2D square lattice (lattice constant a). We denote the discretized wave function by ψ n , with n = (n x , n y ) ∈ Z 2 labeling the lattice points at n x e x + n y e y . For ease of notation we will henceforth set v F , , and a to unity. Staggered discretization a la Stacey means that the wave function and its spatial derivatives are evaluated on a displaced lattice with sites at the center of the unit cells of the original lattice (see Fig. 1). The discretization rules are: In distinction to Susskind staggering, the same discretization applies to each spinor component.
In momentum representation, ψ(k) = n ψ n e −ik·n , the discretized Dirac equation with the k-dependent operators The dispersion relation has a single Dirac point at k = 0. The Dirac point at the edge of the Brillouin zone has been converted into a pole by the Stacey discretization. symmetry symplectic chiral particle-hole class

Including a disorder potential
We break translational invariance by including in the Dirac equation a spatially dependent scalar potential V σ 0 , vector potential A x σ x + A y σ z , and mass M σ z , The electron charge e is set to unity in what follows. The Pauli matrices σ = (σ x , σ y ) and σ z act on the spin degree of freedom, with σ 0 the 2 × 2 unit matrix.
On the surface of a topological insulator the mass term represents a perpendicular magnetization. Alternatively, we can consider a 2D topological superconductor with chiral p-wave pair potential, described by the Bogoliubov-de Gennes (BdG) Hamiltonian The Pauli matrices now act on the electron-hole degree of freedom, electrons and holes are coupled by the pair potential ∝ v ∆ . Since this coupling is linear in momentum k, the quadratic kinetic energy k 2 /2m can be neglected near k = 0. The difference V − E F of electrostatic potential V and Fermi energy E F then plays the role of the mass term M in Eq. (2.5). The low-energy physics of the problem is governed by three symmetry relations, the chiral symmetry (1.2), the symplectic symmetry (1.4), and the particle-hole symmetry (1.5). Chiral symmetry is preserved by A and broken by V or M . Symplectic symmetry is preserved by V and broken by M or A. If at least two of the three potentials V, M, A are nonzero all symmetries of the Dirac Hamiltonian are broken. Finally, if V = 0, A = 0 while M = 0 the particle-hole symmetry (1.5) remains. Table 1 summarizes the symmetry classification [15].
The inclusion of the vector potential requires a separate consideration, in order to preserve gauge invariance. We delay that to Sec. 4, at first we only include V and M .
To incorporate the spatially dependent terms in the discretization scheme we write the operators U and Φ in the position basis. In view of the identity we have Ω ss = n ss |n n| + s|n n + e x | + s |n n + e y | + |n n + e x + e y | .
For later use we also define the factorization Φ = Φ x Φ y , with commuting operators Φ x , Φ y given by (2.10) In these equations the ket states |n refer to sites on the displaced lattice (open lattice points in Fig. 1), while the bra states n| refer to sites on the original lattice (closed lattice points). The inner product is defined such that the two sets of eigenstates of position are orthonormal, n |n = δ n,n .
We define the potential and mass operators, Eq. (2.12) is a generalized eigenvalue problem, with operators on both sides of the equation. Neither operator is Hermitian. This is problematic in a numerical implementation, and we will show in the next section how to resolve that difficulty.

Symmetrization of the generalized eigenproblem
We wish to rewrite Eq. (2.12) in the form Hψ = EPψ, with Hermitian H and Hermitian positive definite P. Such a symmetrization of the generalized eigenvalue problem allows for a stable and efficient numerical solution [27][28][29]. Moreover, it guarantees real eigenvalues E and eigenvectors ψ E that satisfy the orthogonality relation We multiply both sides of Eq. (2.12) by Φ † and note that Φ † U is a Hermitian operator. In position basis it reads We thus arrive at the generalized eigenproblem In the translationally invariant case the operators H and P are given by H = 1 2 σ x (1 + cos k y ) sin k x + 1 2 σ y (1 + cos k x ) sin k y , P = 1 4 (1 + cos k x )(1 + cos k y ). Both operators are Hermitian and P is also positive semi-definite. Moreover, P is positive definite if the edges of the Brillouin zone (k x or k y equal to ±π) are excluded from the spectrum. To ensure that, we can choose an odd number N x , N y of lattice points with periodic boundary conditions in the x-and y-directions (or alternatively, even N x , N y with antiperiodicity). By way of illustration, we work out the expectation value (3.5) In the translationally invariant case, the effective Hamiltonian reduces simply tõ H = 2σ x tan(k x /2) + 2σ y tan(k y /2).

Locally conserved particle current
In real space the effective Hamiltonian (3.5) produces infinitely long-range hoppings, as in the slac fermion discretization [10,16]. The transformation to the generalized eigenproblem (3.2) restores the locality of the hoppings. One might wonder whether there is a physical content to this mathematical statement. Yes there is, as we show in this section the Stacey discretization allows for the construction of a locally conserved particle current.
We define the particle number corresponding to the density operator With reference to the two staggered grids in Fig. 1, the particle density on an open lattice point n is given by the norm squared of the average of the wave function on the four adjacent closed lattice points, ψ|ρ(n)|ψ = | 1 4 (ψ n + ψ n+ex + ψ n+ey + ψ n+ex+ey )| 2 . The current density operator is given by or equivalently, in terms of the operators Φ x , Φ y defined in Eq. (2.10). The current density in the state ψ then takes the form ψ|j x (n)|ψ = 1 4 (ψ n + ψ n+ey ) † σ x (ψ n + ψ n+ey ), ψ|j y (n)|ψ = 1 4 (ψ n + ψ n+ex ) † σ y (ψ n + ψ n+ex ). (4.6) The current density at an open lattice point is evaluated by averaging the wave function at the two nearby closed lattice points connected by an edge perpendicular to the current flow.
The local conservation law is derived in App. B. Knowledge of the current operator allows us to introduce the vector potential operator A = n A n |n n| such that lim A→0 ∂H ∂A n = j(n).
In App. C we check that the Hamiltonian (4.9) is gauge invariant to first order in A.
Higher order terms are nonlocal and we will not include them.

Spectral statistics
We have tested the validity and capability of the generalized eigenvalue problem by comparing the spectral statistics with predictions from random-matrix theory (RMT). Similar tests for different methods to place Dirac fermions on a lattice have been reported in the particle physics literature [30][31][32].
We have solved the generalized eigenproblem on a square lattice of size N x × N y . Antiperiodic boundary conditions in the x-and ydirection account for the π Berry phase accumulated by the spin when it makes one full rotation. The dimensions N x , N y are even to ensure a positive definite Φ (no zero-mode in the spectrum). The spectrum was calculated for 5 · 10 4 realizations of a random disorder, chosen independently on each site from a uniform distribution in the interval (−δ, δ).
To access the four symmetry classes from Table 1 we took • A x , A y ≡ 0 and random V, M with δ = 15/ √ 2 for class A; • M, A x , A y ≡ 0 and random V with δ = 15 for class AII; • V, M ≡ 0 and random A x , A y with δ = 1 4 √ 2 for class AIII; • V, A x , A y ≡ 0 and random M with δ = 15 for class D.  Table 1. The red dashed line is the prediction (5.2) from random-matrix theory in the presence of symplectic symmetry (β = 4) and in its absence (β = 2).
The relatively weak disorder in class AIII was chosen in view of the linearization in the vector potential. For that case we took N x = N y = 150, in the other symmetry classes with stronger disorder we took N x = N y = 100. Symmetry class D is insulating for weak disorder in the mass M ∈ (−δ, δ), it undergoes a metal-insulator transition at δ c = 3.44 [7]. This is the thermal metal phase of a topological superconductor [33]. The thermal metal can be reached by vortex disorder, as in the network model studied in Ref. 34, or it can be reached by electrostatic disorder in the BdG Hamiltonian (2.6), as in the tight-binding models studied in Refs. 7, 35. Here we follow the latter approach, taking δ = 15 much larger than δ c , so that we are deep in the metallic regime.
In Fig. 2 we show the probability distribution of the level spacing δE in the bulk of the spectrum, far from E = 0, where the average spacing E is energy independent. We compare with the Wigner surmise from RMT [36], P (s) = 32 π 2 s 2 e −4s 2 /π in class A, AIII, D, with s = δE/ δE . The characteristic difference between the two distributions is the decay ∝ s β for small spacings, with β = 4 in the presence of symplectic symmetry, while β = 2 in its absence. (The case β = 1 of RMT is not realized in a spin-full system.) In Fig. 3 we make a similar comparison for the density of states near E = 0. In class A and AII the ensemble averaged density of states ρ(E) is flat in a broad energy range around E = 0. Chiral symmetry in class AIII introduces a linear dip in the density of states, while particle-hole symmetry in class D introduces a quadratic peak. The RMT Figure 3: Density of states in the four symmetry classes, calculated numerically from the discretized Dirac Hamiltonian (blue solid lines) and compared with the RMT prediction (5.3) (red dashed lines). Chiral symmetry introduces a linear dip (class AIII), while particle-hole symmetry introduces a quadratic peak (class D).
predictions are [37,38] with ε = E/ δE . The mean level spacing δE is computed away from E = 0. The good agreement between the numerical results from the disordered Dirac equation and the RMT predictions, evident in Figs. 2 and 3, is reached without any adjustable parameter. Remaining discrepancies are likely due to a dynamics that is not fully chaotic. (In particular, incipient localization can explain the shift to smaller spacings noticeable in Fig. 2.) The computer code to reproduce this data is provided [39].

Conclusion
In conclusion, we have developed and implemented a lattice fermion Hamiltonian that, unlike the familiar Wilson fermion and Susskind fermion Hamiltonians [12,13], preserves both chiral symmetry and symplectic symmetry while avoiding fermion doubling. Our approach is a symmetrized version of Stacey's generalized eigenvalue problem [11], which allows for the construction of a locally conserved particle current. To demonstrate the universal applicability of the lattice fermion Hamiltonian we have shown how it can reproduce the characteristic spectral statistics for each of the four symmetry classes of Dirac fermions.
We mention three topics for further research. Firstly, we have only succeeded in including the vector potential in a gauge invariant way to first order, so for a flux through a unit cell that is small compared to the flux quantum. Is it possible to remove this limitation? Secondly, can we extend the approach to discretize time as well as space? And thirdly, can we incorporate boundary conditions without breaking the local current conservation?

A Susskind discretization breaks symplectic symmetry
The staggered discretization of the 2D Dirac equation a la Susskind [13] produces a conventional eigenvalue problem, with a local Hamiltonian. There is a single Dirac cone in 1D but there are 2 Dirac cones in 2D. Chiral symmetry is preserved, but symplectic symmetry is broken. To contrast this with the symplectic-symmetry-preserving single-cone Stacey discretization used in the main text, we give a brief description of the Susskind discretization, first in 1D and then in 2D.
In 1D the staggering refers to the prescription that the derivative of the A component of the spinor ψ = (ψ A , ψ B ) is calculated at x = n + 1/2, while the derivative of the B component is calculated at x = n − 1/2. Hence the term k x σ x in the Dirac Hamiltonian is substituted by The exponential e ∂x , with ∂ x = ∂/∂x, is the translation operator: e ∂x ψ(x) = ψ(x + 1).
In momentum representation, ∂ x → ik x , the discretized Hamiltonian reads The corresponding dispersion relation has a single Dirac cone at k x = 0 in the Brillouin zone −π < k x ≤ π. The 2D generalization is The resulting dispersion relation, We take for O the density operator (4.2), soρ(n) = |n n|. This projector commutes with the operators V and M inH, what remains is the commutator with U Φ −1 : In the last equality we used that Φ † U = U † Φ.
In terms of the current operator

C Gauge invariant vector potential
To include the vector potential A(r) in a gauge invariant way in the discretized Dirac equation, we follow the procedure of minimal coupling: We first discretize without a vector potential, then perform a U(1) gauge transformation on the lattice, and finally replace the gradient of the phase field by the vector potential. We define the gauge field operator where the line integral of the vector potential is taken along a lattice bond. With this prescription the substitution can also be applied to vector potentials that do not derive from a gauge field.