Dispersion relation and spectral range of Karsten-Wilczek and Borici-Creutz fermions

We investigate some properties of Karsten-Wilczek and Borici-Creutz fermions, which are the best known varieties in the class of minimally doubled lattice fermion actions. Our focus is on the dispersion relation and the distribution of eigenvalues in the free-field theory. We consider the situation in two and four space-time dimensions, and we discuss how properties vary as a function of the Wilson-like lifting parameter $r$.


Introduction
The choice of any lattice fermion action is a bit of a compromise. Ideally, one would want to realize ultra-locality, chiral symmetry, and absence of doublers (in addition to a correct continuum limit, of course). But these are precisely the ingredients which, according to the Nielsen-Ninomiya theorem, cannot possibly coexist [1,2]. Furthermore, in real life computational expedience is an important criterion. It follows that the choice of a lattice action which is well-suited to the specific needs of a planned numerical investigation is an important decision which impacts the subsequent analysis of the lattice data in a profound manner.
Staggered fermions and Wilson fermions represent two popular choices in this context. Staggered fermions put a focus on ultra-locality and chiral symmetry, at the expense of having 4 species in the continuum [3]. Wilson fermions, on the other hand, prioritize ultra-locality and absence of doublers, at the expense of breaking chiral symmetry [4].
Staggered fermions seem ideally suited to study theories with four degenerate fermions (or a multiple thereof) in the continuum. The details of taste-breaking [5][6][7][8][9] and the noncommutativity of the continuum limit (a → 0) and the chiral limit (m → 0) [10] impose practical difficulties, but to the best of our knowledge there is no concern about the theoretical soundness of this formulation of QCD with N f ∈ 4 N dynamical flavors.
Things are different, if one desires to study QCD with fewer than 4 degenerate fermions, such as real-life QCD where m d , m u , m s , m c are pairwise different. Given the continuum argument that N f degenerate fermions would raise the functional determinant of a single species to the N f -th power, Marinari, Parisi and Rebbi suggested by means of "reverse engineering" that one would take the square-root of the staggered determinant to simulate 2 degenerate fermions or the quarter-root for a non-degenerate fermion species [11].
In practice, it seems the approach of "rooting" the staggered determinant yields convincing numerical results for real-life QCD, with small statistical error bars even at physical values of the quark masses m d , m u , m s , m c [12][13][14], and with some field-theoretic underpinning through rooted staggered chiral perturbation theory [15][16][17]. Nevertheless, this approach has been criticized [18,19] on the grounds of the argument that -in the presence of the cut-off -no local theory can be constructed (or has been constructed) that would implement exactly (i.e. down to machine precision) the square-root of the staggered determinant as the determinant of a valid 2 species formulation. There have been several reviews of the issue at lattice conferences [20][21][22][23][24] which essentially collected pieces of evidence in favor of the approach. But it holds true that no strict mathematical proof has been found, and no side has been able to convince the opponent.
Given this situation it is natural to ask whether a local formulation with chiral symmetry and just 2 species (the minimum required by the Nielsen-Ninomiya theorem) would shed some light on the issue. Since staggered fermions emerge from the naive formulation by an ingenious procedure of "thinning out" the degrees of freedom by a factor 2 d/2 in d space-time dimensions, one might dream of a similarly ingenious second step that would reduce the degeneracy from 4 to 2 in d = 4 dimensions. However, the eigenvalue spectrum of staggered fermions on interacting backgrounds shows a 4-fold near-degeneracy (i.e. no exact degeneracy) [25,26], and this means that such a second reduction step cannot possibly take place.
However, there is a better approach. It is based on adding an extra term (of mass-dimension five) to the naive action which lifts 14 of the 16 species in d = 4 dimensions, albeit with the important difference to the Wilson term that it does not break chiral symmetry. Such "minimally doubled fermions" have been proposed by Karsten [27] and Wilczek [28], and later by Creutz [29] and Borici [30]. More recently, yet another variety with "twisted ordering" has been proposed by Creutz and Misumi [31]. From the viewpoint of computational efficiency, all these formulations augur well, since they are ultra-local with on-axis links only. In the literature issues of mixing with lower-dimensional operators have been addressed, and how one can defeat them with appropriate tuning strategies [32][33][34][35][36][37]. Also the consistency of these formulations with the index theorem has been verified [38]. Furthermore, some minimally doubled actions have been shown to possess an extra "mirror fermion" symmetry [39], and it has been demonstrated that the Karsten-Wilczek determinant is invariant under all of the discrete symmetries [40]. Still, some basic features of minimally doubled fermion actions have hardly been explored, for instance the respective quark-level free-field dispersion relations (including the cut-off effects on a heavy quark mass) and spectral bounds.
In this article we try to fill some of these gaps for Karsten-Wilczek (KW) and Borici-Creutz (BC) fermions. To understand how these formulations differ from the naive one we think it is useful to introduce a lifting parameter r, similar to what is done for Wilson fermions. Hence for r = 0 we start with the naive action, and we expect to see a cascade of species reductions as r increases, eventually realizing 2 species at r = 1. Since chiral symmetry holds throughout, the Nielsen-Ninomiya theorem demands that the reduction proceeds in steps of (integer multiples of) 2. The remainder of this article is organized as follows. In Sec. 2 we give a brief review of the situation with naive and Wilson fermions, mostly to specify our notation. The results for KW fermions are presented in Sec. 3, and those for BC fermions in Sec. 4. We summarize our findings in Sec. 5, and more lengthy calculations are arranged in appendices A-E.

Naive and Wilson fermions
Throughout this article ∂ µ and ∂ * µ denote the discrete forward an backward derivative, respectively, and ∇ µ = (∂ µ + ∂ * µ )/2 is the symmetric derivative. These operators are gauged in the obvious manner; for instance the symmetric covariant derivative is where U µ (x) is the parallel transporter from x +μ to x, andμ denotes a times the unit-vector in direction µ. Similarly, µ = ∂ * µ ∂ µ = ∂ µ ∂ * µ denotes the second discrete derivative, that is in the presence of a gauge field U µ (x). With this notation the naive Dirac operator is defined as where the anti-hermitean behavior ∇ † µ = −∇ µ and the anti-commutation property {γ µ , γ 5 } = 0 of the hermitean γ-matrices imply that the naive Dirac operator is γ 5 -hermitean, i.e. γ 5 D nai γ 5 = D † nai . In the free-field limit this operator assumes a diagonal form in momentum space, which again highlights the anti-hermitean nature of the derivative (momentum) term. The Wilson Dirac operator follows by adding a hermitean and positive semi-definite term of dimension 5 to the naive Dirac operator The hermitean behavior † µ = µ together with the properties used in the naive case imply that the Wilson (W) operator is γ 5 -hermitean, i.e. γ 5 D W γ 5 = D † W . An unpleasant feature is that the term µ µ (x, y) mixes (on interacting gauge backgrounds) with the identity. As a result, the bare mass m in (5) gets renormalized, and chiral symmetry is broken [41]. In the free-field limit the Wilson operator assumes a diagonal form in momentum space, which again highlights the anti-hermitean and hermitean positive semi-definite nature of the two terms, respectively. Specifically for r = 1 the 15 unwanted species do not propagate in any of the on-axis directions [41].   Wilson, L/a=32, r=1 With these expressions in hand, we are in a position to study the quark-level free-field dispersion relations. Both for naive and Wilson fermions, the energy aE can be given as an analytic function of the spatial momentum a p, see App. A for a brief account of this standard calculation. The results are shown in Fig. 1 for the naive formulation and in Fig. 2 for the Wilson action at r = 1. In either figure the situation at am = 0 and at am = 0.5 is compared to the respective continuum dispersion relation. Apart from the unwanted zeros (or minima) at a| p| = π, √ 2π, √ 3π the naive action features well for small enough momenta. In particular at a| p| = 0 it features better than the Wilson action, since the gap to the continuum curve is smaller. This can be understood on analytical grounds, too. The energy at zero momentum is nothing but the heavy quark mass. As detailed in App. B, it is known to be afflicted with cut-off effects O((am) 2 ) in the naive case, but O(am) in the Wilson case.
From expression (4) or (6) one finds the eigenvalues of the free-field operators. Since each γ µ (µ = 1, .., 4) has eigenvalues ±1, it follows that the upper end of the massless naive eigenvalue spectrum is realized for ap µ along the hyperdiagonal (1, 1, 1, 1), or flipped versions thereof, so |Im(λ nai )| ≤ 2 (7) and similarly it follows that the Wilson eigenvalue spectrum is contained in the rectangle The complex eigenvalue spectrum of the Wilson operator at am = 0 is shown in Fig. 3. The symmetry about the real axis reflects the pairing property imposed by the γ 5 -hermiticity. As is well known, the five branches in the Wilson eigenvalue spectrum correspond to species with multiplicities 1, 4, 6, 4, 1, and chiralities +, −, +, −, +, respectively [41]. In total 8 species thus have correct chirality, and 8 have opposing chirality. The free-field eigenvalue spectrum of the naive (or staggered) operator follows by horizontally projecting the λ W onto the imaginary axis (and reducing degeneracies by a factor 4 in the staggered case). Under this operation the separation between right-chirality and opposite-chirality species gets lost, and this feature will carry on to minimally doubled actions.

Karsten-Wilczek fermions
The Karsten-Wilczek proposal is to restrict the Wilson term in (5) to the spatial components with an extra factor iγ 4 to make it anti-hermitean and anti-commuting with γ 5 [27,28]. As a result the Karsten-Wilczek (KW) operator is γ 5 -hermitean, i.e. γ 5 D KW γ 5 = D † KW . An issue discussed in the literature is that γ 4 i i (x, y) mixes (on interacting backgrounds) with γ 4 [32][33][34][35][36][37]. In the free-field limit the KW operator assumes a diagonal form in momentum space, which again highlights the anti-hermitean nature of either term. This formulation was shown to have 2 species for r = 1 in the original works [27,28], but how this number decreases from 16, at r = 0, to the minimally doubled value has, to the best of our knowledge, not been investigated. We find that the number of species is reduced in three steps (in d = 4 dimensions). At r = 1/6, 1/4, 1/2 the number of species is reduced by 2, 6, 6, respectively, so the species chain is 16 → 14 → 8 → 2. See App. C for details, e.g. the situation with d = 4. Starting from eqn. (10) one can work out the free-field dispersion relation of KW fermions, see App. A for details. For a given momentum configuration a p the Euclidean energy aE is, in general, complex valued, and its real part is plotted in Fig. 4 for r = 1. Again, two values of the quark mass are used, am = 0 and am = 0.5. In either case the KW dispersion relation follows the continuum curve faithfully, out to momentum values a| p| 1. In particular at a| p| = 0 it features much better than the Wilson action, reminiscent of the naive action. This is not a coincidence, since the rest energy has the same functional dependence on am as the naive action, see App. B for details. In other words, cut-off effects on this quantity start at O((am) 2 ) only, unlike the O(am) signature of Wilson fermions.
From eqn. (10) one finds the eigenvalues of the free-field KW operator. The result for r = 1 and am = 0 is shown in Fig. 5. On a 32 4 lattice one finds 4 · 32 4 purely imaginary and γ 5 -paired eigenvalues, as expected. Plotting the eigenvalues against the index means that the inverse slope encodes for the density of the λ KW /i on the imaginary axis.
It is instructive to repeat this for a series of r values; the result is shown in Fig. 6, with vertical lines marking the abscissa values r = 1/6, 1/4, 1/2 where the number of species changes.   The spectral range is seen to increase with growing r. In addition, the low-energy end of the eigenvalue spectrum seems unstable for small r, but stable in a broad range around r = 1. Given Fig. 6, one may wonder about the existence of an analytic function which describes the upper end as a function of r. In App. E we derive the free-field spectral bound in d = 4 dimensions. Hence at r = 1 the imaginary parts λ KW /i cover the range [−7, 7], to be compared to [−2, 2] for naive and staggered fermions. On the other hand, the smallest non-zero KW eigenvalue is found in essentially the same place 1 as the smallest staggered eigenvalue, see Fig. 7. This amounts to an enhancement of the condition number of D † D, compared to the staggered formulation at the same am, by a factor up to 3.5 2 = 12.25 (in the chiral limit).
In App. C we discuss, for d = 4, how the number of species is reduced by 6 at r = 1/ √ 3, and by 8 at r = 1/ √ 2; so the species chain is 16 → 10 → 2. Of course, the number of species is unchanged by a sign flip of r.
In two (Euclidean) space-time dimensions it is customary to use γ 1 = σ 1 , γ 2 = σ 2 . Similarly, the chirality operator is defined as the operators (9) and (16) are seen to take the simple form which shows that the BC operator is not a symmetrized form of the KW operator; it has an extra term. This explains why we deviate, in the sign of the mass-dimension 5 term in eqns. (16,18), from the literature. With our convention the joint terms in eqns. (23,24) have like sign. Starting from eqn. (18) one can work out the free-field dispersion relation of BC fermions, see App. A for details. For a given momentum configuration a p the Euclidean energy aE is, in general, complex valued, and its real part is plotted in Fig. 8 for r = 1. Again, two values of the quark mass are used, am = 0 and am = 0.5. In either case the BC dispersion relation follows the continuum curve reasonably well, out to momentum values a| p| 1 2 , a range slightly narrower than what was found in the KW case. Specifically at a| p| = 0 it features much better than the Wilson action, though a little worse than the KW action. In App. B the rest energy  Borici-Creutz, L/a=32, r=1 Figure 9: Free-field eigenvalue spectrum of the BC operator at r = 1. The imaginary part is plotted against the index.
of a heavy BC fermion is found to have a contribution ∝ 19 96 (am) 2 , nearly as good as ∝ 1 6 (am) 2 of the KW action (both values are for r = 1 in d = 4 dimensions). Unlike the KW energy, the BC energy has an imaginary contribution ∝ 1 4 am, which is not a desirable property. What is particularly disconcerting in the free-field dispersion relation of a BC fermion at heavy quark mass is that the global minimum is not necessarily at a| p| = 0. For am = 0.5 the effect happens to be numerically small, but a spontaneous breaking of translation invariance (even if confined to Weiss-type sub-domains of the Brillouin zone) would cause headache.
From eqn. (18) one finds the eigenvalues of the free-field BC operator. The result for r = 1 and am = 0 is shown in Fig. 9. Similar to the KW case, eigenvalues come in γ 5 -pairs and are purely imaginary. The only difference to Fig. 5 is that the BC range is slightly narrower.
It is instructive to repeat this for a series of r values; the result is shown in Fig. 10, with vertical lines marking the abscissa values r = 1/ √ 3, 1/ √ 2 where the number of species changes. The spectral range is seen to increase with growing r. In addition, the low-energy end of the eigenvalue spectrum seems far less stable than in the KW case, the canonical choice r = 1 seems to represent a small island of stability.
Given Fig. 10, one may wonder about the existence of an analytic function which describes the upper end as a function of r. In App. E we derive the free-field spectral bound Borici-Creutz versus staggered, L/a=32, r=1 Figure 11: Sorted free-field eigenvalues (with a 2-fold degeneracy removed) of the BC operator at r = 1 plotted versus sorted staggered eigenvalues (with a 4-fold degeneracy removed). In both cases the imaginary part Im(λ) = λ/i at am = 0 is used, and plenty of degeneracies remain. The dotted line shows the identity for comparison.
smallest non-zero BC eigenvalue is found in essentially the same place 2 as the smallest staggered eigenvalue, see Fig. 11. This amounts to an enhancement of the condition number of D † D, compared to the staggered formulation at the same am, by a factor up to (1 + √ 2) 2 = 5.8284 (in the chiral limit). Hence, the slowdown (relative to the staggered formulation) implied by the larger condition number is not as pronounced as in the KW case, but still significant.

Summary
In this paper we tried to fill some of the most obvious gaps in the knowledge about the two most popular minimally doubled fermion actions, namely the formulations due to Karsten-Wilczek (KW) and Borici-Creutz (BC), respectively. The gaps concern the eigenvalue spectra and the dispersion relations (including the leading cut-off effects on the heavy fermion mass) in the free-field limit. We studied these issues as a function of the lifting parameter r, in order to see how the number of species gets reduced from 16 (at r = 0) to 2 (at r = 1).
Regarding the eigenvalue spectra we find an extension, relative to the staggered/naive one, by a factor 3.5 for KW fermions, or 1 + √ 2 2.4142 for BC fermions (both at r = 1). This leads to an enhancement of the condition number of D † D (as relevant for generating dynamical ensembles) by a factor up to 12.25, or 5.8284, respectively, compared to the staggered/naive case. This, together with the matrix size being a factor 4 larger than for staggered fermions, limits our optimism regarding the computational efficiency of these two formulations.
In addition, we studied the dispersion relations. On the one hand, we find that the KW operator features very well in this respect. It follows the continuum dispersion relation more closely than the Wilson operator. In particular at a p = 0 the cut-off effects on the heavy quark mass start at O((am) 2 ), just like the naive/staggered action, not at O(am) like the Wilson operator. On the other hand, the dispersion relation of the BC operator in d = 4 dimensions shows some more problematic features, including a funny behavior at small a| p| and an imaginary part of the heavy quark rest mass which starts at O(am).
Obviously, there remain many unexplored issues with these fermion formulations. We think it would be interesting to study the behavior of small eigenvalues on interacting backgrounds (especially some with non-zero topological charge), and how they implement the constraints imposed by the Nielsen-Ninomiya theorem. Our Figs. 7 and 11 are inspired by Figs. 6 and 7 of Ref. [26], and we hope to see the "fingerprint property" of low-energy fermion eigenvalues 3 confirmed with the KW and BC formulations, too. Also some more light on the mixing pattern with lower-dimensional operators (beyond what was found in [32][33][34][35][36][37]) might prove useful. Overall, we feel a collaboration aiming for exploratory large-scale production runs with minimally doubled fermions would be well advised to give first priority to the KW formulation.

A Dispersion relations A.1 Naive fermions
The naive operator and its Green's function take the form The dispersion relation follows from searching for zeros of the denominator with p 4 → iE, so means that the physical solution is given by the positive root A. 2

Wilson fermions
The Wilson operator and its Green's function take the form withp µ = 1 a sin(ap µ ) andp µ = 2 a sin( apµ 2 ). It follows that or ar 2p 2 = dr a − r a µ cos(ap µ ), and searching for a zero of the denominator with p 4 → iE yields For r = 1 the identity cosh 2 − sinh 2 = 1 turns this into a linear equation in cosh(aE) which one solves for aE > 0 by means of acosh(x) = ln(x + √ x 2 − 1) for x > 1. For r = 1 one stays with a quadratic equation in cosh(aE) which one addresses by first solving for a real positive cosh(aE) and then inverting the cosh.

A.3 Karsten-Wilczek fermions
The KW operator and its Green's function take the form where in the last step specific properties of the Dirac-Clifford algebra were used. Searching for a zero of the denominator with ar which does not necessarily yield a real solution for E. In such a situation one should go for a complex E, and treat its real part as the "energy" of the respective mode. In other words yields a complex sinh(aE), and through the asinh function the definition of a complex aE is obtained, whose positive real part is plotted against

A.4 Borici-Creutz fermions
The BC operator and its Green's function take the form and our task is to further simplify the denominator. The first term is symmetric inp ρ ↔p σ ; it may be rewritten as 1 2 ρ,σ {γ ρ , γ σ }p ρpσ = λp 2 λ , where the Dirac-Clifford property of the γ-matrices has been used. For exactly the same reason the fourth term may be rewritten as λ , where the Dirac-Clifford property of the γ -matrices has been used. The two cross-terms are a bit more tricky to deal with. It proves useful to notice that the second term can be inflated to look like ar with the substitution p 4 → iE for aE.
which the substitution then brings into the form (with i, j running from 1 to d − 1) In d = 2 space-time dimensions this expression simplifies to (each sum contains a single term) while in d = 4 space-time dimensions one finds In the special case r = 1 the d = 2 version simplifies to These equations look complicated, and this is why we shall work our way backwards, from the simplest case to the more complicated case. A peculiar feature of the d = 2, r = 1 case is that the equation is linear in sinh(aE) and cosh(aE). This suggests multiplying eqn. (48) with exp(aE) to obtain (50) which is a quadratic equation in e aE . Evidently, this means that we should go for the two complex e aE as function of ap 1 , to obtain a complex aE whose positive real part is plotted against |p 1 |. By contrast, the d = 4, r = 1 case has a mixed term in sinh(aE) cosh(aE). The hyperbolic semi-angle substitution t = tanh(aE/2), whereupon sinh and upon multiplying this equation with (1 − t 2 ) 2 one finds the (possibly modified) condition which amounts to a fourth-order polynomial in t.
For generic r we resort to the hyperbolic semi-angle substitution, regardless of the spacetime dimension. For d = 2 we obtain the relation and upon multiplying this equation with (1 − t 2 ) 2 one finds the (possibly modified) condition which amounts to a fourth-order polynomial in t. Note that for r 2 = 1 the second term in eqn. (53) vanishes. It is then sufficient to multiply the equation with 1 − t 2 , and one ends up with a quadratic polynomial in t (equivalent to the procedure used above). In other words, after setting r = 1 and dropping a factor 1 − t 2 eqn. (54) is equivalent to eqn. (50). For d = 4 the same semi-angle substitution yields and upon multiplying this equation with (1 − t 2 ) 2 one finds the (possibly modified) condition which amounts to a fourth-order polynomial in t. Upon setting r = 1 eqn. (56) simplifies to eqn. (52) without further ado.
Using the built-in capabilities of a computer algebra program or a numerical package like matlab/octave, it is straight-forward to find all (in general complex-valued) solutions to a fourth-order polynomial with given numerical coefficients. In this spirit we evaluate, for a given (p 1 , p 2 , p 3 ) configuration, the four solutions t and apply aE = 2 atanh(t) to obtain the energies. The one with the smallest positive real part is interpreted as the energy of the fermion in that momentum configuration, its imaginary part gives the damping of the pertinent mode. This is the numerical basis of all dispersion relations shown in this article. On the analytical side, one may proceed one step further upon expanding the physical solution in powers of am. This yields results relevant to assess the suitability of these actions for heavy-quark physics, as discussed in the main part of the article and App. B.
Hence, the rest-mass of a fermion in the naive discretization has cut-off effects O((am) 2 ).

B.2 Wilson fermions
At a p = 0 the Wilson dispersion relation for arbitrary d and r = 1 simplifies to cosh(aE) = 1 2(1 + am) which is solved if exp(aE) = 1 + am, that is for aE = log(1 + am). The series expansion shows that such cut-off effects scale like O(am). For arbitrary d and generic r one starts from the quadratic equation (r 2 − 1) cosh 2 (aE) − 2r(r + am) cosh(aE) + 1 + (r + am) 2 = 0 whereupon cosh(aE) = r(r + am) ± 1 + 2ram + (am) 2 out of which only the second solution (with negative sign) is physical, since it is the one which agrees, in the limit r → 1, with the solution found in this special case. This yields the expansion which, again, in the special case r = 1 is found to agree with the previous expansion. The lesson is that cut-off effects of Wilson fermions are linear in am. It is impossible to get rid of this undesirable term through a clever choice of r, since for r = 0 we are back to 2 d species.

B.3 Karsten-Wilczek fermions
At a p = 0 the KW dispersion relation simplifies to 0 = − sinh 2 (aE) + (am) 2 and thus to the form (57) of the naive action. Accordingly, the expansion of the rest energy of a static KW fermion in powers of am agrees with (58). Hence, the KW action yields a 2 species formulation which maintains the desirable heavy-quark features of the naive discretization.

B.4 Borici-Creutz fermions
but only the first solution (with positive sign) is physical, since it is the one which agrees, in the limit r → 1 with the solution (64) found previously. It expands as and a quick check reveals that each coefficient in the r = 1 expansion is recovered in that limit. For r = 1 and d = 4 the solution of eqn. (63) can only be given as the logarithm of the roots of the polynomial (ir + r 2 − 1)z 4 + (−2ir − 4r 2 )z 3 + (4m 2 + 6r 2 + 2)z 2 + (2ir − 4r 2 )z − 1 − ir + r 2 = 0, and a power expansion yields which, for r → 1, would indeed simplify to (65). In short, we find that in d = 2 dimensions the rest-mass of a BC fermion has discretization effects O((am) 2 ) for generic r. For r 2 = 4/3 they are even pushed to O((am) 4 ). By contrast, in d = 4 dimensions the rest mass of a BC fermion has O(am) cut-off effects, but this order affects only the imaginary part. Quite generally, it seems that in d = 4 dimensions the real part of E/m is even in r and am, while the imaginary part is odd in r and am.
Another way to see the difference between the cases d = 2 and d = 4 is to apply the hyperbolic semi-angle substitution to eqn. (63). Multiplying it with (1 − t 2 ) 2 yields where t = tanh(aE/2). Specifically for d = 2 the troublesome cubic term is gone and the equation is bi-quadratic, while for d = 4 one ends up with which is a genuine fourth-order equation in t.

C.1 Naive fermions
The denominator of G nai at am = 0 is a 2p2 = µ sin 2 (ap µ ). It has 16 zeros in the Brillouin zone, one at ap µ ∈ {0, π} for each µ, if the range in each direction is taken to be ] − π 2 , 3π 2 ].

C.2 Wilson fermions
The denominator of G W at am = 0 is a 2p2 + ( a 2 r 2p 2 ) 2 ; evidently it is only zero if a 2p2 = 0 and a 2p2 = 0 hold simultaneously. The first term has 16 zeros in the Brillouin zone, the second one only one, at (0, 0, 0, 0). The Wilson term thus lifts 15 of the 16 species of the naive action to a level 2r/a, 4r/a, 6r/a and 8r/a, with degeneracies 4, 6, 4, and 1, respectively.
In view of a similar discussion below for BC fermions, it is perhaps useful to illustrate the solutions to the system (71) in the (r, ap 4 ) plane, see Fig. 12. The degeneracies and multiplicities of the modes are given in the legend and the caption. The main mode (0, 0, 0, 0) is labeled "survivor (1, +)", since it is non-degenerate with correct chirality. The doubler mode (0, 0, 0, ±π) is labeled "survivor (1, −)", since it is non-degenerate with opposite chirality.

C.4 Borici-Creutz fermions
and thus one solution, ap = 0, is independent of r. To the second square bracket we apply the trigonometric semi-angle substitution t = tan(ap/2) with sin(ap) = 2t/(1 + t 2 ) and cos(t) = (1 − t 2 )/(1 + t 2 ). Upon multiplying the result with 1 + t 2 , the second factor becomes and this yields the two-fold zero t = −1/r, hence ap = −2 arctan(1/r). For the non-symmetric modes it is useful to notice that (72) is the sum of two squares and one can thus reformulate the condition as a system of two coupled equations The aforementioned trigonometric semi-angle substitution turns this into which, after multiplication by (1 + t 2 1 )(1 + t 2 2 ), leads to the conditions There are four real solutions, {t 1 = 0, t 2 = 0}, {t 1 = −1/r, t 2 = −1/r}, where the first two are again symmetric in p 1 ↔ p 2 , and the last two interchange under t 1 ↔ t 2 .
For the square-root in (79) to be real, one needs 1 − 2r 2 − 3r 4 ≥ 0, and this means r 2 ≤ 1/3. At r = 1/ √ 3 the solutions become t 1 = t 2 = − √ 3, and thus coincide with the symmetric solution, t = −1/r = − √ 3 at this point. In short, for 0 < r < 1/ √ 3 we have a 4 species action (with two symmetric and two non-symmetric modes), while for 1/ √ 3 < r the BC action in d = 2 dimensions encodes for 2 species (which live on the diagonal of the Brillouin zone).
The second version is useful for solutions with 2-to-2 momentum pairing. Without loss of generality we assume p 1 = p 2 ≡ p, p 3 = p 4 ≡ q, so eqn. (84) takes the form 0 = sin(ap) + r{1 − cos(aq)} 0 = sin(aq) + r{1 − cos(ap)} (93) but this is identical to the system (76) for BC fermions in d = 2 dimensions. It follows that {t = 0, u = 0} and {t = −1/r, u = −1/r} are the symmetric solutions, and are the non-symmetric ones. Evidently, the second solution emerges from the first one by interchanging t ↔ u. The square-root is real for |r| ≤ 1/ √ 3, and at this point the nonsymmetric solutions take the form t = −1/r = − √ 3, u = −1/r = − √ 3, which means that they merge into the symmetric solution. We recall that the choice p 1 = p 2 ≡ p was one out of three possibilities, hence we have six rather than two non-symmetric solutions.
The first version of the system was not used at all. It seems (83) would be most useful for finding a totally unsymmetric mode, i.e. one with pairwise unequal p 1 , p 2 , p 3 , p 4 . Apart from having already found 2 + 8 + 6 = 16 solutions (for r small enough), permutations demanded by invariance under exchange of any axes would beef up a totally unsymmetric solution to 4! = 24 solutions, and that is too many of them.
Overall, we thus arrive at the following picture for BC fermions in d = 4 dimensions. For infinitesimally small r there is an invariant solution, ap = (0, 0, 0, 0), a symmetric solution, ap = − arctan(1/r)(2, 2, 2, 2), eight solutions with 3-to-1 momentum pairing of the type (90), and six solutions with 2-to-2 momentum pairing of the type (94). The first two solutions have the correct chirality on topologically charged backgrounds, the 3-to-1 paired solutions all have opposite chirality, and the 2-to-2 paired solutions have correct chiralities again. At r = 1/ √ 3 a dramatic merger and exchange of chiralities takes place, since the 2-to-2 paired solutions cease to exist, but their chiralities are transferred to the remaining modes in the sense that the 3-to-1 paired solutions fall into two sub-categories (the four of them who passed through the central point now have correct chirality), and the symmetric solution gets flipped to opposite chirality at this point. In short, for r slightly above 1/ √ 3 the BC action has 10 species (five of each chirality). At r = 1/ √ 2 the second change takes place, since the eight solutions with 3-to-1 pairing annihilate each other (in two different points of the Brillouin zone, and they can do so, since four of them have correct chirality, and four of them have opposite chirality). Slightly above this value of r the BC action is minimally doubled, i.e. has one species of each chirality.
Following a similar attempt in the KW case, we try to illustrate the various modes in the (r, ap) or (r, aq) plane in Fig. 15. At r = 1/ √ 3 the dramatic merger and exchange of chiralities takes place, as discussed above. At r = 1/ √ 2 the second reduction in the number of species takes place, since at this point all 3-to-1 paired solutions annihilate each other. The trivial solution (0, 0, 0, 0) has correct chirality, symmetric solution has correct chirality for r < 1/ √ 3 and opposite chirality for r > 1/ Here p is the three-fold momentum, and q is the single momentum. The pertinent contours, with momentum range ] − 5 4 π, 3 4 π[ for both ap and aq, are shown in Fig. 16. For small r one sees the trivial solution t = u = 0, the symmetric solution t = u = −r, as well as the upper line of (90). At r = 1/ √ 3 the pole above the diagonal hits the symmetric solution, and the dramatic exchange of chiralities (which involves solutions which are not visualized in Fig. 16) takes place. And at r = 1/ √ 2 the two poles below the diagonal annihilate each other, and after this point the BC action is minimally doubled.
In d = 2 dimensions, eqn. (44) simplifies to and with p 1 = (q + iE)/ √ 2 and p 2 = (−q + iE)/ √ 2 it takes the form which is a quartic equation in e aE/ √ 2 . Specifically at q = 0 it simplifies to and upon setting r = 1 it further simplifies to This equation has formally four solutions which, in the limit r → 1, is seen to coincide with the previous expansion.
In short, for BC fermions with propagation in the hyperdiagonal direction we find similar properties than with the standard propagation along the d-th axis. In d = 2 and d = 4 dimensions the rest mass of a BC fermion with diagonal propagation direction has O(am) cutoff effects, but this order affects only the imaginary part. The coefficient of the O((am) 2 ) cut-off effects is −[3r 2 + 1]/12 in d = 2 dimensions and −[3r 2 + 1]/24 in d = 4 dimensions. Again, it seems that the real part of E/m is even in r and am, while the imaginary part is odd in r and am. Overall, we do not see any compelling advantage of the (++) or (++++) propagation direction over the standard propagation in the d-th direction. A peculiarity of d = 2 space-time dimensions is that the propagation direction can be chosen orthogonal to the hyperdiagonal direction, and in this case the O(am) cut-off effects disappear, and for r 2 = 4/3 the leading cut-off effect in the heavy-quark mass is pushed to O((am) 4 ).

E Spectral bounds E.1 Karsten-Wilczek operator
Plugging in the momenta in eqn. (10) at m = 0 yields and with ω KW ≡ max(λ KW /i) it follows that the symmetry among the spatial axes implies for some appropriately chosen momentum configuration. Obviously, setting ap d = π/2 helps to reach the maximum. Using 2 sin 2 (ap/2) = 1 − cos(ap) we thus need to maximize over p ∈ [−π, π]. We need to keep in mind that at either endpoint we have the value and equating this with the endpoint value shows that the switching beween the two solutions happens at r = 1/3. Plugging the d = 2 result into the general expression yields and equality with the endpoint value is reached at r = 1/2. To summarize, in d = 4 dimensions we find the spectral bound (11). In d = 2 dimensions |Im(λ KW )| ≤ 2/(1 − r) r ≤ 1/2 1 + 2r r ≥ 1/2 (123) and the bound for large r generalizes to 1 + 2(d − 1)r in d dimensions. For r → 0 the general result tends to √ d, which is the upper bound of the staggered free-field eigenvalue spectrum.