Pseudoscalar Meson in Two Flavors QCD with the Optimal Domain-Wall Fermion

We perform hybrid Monte Carlo (HMC) simulatons of two flavors QCD with the optimal domain-wall fermion (ODWF) on the 16^3 x 32 lattice (with lattice spacing a ~ 0.1 fm), for eight sea-quark masses corresponding to pion masses in the range 228-565 MeV. We calculate the mass and the decay constant of the pseudoscalar meson, and compare our data with the chiral perturbation theory (ChPT). We find that our data is in good agreement with the sea-quark mass dependence predicted by the next-to-leading order (NLO) ChPT, and provides a determination of the low-energy constants \bar{l}_3 and \bar{l}_4, the pion decay constant, the chiral condensate, and the average up and down quark mass.

Lattice QCD with exact chiral symmetry [1,2] is an ideal theoretical framework to study the nonperturbative physics from the first principles of QCD. However, it is rather nontrivial to perform Monte Carlo simulation such that the chiral symmetry is preserved at a high precision and all topological sectors are sampled ergodically.
Since 2009, we have been using a GPU cluster (currently constituting of 250 Nvidia GPUs) to simulate unquenched lattice QCD with the optimal domain-wall fermion (ODWF) [3,4]. Mathematically, ODWF is a theoretical framework which preserves the chiral symmetry optimally with a set of analytical weights, {ω s , s = 1, · · · , N s }, one for each layer in the fifth dimension [3]. Thus the artifacts due to the chiral symmetry breaking with finite N s can be reduced to the minimum, especially in the chiral regime. The 4-dimensional effective Dirac operator of massless ODWF is which is exactly equal to the Zolotarev optimal rational approximation of the overlap Dirac operator. That is, is the optimal rational approximation of (H 2 w ) −1/2 [5,6]. Recently we have demonstrated that it is feasible to perform a large-scale unquenched QCD simulation which not only preserves the chiral symmetry to a good precision, but also samples all topological sectors ergodically [7]. To recap, we perform HMC simulations of 2 flavors QCD on a 16 3 × 32 lattice, with ODWF at N s = 16 and plaquette gauge action at β = 5.95. Then we compute the low-lying eigenmodes of the overlap Dirac operator, and use its index to obtain the topological charge of each gauge configuration, and from which we compute the topological susceptibility for 8 sea-quark masses, each of 300 configurations. Our result of the topological susceptibility agrees with the sea-quark mass dependence predicted by the NLO ChPT [8], and provides the first determination of both the pion decay constant and the chiral condensate simultaneously from the topological susceptibility. In this paper, we use the same set of gauge configurations to compute the mass and the decay constant of the pseudoscalar meson, and compare our results with the NLO ChPT. We find that our results are in good agreement with the sea-quark mass dependence predicted by NLO ChPT [9], and from which we obtain the low-energy constants F , Σ,l 3 andl 4 . With the low-energy constants, we determine the average up and down quark mass m MS ud (2 GeV), and the chiral condensate Σ MS (2 GeV).
First, we outline our HMC simulation of 2 flavors QCD with ODWF. Starting from the ODWF action S =ΨDΨ [3] on the 5D lattice, we separate the even and the odd sites (the so-called even-odd preconditioning) on the 4D lattice, and rewrite D as where m q denotes the bare quark mass, D w denotes the standard Wilson Dirac operator plus a negative parameter −m 0 (Here m 0 = 1.3 in this paper.), and D EO/OE w denotes the part of D w with gauge links pointing from odd/even sites to even/odd sites, and (ω) ss ′ = ω s δ ss ′ , Since det D = det S −1 1 · det C · det S −1 2 , and S 1 and S 2 do not depend on the gauge field, we can just use C for the HMC simulation. After including the Pauli-Villars fields (with m q = 2m 0 ), the pseudo-fermion action for 2 flavors QCD (m u = m d ) can be written as In the HMC simulation [10], we first generate random noise vector ξ with Gaussian distribution, then we obtain φ = C −1 P V Cξ using the conjugate gradient (CG). With fixed φ, the system is evolved under a fictituous Hamiltonian dynamics, the so-called molecular dynamics (MD). In the MD, we use the Omelyan integrator [11], and the Sexton-Weingarten multiple-time scale method [12]. The most time-consuming part in the MD is to compute the vector η = (CC † ) −1 C P V φ with CG, which is required for the evaluation of the fermion force in the equation of motion for the conjugate momentum of the gauge field. Here we take advantage of the remarkable floating-point capability of the Nvidia GPU, and perform the CG with mixed precision [13]. Moreover, the computations of the gauge force and the fermion force, and the update of the gauge field are also ported to the GPU. In other words, almost the entire HMC simulation is performed within a single GPU.
Furthermore, we introduce an auxillary heavy fermion field with mass m H (m q ≪ m H < 2m 0 ), similar to the case of the Wilson fermion [14]. For two flavors QCD, the pseudofermion action (with C H ≡ C(m H )) becomes, which gives exactly the same fermion determinant of (1). Nevertheless, the presence of the heavy fermion field plays a crucial role in reducing the light fermion force and its fluctuation, thus diminishes the change of the Hamiltonian in the MD trajactory, and enhances the acceptance rate. A detailed description of our HMC simulations will be presented in a forthcoming paper [15]. Simulations are carried out for two flavors QCD on a 16 3 × 32 lattice, for eight sea-quark masses m q a = 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, and 0.08 respectively. For the quark part, we use ODWF with N s = 16, and λ min /λ max = 0.02/6.4. For the gluon part, we use the plaquette action at β = 5.95 and β = 5.90 respectively. For the β = 5.95 ensemble, our result of the topological susceptibility has been presented in Ref. [7]. In this paper, we focus on the mass and the decay constant of the pseudoscalar meson for the β = 5.95 ensemble.
For each sea-quark mass, we perform simulations on 30 GPUs independently, with each GPU generating 400 trajectories (Currently, the statistics has been increasing to 500-600 trajectories.) After discarding 300 trajectories for thermalization, each GPU yields 100 trajectories. Thus, with 30 GPUs running independently, we accumulated total 3000 trajectories for each sea-quark mass. From the saturation of the binning error of the plaquette, as well as the evolution of the topological charge, we estimate the autocorrelation time to be around 10 trajectories. Thus we sample one configuration every 10 trajectories, then we have 300 configurations for each seaquark mass. With a GPU cluster of 250 GPUs, we can simulate 8 sea-quark masses concurrently. It takes about 5 months to complete the simulations for the β = 5.95 ensemble.
We determine the lattice spacing by heavy quark potential with Sommer parameter r 0 = 0.49 fm. The lattice spacing versus the quark mass is plotted in Fig. 1. Using the linear fit, we obtain the lattice spacing in the chiral limit, a = 0.1032(2) fm, which gives a −1 = 1.911(4)(6) GeV, where the systematic error is estimated with the uncertainty of r 0 . For the computation of the valence quark propagator of the 4D effective Dirac operator of ODWF, two methods have been presented in Ref. [16]. For ODWF with two auxiliary layers at s = 0 and s = N s + 1 (with ω 0 = ω Ns+1 = 0), and the quark fields defined in terms of these two auxiliary layers, then the valence quark propagator can be written as where D ′−1 x,s;y,s ′ is the propagator of the ODWF operator, which can be obtained by conjugate gradient (CG) with even-odd preconditioning. On the other hand, if one does not introduce the two auxiliary layers at s = 0 and s = N s + 1 with ω 0 = ω Ns+1 = 0, then one can obtain the valence quark propagator by solving the following linear system (with even-odd preconditioned CG), where B −1 x,s;x ′ ,s ′ = δ x,x ′ (P − δ s,s ′ +P + δ s+1,s ′ ) with periodic boundary conditions in the fifth dimension. Then the solution of (2) gives the valence quark propagator In this paper, we compute the valence quark propagator with the point source at the origin, and with parameters exactly the same as those of the sea-quarks.
To measure the chiral symmetry breaking due to finite N s , we compute the residual mass where J 5 (x) =ψ x,n P − ψ x,n+1 −ψ x,n+1 P + ψ x,n , n ≡ N s /2, denotes the chiral density at the central layer in the 5-th dimension, (D c + m q ) −1 denotes the valence quark propagator with m q equal to the sea-quark mass, tr denotes the trace running over the color and Dirac indices, and the subscript {U } denotes averaging over an ensemble of gauge configurations. It turns out that, after averaging over an ensemble of a few hundreds of independent gauge configurations, m res seems to be insensitive to the location of the origin x µ = (0, 0, 0, 0). Thus (3) gives a reliable estimate of the chiral symmetry breaking due to finite N s . The derivation of (3) will be given in a forthcoming paper [17]. In Fig. 2, we plot the residual mass versus the sea quark mass. Using the linear fit, we obtain the residual mass in the chiral limit, m res a = 0.00040(4), less than 5% of the lightest sea quark mass. Thus we confirm that the chiral symmetry is preserved to a good precision in our simulation. In the following, it is understood that each bare sea-quark mass m q is corrected by its residual mass, i.e., m q → m q + m res . Using the valence quark propgator with quark mass equal to the sea-quark mass, we compute the timecorrelation function of the pseudoscalar interpolator where the trace runs over the Dirac and color space. Then the ensemble average C π (t) of each m q is fitted to the formula to extract the pion mass M π a and the decay constant where the excited states have been neglected in (4). In Fig. 3, we plot M 2 π /m q and F π versus m q respectively. Here we have made the correction for the finite volume effect using the estimate within ChPT calculated up to O(M 4 π /(4πF π ) 4 ) [18], since our simulation is done on a finite volume lattice with M π L ∼ 2.0 for the lightest sea quark, and its finite volume effect cannot be neglected.
Taking into account of the correlation between M 2 π /m q and F π for the same sea-quark mass, we fit our data to the formulas of NLO ChPT [9] where Λ 3 and Λ 4 are related to the low energy constants l 3 andl 4 as follows.
The strategy of our data fitting is to search for the values of the parameters Σ, F , Λ 3 and Λ 4 such that they minimize where C i is the 2×2 covariance matrix for M 2 π /m q and F π with the same sea-quark mass, and the matrix elements of C i are estimated using the binning method followed by the jackknife.
where the systematic errors are estimated by varying the number of data points from 8 (M π ≤ 580 MeV) to 5 (M π ≤ 480 MeV).
For completeness, we re-plot everything in Fig. 3 as a function of M 2 π , as shown in Fig. 4.  With the fitted parameters, we can use the NLO ChPT formula (5) to solve for the physical bare quark mass, with the physical inputs M π ≃ 0.135 GeV. On the other hand, we can also use (6) to solve for the physical bare quark mass, with the physical input F π = 0.093 GeV. However, the solutions of these two equations are different. Thus, instead of adopting either one of these solutions, we use the physical ratio  (6) and (5), we obtain the pion decay constant and the pion mass at the physical point, Since we have used the physical ratio 1.45 as the input, in principle, we can only regard either (8) or (9) as our predicted physical result. In order to convert the chiral condensate Σ and the average m u and m d to those in the MS scheme, we calculate the renormalization factor Z MS s (2 GeV) using the non-perturbative renormalization technique through the RI/MOM scheme [19], and our result is [20] Z MS s (2 GeV) = 1.244(18)(39).
Then the values of Σ and the average of m u and m d are transcribed to m MS ud (2 GeV) = 4.06(10)(12) MeV. (11) Since our calculation is done at a single lattice spacing, the discretization error cannot be quantified reliably, but we do not expect much larger error because our lattice action is free from O(a) discretization effects.
We also investigated to what extent our results of the low-energy constants depending on the chiral symmetry of the valence quark propagators. We repeated above analysis with valence quark propagators computed with N s = 32 and λ min /λ max = 0.01/6.4, which has the residual mass m res a = 0.000191 (12) in the chiral limit. The low-energy constants turn out to be in good agreement with those in (7).
Moreover, our present results of the chiral condensate (10) and the pion decay constant (8) are in good agreement with our recent results extracted from the topological susceptibility [7], To compare our results with those obtained by other lattice groups, we rely on the recent review [21]. In general, our results of the SU (2) low-energy constants, the chiral condensate and the average up and down quark mass are compatible with those obtained by other lattice groups.
Currently, for the 16 3 ×32×16 lattice, we are increasing the statistics of both β = 5.95 and β = 5.90 ensembles, targeting at a total of 8,000-10,000 trajectories for each sea-quark mass in each ensemble. The increased statistics will reduce the statistical error. Furthermore, we have been performing simulations on two larger lattices, 20 3 × 40 × 16 and 24 3 × 48 × 16, for β = 5.95 and β = 5.90 respectively. This will allow us to study the volume dependence as well as the discretization error.
To conclude, our results of the mass and the decay constant of the pseudoscalar meson are in good agreement with the sea-quark mass dependence predicted by the next-to-leading order (NLO) ChPT, and provide a determination of the low-energy constantsl 3 andl 4 , the pion decay constant, the chiral condensate, and the average up and down quark mass. Together with our recent result of the topological susceptibility [7], these suggest that the nonperturbative chiral dynamics of the sea quarks are well under control in our HMC simulations. Moreover, this study also shows that it is feasible to perform largescale simulations of unquenched lattice QCD, which not only preserve the chiral symmetry to a high precision, but also sample all topological sectors ergodically. This provides a new strategy to tackle QCD nonperturbatively from the first principles.