New One-Flavor Hybrid Monte Carlo Simulation Method for Lattice Fermions with gamma-five Hermiticity

We propose a new method for Hybrid Monte Carlo (HMC) simulations with odd numbers of dynamical fermions on the lattice. It employs a different approach from polynomial or rational HMC. In this method, gamma-five hermiticity of the lattice Dirac operators is crucial and it can be applied to Wilson, domain-wall, and overlap fermions. We compare HMC simulations with two degenerate flavors and (1 + 1) degenerate flavors using optimal domain-wall fermions. The ratio of the efficiency, (number of accepted trajectories) / (simulation time), is about 3:2. The relation between pseudofermion action of chirally symmetric lattice fermions in four-dimensional(overlap) and five-dimensional(domain-wall) representation are also analyzed.


Introduction
In hybrid Monte Carlo simulations [1], the positive-definiteness of the action is essential to consider it as the statistical weight. When a lattice Dirac operator D is given, a positive-definite action of two degenerate flavors is easily constructed by using the hermitian conjugate of D, namely D † D. The major difference of two flavor simulations and (2 + 1) flavor simulations is that one cannot easily write down the pseudofermion action for the one flavor sector. Rational or polynomial HMC methods [2,3], which approximates the square root of D † D, are mostly used for odd flavor simulations.
In this paper, we provide a pseudofermion action for the one flavor sector of the lattice fermions with γ 5 hermiticity without invoking the square root approximation for D † D. The main idea is very simple. For any lattice Dirac operators D with γ 5 hermiticity, P + DP + and P − (1/D)P − are hermitian, and one can construct a one-flavor pseudofermion action using these. The resultant action has the same determinant as D without any approximations. The non-trivial parts are to first check the positive-definiteness of the pseudofermion action and the discussion of how to obtain the pseudofermion action when there are mass preconditioners like the one in Hasenbusch method.
The construction of the paper is as follows. In Sec. 2 -5, the application to Wilson fermions, Wilson fermions with mass preconditioner, domain-wall type fermions and overlap fermions are demonstrated, respectively. In Sec. 6, the relation between the pseudofermion action in four and five dimensional representations are presented. Numerical results are given in Sec. 7, while a summary and conclusion are provided in Sec. 8.

Wilson Fermions
In this section, we derive a pseudofermion action which will yield the same determinant as the Wilson-Dirac operator, where 1 2×2 is two-by-two unit matrix and σ i (i = 1, 2, 3) are Pauli matrices. Throughout this work, we use chiral basis for gamma matrices. To get positive-definite pseudofermion action, we use a determinant relation which is valid for general matrix, where P ′ +/− are the projectors which reduce two chiral components into one chiral sector, P ′ + = (1 0), P ′ − = (0 1). The operators (P ′ − D W (m) −1 P ′ † − ) and (P ′ + D W (m)P ′ † + ) are hermitian due γ 5 symmetry. When the inverse of (W + m) is well-defined, the determinant of (1) is written as such that W H (m) is the Schur complement of D W (m), i.e., Thus, the pseudofermion action for one-flavor Wilson fermions can be written as where H W (m) = γ 5 D W (m), Φ 1 is a pseudofermion field without a Dirac index, and Φ 2 is a pseudofermion field with two spinor components. The positive-definiteness of the operator W H (m) are discussed as follows. We note that for any background gauge field, the eigenvalues of W and (σ · t) satisfy the inequalities 1 : 0 ≤ λ(W) ≤ 8, and |λ(σ · t)| ≤ 4. It then follows that W H (m) is positive-definite for m > 4. Now, consider decreasing m from 4 to a smaller value. If the operator W H (m) is not positive for some m, e.g. m ′ , then there must exist values of m at which det(W H (m)) is zero or singular in the region m ′ < m < 4. Among these values, if we denote the largest one as m cr , then for m > m cr , the positive-definiteness of W H (m) is assured.
The value m cr corresponds to the opposite sign of the smallest eigenvalue of W. To see this, we use a relation 2 , From this relation, one can see that the value of m which makes det (W H (m)) singular (i.e. det(W + m) = 0) is larger than the value of m which satisfies det (D W (m)) = det (W H (m)) = 0. There are no values of m which satisfy det (W H (m)) = 0 or det (W H (m)) = ±∞ above that. The condition for the positive-definiteness of W H (m) can then be written as The smallest mass one can use in the one flavor method here is restricted by this bound. But it can be relaxed by using Hasenbusch's preconditioner as discussed in the next section.

Generating Pseudo Fermion Fields
To generate the pseudofermion field Φ 2 from a Gaussian random noise field Ξ 2 , we need to take the square root of W H (m), i.e., Φ 2 = √ W H (m)Ξ 2 . This can be approximated by using the rational function of W H (m), 1 Throughout this work, λ(X) means any one of the eigenvalues of an operator X. 2 The proof of this relation is given in the appendix.

2
where p 0 , p l and q l are expressed in terms of Jacobian elliptic functions. At first glance, the operations in (8) look formidable. However, since W H (m) is the Schur complement of D W (m), the inversion of the operator within the summation, in Eq. (8) can be obtained by the inversions of the operator (D W (m) + q l P − ), i.e., Note that one cannot apply the multi-shift solver in Eq. (9), since P − does not commute with D W (m). However, for a given N app number of inversions, the total number of iterations in the solver can be reduced by using the same idea as the chronological inversion method [4]. When one solves a set of linear equations (D W (m) + q l P + )η l = Ξ 2 (1 ≤ l ≤ N app ) with a given Ξ 2 serially from smaller l by using the iterative method, one can set a better initial guess η (0) l for the iterative method to solve (D W (m) + P + )η l = Ξ 2 (2 ≤ l) by using a linear combination of the solutions η j ( j < l) which have already been calculated i.e. η (0) l = j<l c l j η j . The coefficients c l j are determined according to the prescription written in Ref. [4].
In this one flavor method, generating pseudofermions using the approximation given by Eq. (8) makes the simulation not exact. However, without using higher degrees of the approximation, one can make the algorithm exact by adding an accept/reject step after generating the pseudofermion field. In the exact algorithm, the pseudofermion field should be produced according to a probability distribution proportional to e −φ † 1 W H (m) φ . This is obtained by multiplying the operator √ W H (m) to a Gaussian noise field Ξ. In practice, however, the operator one uses in the simulation is approximated by f app (W H (m)) rather than √ W H (m). This leads to a probability distribution that is pro- To adjust the difference, one can add an accept/reject step for φ with the probability This factor should be smaller than 1 and can be enforced by choosing f app (x) > √ x for the whole eigenvalue region of W H (m). The discussion using eigenvectors is given in the appendix.

Wilson fermions with the Hasenbusch method
The idea of the Hasenbusch method [5] is to factorize the determinant det (D W (m 1 )) into a product of determinants, det D W (m 1 ) = det (D W (m 1 )/D W (m 2 )) · det (D W (m 2 )) (10) and the pseudofermion force coming from det (D W (m 1 )/D W (m 2 )) is updated less frequently than the one coming from det (D W (m 2 )) in the molecular dynamics steps. The parameter m 2 is chosen such that the simulation cost is reduced. Furthermore, as already mentioned, in the one flavor method for Wilson fermions presented here, the factorization of the determinant allows us to use a smaller fermion mass in the simulation. For the second factor det(D W (m 2 )), the same argument as the previous section Sec. (2) can be applied. Here, we assume that m 1 < m 2 and that m 2 is large enough such that the positive-definiteness of W H (m 2 ) is always assured. Now, we consider how to treat the first factor det(D W (m 1 )/D W (m 2 )). Naively, one might consider the pseudofermion action like, To remedy this situation, we rewrite the determinant det(D W (m 1 ))/(D W (m 2 )) by using an operator which has different masses for the chirality plus and chirality minus sector, i.e., One may note that the first and the second factor on the right hand side are the same as the determinant of, 3 respectively. Then, the pseudofermion action is written as where Φ 3 and Φ 4 are pseudofermion fields with two spinor components. The condition for the positive-definiteness is given as follows. When m 1 = m 2 , the operators in Eq (13) are the identity operators and their positive-definiteness is trivial. When one decreases m 1 with m 2 fixed, the positivedefiniteness will be lost only if either of the determinants in Eq. (12) becomes zero or singular. Thus, up to the largest m 1 which satisfies det(D W (0) + m 1 P + + m 2 P − ) = 0, the positive-definiteness of the operators in Eq (13) are assured. Note that this mass value is smaller than the limit given in Eq (7). This can be seen by using the eigenvalue relation which is similar to Eq (6), The largest m 1 which yields det(D W (0)+m 1 P + +m 2 P − ) = 0 is larger than the largest m which yields det D W (m) = 0. This can be understood from and the properties of spectral flows of H W (m) [6]. Moreover, for some gauge configurations, det(D W (m)) may not be zero for any value of m. But there must exist m 1 which satisfies det(D W (0) + m 1 P + + m 2 P − ) = 0 for any gauge configurations.

Domain-wall fermions
The domain-wall type fermion operator [7] can be expressed as, with L(m) = P + L + (m) + P − L − (m) such that, where m is the (bare) fermion mass, and m 0 ∈ (0, 2) is a parameter called the "domain-wall height". The label s is the coordinate in the fifth dimension. Throughout this work, we assume that the number of sites in the fifth dimension N s is even. The constant c and the diagonal matrix ω = diag{ω s } specify the type of domain-wall fermion. The operator D dwf (m) is the conventional domain-wall fermion for c = 0 and ω s = 1. It becomes an optimal domain-wall fermion when c = 1 and ω s 's are tuned such that maximal chiral symmetry is obtained.
To obtain the pseudofermion action for one flavor, we modify the operator D dwf (m) by multiplying 1/(1 + cL) from the right. with In the following, we suppress the argument c of M ± (m, c) for simplicity. Using the Schur decomposition, the determinant of the operator D ′ dwf (m) can be written as where Here R 5 is the reflection operator in the fifth dimension, (R 5 ) s,s ′ = δ s,N s +s ′ −1 , which is introduced such that W H (m) and W H (m) are hermitian. For optimal domain-wall fermions, one still can choose ω s 's which maintain maximal chiral symmetry and satisfies ω N s +1−s = ω s . After incorporating the contributions of the Pauli-Villars fields, the fermion determinant for domain-wall fermions becomes, In principle, one can use a pseudofermion action like, But it is known that in the two-flavor simulation of domain-wall fermions, it is effective when a pseudofermion action uses a single set of pseudofermion field to estimate both the light fermion and Pauli-Villars terms [8]. Then, we use the same we used in Hasenbusch method for Wilson fermions.
Using the Schur decomposition of [D dwf (1) − (M + (1) − M + (m)) P + ], we obtain the relation where The properties of these matrices are given in the appendix. Using (26), we can write the inverse of (24) as and it can be used to construct the pseudofermion action. Using (B.6), we can simplify (28) to where, The constant g ′ and vector v are given in Eq. (B.6) and (B.4). In the following, the arguments of g ′ are suppressed for simplicity (g ′ = g ′ (m, 1, c)). Note that the five-dimensional matrix in Eq. (28) is reduced to a four-dimensional matrix in this expression. Thus we can write the pseudofermion action for one-flavor domain-wall fermions as where Φ 1 and Φ 2 are the pseudofermion fields (on the four-dimensional lattice) with two spinor components. Now we assert that the operators in (29) are positive-definite for 0 < m ≤ 1. At m = 1, they are equal to the identity operator, and thus are positive-definite. As m is decreased, the operators in (29) will cease to be positive-definite only

Generating Pseudo Fermion Field
We now discuss how to approximate the inverse square root of the operators A and B when one generates the pseudofermion fields Φ 1 and Φ 2 from the Gaussian noise. Focusing on Eq. (35), we start from the inverse relation of five-dimensional operators, Multiplying by S † from the left and S from the right, one obtains Here, M ss ′ = 0 for s < s ′ . Then, the relation, Eq. (34), also holds for the sub-block (s, s ′ ) = (N s , N s ), One can obtain a similar equation for W H (m) and ∆ − . For the square root of a general positive-definite operator A, the rational approximation can be used.
In the case of A, one has to calculate Each term of the summation in Eq. (36) is 3 Recall that M(m) for a conventional domain-wall fermion is a difference operator with anti-periodic boundary condition. 6 Then, to obtain φ 1 one needs to calculate The same method can be used for the operator B.

Overlap Fermions
In this section, we construct the pseudofermion action for one flavor overlap fermions [9], This operator satisfies Ginsparg-Wilson Relation(GWR) [10], The overlap fermion also possesses γ 5 hermiticity, and so the application is similar to Wilson fermions. We begin by breaking this operator into its chiral components.
The determinant is written as, By using GWR, one can show that the operators on the right hand side are positive-definite for 0 < m ≤ 1, provided that D ov (m) itself is well defined. The pseudofermion action is written as Practically, one has to use a reflection/refraction [11] or topology fixing term [12] to treat or avoid singularities in D ov (m) related to the topological change. The authors of [13] proposed a simulation method for one-flavor with overlap fermions. We now highlight the differences between their work and our work. In Ref. [13], the authors use GWR, and factorize det D ov (m) † D ov (m) into two parts. det The difference of the first factor and second factor of the right-hand side comes from the topological zero-mode of D ov (m). Then, the one-flavor determinant is written as Here, N +/− are the numbers of topological zero-modes with a definite chirality. By using GWR, one obtains the relation between these operators and Schur complement, In other words, in Eq. (46), the determinant is factorized as, Later, we show that the second factor corresponds to the A term in Eq. (31). Due to the cancellation between the numerator and the denominator, the force for the A term in HMC is smaller than that of the B term. When one performs a HMC simulation with a topological fixing term, it is apparent that using the pseudofermion action with the factorization in Eq. (49) is more effective than the factorization in Eq. (51). For the five-dimensional representation Eq. (31), using only the B term and choosing the fermion mass parameter m ′ to satisfy (1 − m ′ ) = (1 − m 2 )/m, one can perform the HMC simulation respecting the factorization in Eq. (49).

The relation of the pseudofermion action with four and five dimensional representations
The ratio of determinants, Eq. (24), is equivalent to the determinant of the effective four-dimensional operator, where the function f (x) is polar or a rational function which approximates the sign function. The form of the function f is determined by ω. The operator H kernel is defined as, By using the Schur decomposition, the one flavor pseudofermion action is written as Examining two cases, c = 0 and c = 1, we show that the five-dimensional expression Eq. (31) reproduces Eq. (54). This equivalence is not only of theoretical value, but aids in practical simulations. When one generates φ 1 and φ 2 , from Eq. (31), one has to know the eigenvalue spectrum of A and B, a priori, in order to apply the rational approximation for the square root function. However, in the four-dimensional case, the spectrum is already known.
In this case, the propagator of the five-dimensional fermion at the boundary s = 1 or N s yields the propagator of the four-dimension fermion [14], 1 Here, B is defined as This relation holds only for c = 0. Note that this operator satisfies chiral symmetry D ch γ 5 + γ 5 D ch = 0 in the limit N s → ∞, in which f (x) becomes the sign function. From Eq. (55) and Eq. (56), After incorporating the pseudofermion field, the left-hand side becomes the first term of Eq. (54). The right-hand side equals to the A term of Eq. (31) by substituting g ′ = (1 − m) and v † = (0, · · · , 0, 1). Iter are the average CG iterations for one trajectory (for generating initial pseudofermion fields, molecular dynamics, and their sum respectively).
In this case, there are relations for D 4d and five-dimensional operators [15], By using the relations , one can show that Eq. (60), Eq. (61) and projector P +/− reproduce the operators in Eq. (31).

Numerical tests
We compare the efficiency of the HMC simulation for two-flavor and (1 + 1)-flavor QCD with domain-wall type fermion with c = 1 and ω s = 1, on a 12 3 × 24 × 16(N s ) lattice. For the gluon action, we use Iwasaki gauge action at β = 2.30. In the molecular dynamics, we use the Omelyan integrator [16], and the Sexton-Weingarten method [17]. The pseudofermion action for the two-flavor simulation is the one with even-odd preconditioning which is described in Ref. [18]. The time step for the gauge field, (∆τ Gauge ), is the same for both two-flavor and (1 + 1)-flavor cases, while the time step (∆τ PF ) for the pseudofermion fields in the (1 + 1)-flavor case is four times larger than that for the two-flavor case. The acceptance rate is roughly the same for both cases. We use conjugate gradient (CG) with mixed precision for the inversion of the quark matrix (with even-odd preconditioning). The length of each trajectory is set to two. After discarding 300 trajectories for thermalization, we accumulate 100 trajectories for the comparison of efficiency. Our results are given in Table.1. We see that the acceptance rate is almost the same for (1 + 1)-flavor and two-flavor simulations. If the auto-correlation time is the same, then the efficiency of HMC can be estimated by the total acceptance divided by the CG iteration number, and the efficiency ratio for two-flavor and (1 + 1)-flavor is about 3 : 2.

Concluding remarks
In this work, we presented one-flavor method for HMC. The largest difference of the method presented here from RHMC and PHMC is that the pseudofermion action yields the one-flavor determinant without any approximations. For the overlap fermion, the difference from the method in Ref [13] is that this can be used even if GWR is not exact. For the lattice fermions with γ 5 symmetry, it is always possible to construct real (hermitian) pseudofermion action, but one has to careful about the positive-definiteness, since it depends on the type of lattice fermions used. For chirally symmetric fermions like domain-wall/overlap fermions, positive-definiteness is assured in the entire mass parameter region which can be used in two-flavor simulations. On the other hand, for Wilson fermions, there exist a bound for the smallest mass where the positive-definiteness is assured. If the bound is larger than the mass value which is intended to use, one has to add Hasenbusch mass preconditioner then one can push the mass smaller. Comparison between two-flavor and (1 + 1)-flavor using domain-wall type fermion HMC simulation shows that one can increase the step size of (1 + 1)-flavor simulation while keeping same acceptance ratio. The reason might be that for (1 + 1)-simulation it is effectively the same as using four-dimensional operator and the force from the bulk mode is completely cancelled between the light fermion and Pauli-Villars field, while for the two-flavor pseudofermion action used in the comparison in this work, the cancellation was done only partly. The bottle neck of this method is generating pseudofermion field from Gaussian noise and it should be tuned and improved. Not only HMC itself, one can use these pseudofermion action for reweighting method, for example, to adjust the strange quark mass or to see the effects due to the difference of up down sea quark mass by using existing configurations. This approach for one-flavor simulation should be investigated further numerically and theoretically.

Appendix C. Exactness -Point of View with Eigenvalues
Consider the eigenvalues and eigenvectors of the operator W H .
The Gaussian noise field Ξ, and pseudo fermion field φ, are expressed as linear combination of the eigenvectors ψ i .
The path integral of φ is written as the product of integral w.r.t the coefficient c i .

Appendix D. Linear Algebra
Here, we remind the readers some relations of linear algebra.