Quark masses and low energy constants in the continuum from the tadpole improved clover ensembles

We present the light-flavor quark masses and low energy constants using the 2+1 flavor full-QCD ensembles with stout smeared clover fermion action and Symanzik gauge actions. Both the fermion and gauge actions are tadpole improved self-consistently. The simulations are performed on 11 ensembles at 3 lattice spacings $a\in[0.05,0.11]$ fm, 4 spatial sizes $L\in[2.5, 5.1]$ fm, 7 pion masses $m_{\pi}\in[135,350]$ MeV, and several values of the strange quark mass. The quark mass is defined through the partially conserved axial current (PCAC) relation and renormalized to $\overline{\mathrm{MS}}$ 2 GeV through the intermediate regularization independent momentum subtraction (RI/MOM) scheme. The systematic uncertainty of using the symmetric momentum subtraction (SMOM) scheme is also included. Eventually, we predict $m_u=2.45(22)(20)$ MeV, $m_d=4.74(11)(09)$ MeV, and $m_s=98.8(2.9)(4.7)$ MeV with the systematic uncertainties from lattice spacing determination, continuum extrapolation and renormalization constant included. We also obtain the chiral condensate $\Sigma^{1/3}=268.6(3.6)(0.7)$ MeV and the pion decay constant $F=86.6(7)(1.4) $ MeV in the $N_f=2$ chiral limit, and the next-to-leading order low energy constants $\ell_3=2.43(54)(05)$ and $\ell_4=4.322(75)(96)$.


I. INTRODUCTION
As fundamental parameters of the standard model which are not directly measurable in experiments, the mass of the lightest three flavors can only be determined accurately using lattice quantum chromodynamics (QCD).Lattice QCD offers a non-perturbative approach to solve QCD, the underlying theory of the strong interactions, but a set of complete and accurate ensembles is essential to ensure reliable results.
Thus a natural question is, whether it is possible for the clover fermion to reach a high accuracy determination of the light quark masses.In 2007, Ref. [20] proposed an alternative approach to define light quark masses from the PCAC relation and renormalizes it with tadpole improved 1-loop matching.By utilizing the Schrodinger functional (SF) scheme [21], the PCAC quark mass can be renormalized non-perturbatively, and the calculation with physical pion mass but a single lattice spacing a =0.09 fm and m π L ∼ 2 gives m u,d = 3.12(24)(08) MeV [22][23][24][25] at MS(2 GeV) which is 10% lower than the present lattice average value m u,d = 3.381 (40)MeV [26] with large uncertainty.
A more systematic study using the SF scheme was conducted by the ALPHA collaboration with multiple lattice spacing a ∈ [0.05, 0.086] fm but relatively heavy quark masses m π ≥ 200 MeV, and their determination resulted in m u,d = 3.54 (12) (9)MeV [27].So far, the most precise determination of m u,d = 3.469 (47) (48) [27,28] with clover fermion comes from the BMW collaboration, which utilized multiple lattice spacings a ∈ [0.05, 0.012] fm with the lightest pion mass m π = 131 (2) MeV and renormalized the quark mass using the widely-used RI/MOM scheme [29].
But the systematic uncertainty of using the RI/MOM scheme could be underestimated, as the RI/MOM scheme exhibits poor perturbative convergence for the scalar/pseudoscalar current, leading to sensitivity in the final result due to the estimate of the missing higherorder corrections.Thus, the SMOM scheme [30,31] was proposed to suppress this uncertainty and has been employed in most recent quark mass determinations using chiral fermions.Nevertheless, a recent study [32] at multiple lattice spacings shows that using either the RI/MOM or SMOM intermediate scheme can result in the renormalized scalar current under MS scheme differing by 30% at a ∼ 0.1 fm for the clover fermion.
Additionally, it is worth mentioning that the renormalized quark mass using the RI/MOM scheme with the Twisted-mass fermion [33][34][35] is m u,d = 3.64 (7) (6) MeV, which is approximately 5% higher than the results obtained with chiral fermions that predominantly use the SMOM scheme.
In this work, we conduct a detailed comparison of the renormalization constants (RCs) using the RI/MOM and SMOM schemes.It turns out that the sensitivity of the intermediate schemes can be suppressed to ∼ 5% level, which allows us to provide a relatively precise prediction of the quark mass.Based on the Kaon masses with the QED effect subtracted, we also obtain the up, down, and strange quark masses separately, along with other related quantities.We expect that further improvement in the prediction accuracy can be achieved through calculations on more lattice spacings.

II. SIMULATION SETUP
The results in this work, are based on the 2+1 flavor full QCD ensembles using the tadpole improved tree level Symanzik (TITLS) gauge action and the tadpole improved tree level Clover (TITLC) fermion action.
The TITLS gauge action, denoted as S g , is defined in the following, where N c = 3, and The TITLC fermion action uses 1-step stout smeared link V with smearing parameter ρ = 0.125, S q (m) = x,µ=1,...,4,η=± ψ(x) 1 + ηγ µ 2 V ηµ (x)ψ(x+η μa) where c sw = 1 v 3 0 with v 0 = ⟨ ReTr x,µ<ν P V µν (x) 6Nc Ṽ ⟩ 1/4 , and The parameters utilized for the simulation, encompassing the lattice size ( L3 × T ), gauge coupling ( β), and the lattice spacing (a) determined through the gradient flow [36] with w 0 [37] using the Symanzik action, are outlined in Table I.The dimensionless bare degenerated light and strange quark mass ( mb l,s ), renormalized quark masses (m R l,s ) at MS 2GeV, and the respective pion and kaon masses (m π,K ) are also included in the table.The details of the pseudoscalar meson mass and lattice spacing extraction can be found in the supplemental materials [38].The impact of the mistuning effect of the tadpole improvement factors u 0 and v 0 can also be found there.The ensemble set used in this work is designed to control the variables in the systematic uncertainty estimation.For example, the spatial size L of the C24P29, F32P30, and H48P32 ensembles are all within 1-2% of each other, and the unitary pion masses are also similar with a 10% difference.Thus, they are very suitable for investigating the discretization error of the hadron structure with non-zero given momentum.The pion mass and volume of C32P23 are close to those of F48P21 within 10%, and remaining differences can be further suppressed by interpolation with the other ensembles or by generating a new ensemble C36P21 using interpolated parameters.The other ensembles with larger dimensionless volume, such as F64P14 and/or H64P22, should also be helpful in achieving better control over the discretization error, and will be generated in the future.
For the Clover fermion action, defining the renormalized quark mass m R q from the bare quark mass parameter mb q can be subtle since the critical quark mass mcrti vanishing the pion mass is non-zero.A more practical solution defines it through the PCAC relation [20]: where A µ = ψγ 5 γ µ ψ and P = ψγ 5 ψ.The PCAC quark mass m PC q is then defined through the pion correlation functions: where m PS is the pseudoscalar meson mass.The renormalized quark mass is subsequently defined as m R q = Z A /Z P m PC q .In Fig. 1, we plot the dimensionless PCAC quark mass mPC q ≡ m PC q a as a function of the dimensionless input bare quark mass parameters mb q ≡ m b q a, at three lattice spacings with m π ∼ 300 MeV.The figure also includes linear fits using the following form: where mcrti corresponds to the critical pion mass that makes the pion mass and mPC vanish.= m PC q a v.s. the bare quark mass mb q = m b q a at three lattice spacings.The slope should approach 1 in the continuum limit.The data points correspond to six valence quark masses around the unitary light and strange quark masses in those ensembles.
k m = 1 + O(a 2 , α s , aα s ) approaches 1/Z A determined by non-perturbative RI/MOM renormalization (due to the relation Z m Z P = 1) in the continuum limit, while it is affected by the O(a 2 ) discretization error and O(α s ) loop effects at finite lattice spacing.
Unlike the hadron mass, the determination of physical quark mass on the lattice using discretized actions requires additional renormalization.The RCs defined under the MS scheme, can only be obtained through regularization-independent (RI) schemes such as RI/MOM [29] or SMOM [30,31].These RCs should be independent of intermediate schemes.
For the overlap fermion action, which possesses strict chiral symmetry, the relations Z V = Z A and Z P = Z S = 1/Z m are satisfied strictly.The scheme dependence of RI/MOM or SMOM can be ignored compared to other systematic uncertainties [39].
However, in the case of the clover fermion action, the ratio Z A /Z V can deviate from unity due to the additive chiral symmetry breaking effect generated by additional terms in the action , and it is sensitive to the choice of using RI/MOM or SMOM scheme.Based on the f π at three lattice spacing with m π = 317 MeV, the scheme sensitivity is approximately 1% after a linear a 2 continuum extrapolation, with the discretization error through RI/MOM being 25% smaller than that of SMOM.0.000 0.002 0.004 0.006 0.008 0.010 a 2 (fm 2 ) For the chiral symmetry breaking effect between Z S and Z P , it is valuable to consider the scalar matrix element with the valence quark contribution only, g S,π,ME = ⟨π|S|π⟩ val /⟨π|π⟩, (7) where S = ψψ.We show the renormalized g g S,π obtained from RI/MOM and SMOM schemes in Fig. 2, for m π = 317 MeV at three different lattice spacings.It can be observed that the RI/MOM scheme exhibits a smaller discretization error than that of the SMOM scheme, and the continuum extrapolated values differ from each other by approximately 7.6(2.3)%.
Using the Feynman-Hellman theorem, one can also extract g S,π from the quark mass dependence of m π , as where the factor 1 2 in front of ∂mπ(mq) ∂mq is used to average the contribution from two propagators in the pion correlator.Using the renormalized quark mass m (filled green dots for MOM and red boxes for SMOM) are in the range of 3.9 to 4.2 and then slightly smaller than the linear a 2 extrapolated value g S,π,ME =4.35 (9) (green band) using the RI/MOM scheme but consistent with the SMOM value 4.02(6) (red band) .Even more, g S,π,FH is consistent with the g S,π,ME using the RI/MOM scheme at each lattice spacing within two sigma, but have significant difference from the g S,π,ME using the SMOM scheme.Thus the deviation between the g S,π,FH and g S,π,ME using the RI/MOM scheme would be only a systematic uncertainty due to the linear a 2 extrapolation.Thus the renormalized g S shall have about 7% systematic uncertainty with present data, and more reliable continuum extrapolation with data at more lattice spacing is essential to obtain accurate prediction on g S .
The renormalization constants for various quark bilinear operators are detailed in the supplemental materials [38], along with a discussion on the discretization error from different renormalization methods of quark field and mass.

III. RESULTS
Using the lattice spacing shown in Table I, we find that the unitary pion mass on the ensemble C48P14 at a = 0.1053(2) fm is 135.5(1.6)MeV, which perfectly agrees with the physical neutral pion mass m π 0 of 134.98 MeV within 1% statistical uncertainty.The charged pion mass m π ± 139.57MeV receives the QED correction 4.53 (6) MeV [40] and then the subtracted pure QCD m π ± is consistent with that of m π 0 within the uncertainty.
The corresponding renormalized light quark mass and pion decay constants can also be determined as: m MS(2 GeV) l (a = 0.105 fm) = 3.64(8) (11) MeV, f π (a = 0.105 fm) = 121.9(5)MeV, (9) where the second uncertainty of m q comes from that of the renormalization constant.Based on the continuum extrapolation with a 317 MeV pion mass, the pion decay constant can change by approximately 7%, and then agree with the present PDG value 130.4(0.2)MeV[41] after the continuum extrapolation.
In order to process this continuum extrapolation systematically, we calculate the quark propagators with unitary light quark mass and also 2 partially quenched quark masses with the constraint m π L > 3.5, on each of the 11 ensembles.Then we use the NLO partially quenched χPT form [42] to describe the pion masses and decay constants with different valence and sea quark masses, in addition to extra parameters c m/f,a/l for the finite lattice spacing/volume corrections.
Since the statistics on each ensemble are different, we perform 4000 bootstrap re-samplings on each ensemble and conduct the correlated global fit based on these bootstrap samples.In such a strategy, the correlation between different data points in the same ensemble is included automatically, and that between different ensembles vanishes within the statistical uncertainty of the re-sampling.The lattice spacing and renormalization constants are sampled for each bootstrap sample using a Gaussian distribution with their uncertainties as the width of the distribution.[26] and/or PDG [41].The difference between two schemes is considered as systematic uncertainty in the combined determination.varies with the quark mass, at three lattice spacing (colored data points and corresponding bands with dashed line for 0.0519 fm, dash-dotted line for 0.0775 fm, and dotted line for 0.1053 fm) and also continuum (gray band).
To illustrate the lattice spacing dependence and the unitary quark mass dependence, we subtract the partially quenching effect using bootstrap samples of the fit parameters from the original data points, and show the ratio (m π ) 2 /m MS(2 GeV) q (upper panel) and also f π (lower panel) at different quark masses m q .The corrected data points at different lattice spacings use different symbols: red crosses for a=0.105 fm, blue triangles for a=0.077 fm, and filled green dots for a=0.052 fm.The bands with similar color represent the fitted band at the corresponding lattice spacing, and the gray band shows the final prediction in the continuum and infinite volume limit.It is observed that the continuum extrapolation pushes f π to be obviously higher, while the impact on the (m 2 π /m q ) ratio and, consequently, m q is much weaker.
The m π and f π with unitary valence and sea quark masses have the following parameterization, where . Σ, F and ℓ 3,4 are low energy constants.Our determination of those constants are also collected in Table II, consistent with the current N f = 2 + 1 FLAG average but have smaller uncertainties except F .
In this work, we use the m K ± and m K 0 with the constraint m phys u + m phys d = 2m phys l , to determine the up, down and strange quark masses m u,d,s .The QED correction on the kaon mass is subtracted based on the literature [34] under the renormalization scheme m MS q,QCD+QED (2GeV) = m MS q,QCD (2GeV).On each ensemble, we calculate the strange quark propagators with a unitary strange quark mass m v s = m s s , and also 2 partially quenched quark masses m v s ∼ 100 MeV.We construct the Kaon correlation functions with three strange quark masses and three light quark masses used in the pion case.The 3 × 3 partially quenched Kaon masses on all the ensembles are fitted with the following form proposed in a recent work [35], Based on the fit of m 2 K , the total strange quark mass dependence b v s + b s s = 2.37(08) is consistent with the leading order light quark mass dependence b v l + b s l = 2.59(95), and the coefficient of the non-linear quark mass dependence c K l = 0.07 (13) can not be determined based on current statistics.In Fig. 4, we show the corrected kaon mass m cr K and decay constant f cr K with the light quark mass m l corrected to its physical value m phys l .The finite volume and partially quenched effects are also subtracted.We can found that f K also exhibits a strong lattice spacing dependence, similar to the f π case, while the Kaon mass is insensitive to the lattice spacing.
As illustrated in Fig. 3 and Fig. 4, all the global fits of the pseudoscalar meson mass and decay constant provide reasonable χ 2 /d.o.f.More information on the global fit can be found in the supplemental materials [38].
The physical quark masses m u,d,s and also corresponding f π,K using m phys l and intermediate RI/MOM scheme, are collected in Table II.In addition, Table II shows the global fit results using the Z A,P through the SMOM scheme for comparison.As we can see from the continuum extrapolation tests using a 317 MeV pion mass, the SMOM scheme yields quark masses that are 3-4% lower and decay constants that are ∼2% lower compared to the RI/MOM scheme.However, the ratio of the quark masses or decay constants remains unchanged within the uncertainty as the renormalization constants are cancelled.
Therefore, we consider the result using the RI/MOM scheme as the central value due to its smaller discretization error, and treat the difference between the results obtained using the two schemes as systematic uncertainties.Such a systematic uncertainty can also be considered as an estimate of the residual discretiation error, as the correct continuum limit should be independent of the intermediate renormalization scheme.All our determinations are consistent with the present lattice averages [26] and/or PDG [41] within 1-2 sigma.

IV. SUMMARY
In this work, we determine the up, down, and strange quark masses, along with several low energy constants, using the 2+1 flavor full-QCD ensembles with tadpoleimproved clover and Symanzik actions.The major results are summarized in Table II.
Similar to one of the most precise works [28] with the clover fermion, we have skipped the axial current improvement [43] since the improvement coefficient itself strongly depends on the lattice spacing a and bring the improvement term to be consistent with an O(a 2 ) correction.As evidence, both f π and m q show good consistency with a simple linear a 2 lattice spacing dependence.Then assigning c A ∼ 0.05a/(0.105fm) in the improved axial vector current A imp µ = A µ + c A a∂ µ P [43] can eliminate the discretization error in f π,K , however, this error will be transferred to the quark mass.Thus simulations at more lattice spacings would be a more systematic solution to enhance the accuracy of our predictions in the continuum than the axial current improvement.
On the other hand, the additive chiral symmetry breaking makes the renormalization of the quark mass to be highly non-trivial.Our study suggests that the Z S and Z P obtained through the SMOM scheme are closer than those through the RI/MOM scheme, while the latter one can make the discretization error of both the m R q = Z A /Z P m PC q and g R S = Z S g S to be smaller.Our final prediction of the quark masses are 5.6(2.8)%higher than the current 2+1 flavor lattice averages but consistent with the previous 2+1 and 2+1+1 flavor results using the RI/MOM scheme.At the same time, the RI/MOM scheme can also cause the Feynman-Hellman theorem g S,π ≃ mπ 4mq to be violated by 7(3)% after the linear O(a 2 ) continuum extrapolation.
Using the SMOM scheme can eliminate the violation and bring the quark mass prediction closer to the current 2+1 flavor lattice average.However, the SMOM scheme introduces larger discretization errors for all the renormalized quantities we investigated and causes the decay constants f π,K to be 2-3% smaller than the physical values after the linear O(a 2 ) continuum extrapolation.
The above observations indicate that renormalization is a significant issue that requires careful investigation, and conducting similar calculations using chiral fermions would be essential to gain a better understanding of these violations.At the same time, non-perturbative renormalization should remove all the O(α s ) effects, but not all the cross terms like the residual O(aα s ) effect of the clover action, which can cause the O(a 2 ) continuum extrapolation to fail.Thus, we consider the difference between the results obtained by the two schemes as systematic uncertainties in our final determination of the aforementioned quantities, which are larger than the statistical uncertainties in various cases.We anticipate that additional research utilizing ensembles with a greater number of lattice spacings can encompass both the O(a 2 ) and O(aα s ) terms in the continuum extrapolation, resulting in a more dependable and uniform continuum limit.
It is worth mentioning that in Ref. [28], the tracesubtraction trick S = S− 1  4 Tr[S] is applied into renormalization procedure, and the quark mass is renormalized at RI/MOM 2 GeV, followed by perturbative matching at a much higher scale.This approach is crucial in suppressing their truncation error to the sub-percent level.However, in our case, it appears to be inefficient due to significant non-perturbative effects observed at 2 GeV.We plan to conduct a more systematic investigation once the CLQCD ensembles at more lattice spacings are gen-erated.

APPENDIX A. Simulation details
This section will be organized in the following: The dimensionless joint fit on the pseudo-scalar meson mass, its decay constant, and the corresponding partially conserved axial currennt (PCAC) quark mass will be discussed in Sec.IV A 1; Based on the determination of an uniform lattice spacing at given β detailed in Sec.IV A 2, the mistuning effect of the tadpole improvement factors is not always negligible and will be addressed in Sec.IV A 3.

Dimensionless joint fit
We construct two kinds of the two-point functions for the meson states: where C 2 is independent of the source time slice y 0 after taking expectation value, the Coulomb gauge fixed wall source propagators is defined as S(x; y; m q , U ) ≡ ψ(x; m q , U ) ψ(y; m q , U ) is the quark propagator of the quark field ψ with bare quark mass m q on a given gauge configuration U , and U C is the Coulomb gauge fixed configuration satisfying the gauge fixing con- For the Clover fermion action, the PCAC quark mass m q is then defined through the pion correlation functions [20]: The renormalized quark mass is subsequently defined as m R q = Z A /Z P m PC q .Through a joint fit ( Õ is the dimensionless value of any quantity O), CPP 2,wp ( t) CPP 2,ww ( t)  the PCAC quark mass mPC q , pseudo-scalar mass and decay constant fPS are extracted as fit parameters, alongside an additional unphysical fit parameter Z wp for the Coulomb gauge-fixed wall source.In Fig. 5, we shown the joint fit result of the unitary light quark on the physical point ensemble C48P14.To suppress the statistical uncertainty, we repeated the calculation on n src =48 of 96 times slides on n cfg =203 configurations.The values of n src and n cfg of the other ensembles can be found in Table.I.The ratios (red data points) defined in Eq. ( 17) for mq are consistent with a constant in the region of t > 10, as shown in the top left panel.The situation is similar for fPS (top right panel) and also mPS (bottom left panel).Since the statistics on C48P14 are limited, we performed a joint uncorrelated fit in the range t ∈ [10,40] for Eq.(17)(18)(19) and used bootstrap resampling to estimate the uncertainty, represented by the gray band in the figure.We can see that the fit agrees very well with the data points, and the uncertainties of the fitted bands are comparable to the original data.Thus, it is unlikely that the uncertainty has been underestimated by the uncorrelated fit.
The strange quark mass case on the ensemble C48P14 is shown in Fig. 6.Since the statistical uncertainty decreases with with a heavier quark mass, we only repeated the calculation on 3 time slides for each configuration, resulting in more noticeable fluctuations compared to the light quark case.Through the above fit, we can extract the light and strange quark masses through the pion and η s correlators, respectively: Alternatively, the sum of the light and strange quark masses can also be extracted from the kaon correlator, Thus we should verify two extraction provide consistent result, at least in the continuum.Fig. 7 shows the ratio of these two determinations of m l +m s at three lattice spacings, with m π ∼ 300 MeV and m ηs ∼ 700 MeV for better signals.Two determinations deviate from each other by approximately 0.6% at the coarsest lattice spacing but are consistent within the statistical uncertainty of 0.1% after a linear a 2 extrapolation to the continuum limit.
At three lattice spacings with m π ∼ 300 MeV, we fit the dimensionless PCAC quark mass mPC q = m PC q a with the following form: where mb q = m b q a is the original input quark mass parameters, and mcrti corresponds to the critical pion mass that makes the pion mass and mPC vanish.The parameter k m = 1 + O(a 2 , α s , aα s ) approaches 1/Z A determined by non-perturbative RI/MOM renormalization (due to the relation Z m Z P = 1) in the continuum limit, while it is affected by the O(a 2 ) discretization error and O(α s ) loop effects at finite lattice spacing.
Based on the numerical results listed in Table III, we observe that mcrti remains negative even , in the continuum extrapolation to ensure that mcrti −−→ a→0 0, as predicted by lattice perturbative theory.
On the other hand, the discrepancy between k m and 1 diminishes as the lattice spacing becomes smaller.We will discuss this further in the following section on renormalization.Fig. 8 illustrates the lattice spacing determined by the gradient flow [36] with w 0 [37] at three bare couplings β and m π ∼ 300 MeV, using different gauge action improvement coefficient c 1 in the flow: c 1 = −0.331(Iwasaki action), −0.2 and −1/(12u 2 0 ) for interpolation, −1/12 (Symanzik action), and 0 (Wilson action).The c 1 dependence becomes weaker at smaller lattice spacings, indicating that it is a discretization effect.As the tadpole improvement factor approaches 1 with increasing gradient flow t, implementing tadpole improvement for the action used by the flow only affects the small flow time region and is unnecessary.Thus, we consistently use the standard c 1 = −1/12 in the gradient flow to match the gauge action employed in HMC and mitigate the discretization error.Using the FLAG average value of the gradient flow scale, w 0 = 0.1736(9) fm [26], the lattice spacing a w0 for each ensemble is determined using the Symanzik flow and summarized in Table IV.It is evident that a w0 primarily depends on the gauge coupling β, while the quark mass also exhibits significant effects, as indicated by the precise a w0 results with a statistical uncertainty of approximately 0.1%.Empirically, a w0 can be described by the following parameterization, yielding a χ 2 /d.o.f. of 2.2 (with a corresponding p-value of 0.09, still larger than the standard lower bound of 0.05): Here, m π,phys = 134.98MeV [41] represents the physical pion mass without QED correction, and m ηs,phys = 689.63(18) MeV [49] corresponds to the pseudoscalar meson mass of the strange quark with only connected insertions.The tadpole improvement factors used in the actions are denoted as u I 0 and v I 0 , while u 0 and v 0 are tadpole improvement factors obtained from the generated configurations.The lattice spacing a(β) with physical quark masses at each β, along with the other fitting parameters c l,s,L,u0,v0 are also provided in Table IV.
We observe that the volume dependence is consistent with zero, but there is still a non-vanishing mismatch effect (u I 0 ̸ = u 0 and v I 0 ̸ = v 0 ) in the tadpole improvement factors.It is also interesting to note that the dependence of a w0 on the strange quark mass is weaker compared to that of the light quark mass by a factor of approximately 4, instead of 2 from two light flavors.

Mismatch effect of the tadpole improvement factors
As shown in Table IV, the mismatch effect of the tadpole improvement factors can have a non-zero impact on the determination of a w0 through the gradient flow.Both u 0 and v 0 represent vacuum expectation values and cannot be accurately determined before generating gauge configurations using the realistic Hybrid Monte Carlo (HMC) production.Therefore, for each ensemble, we initiate the HMC with an initial guess for u 0 and v 0 , measure their values on each trajectory until they stabilize, and subsequently restart the HMC production with the updated values of u 0 and v 0 .After several iterations, the input values u I 0 and v I 0 become consistent with u 0 and v 0 , as measured from the configurations, at a level of approximately 0.002%, resulting in an effect of around 0.6% on the lattice spacing.The only exception is the test ensemble C24P34, which exhibits a larger deviation in both u 0 and v 0 ; however, their impact on aw 0 is mainly canceled out due to the opposite signs of c u0 and c v0 .After averaging u 0 and u I 0 to obtain an estimate of the self-consistent tadpole improvement factor ū0 , we observe that the volume dependence is at the level of 0.001% based on the spatial sizes utilized.On the other hand, the quark mass dependence appears to be more pronounced.
In the upper panel of Fig. 9, the value of ū0 exhibits a linear behavior with respect to the combined quark mass 2 mPC l + 0.55(1) mPC s ∝ 0.5 m2 π + 0.14 m2 ηs .This suggests that the dependence on the strange quark mass exhibits a similar suppression as observed in a w0 , as shown in Table IV.In the lower panel, we observe that v 0 exhibits  the same quark mass dependence within the uncertainty.However, such a mismatch will not significantly alter the physical observables.To illustrate this, let's consider the C24P34 ensemble with the largest mismatch effect.We calculate the pion two-point function C P P 2,wp (t; mb l , c sw ) using the same mb l = −0.2770,but with either c sw = 1/(v I 0 ) 3 or 1/v 3 0 .The joint fit parameters defined in Eq. (17)(18)(19) are listed in TableVI.We find that the matrix elements fPS and Zwp remain unchanged within 0.05% uncertainties.However, mPS experiences a 0.6% change, resulting in a 1.1% shift in the quark mass.This can be understood by noting that increasing c sw leads to an increase in mcrti , thereby making the multiplicative renormalizable quark mass mPC l larger for the same mb l .Consequently, the mismatch effect falls within the statistical uncertainty of C24P34 and is an order of magnitude smaller in other ensembles.Hence, we can safely disregard this mismatch effect in the subsequent discussions.

B. Renormalization
Unlike the hadron spectrum, the determination of hadron matrix elements on the lattice using discretized actions requires additional renormalization.The RCs defined under the MS scheme, can only be obtained through regularization-independent (RI) schemes such as RI/MOM [29] or SMOM [30,31].These RCs should be independent of intermediate schemes.
For the overlap fermion action, which possesses explicit chiral symmetry, the relations Z V = Z A and Z P = Z S = 1/Z m are guaranteed.The consistency of using RI/MOM or SMOM schemes has been verified within systematic uncertainties [39].However, in the case of the clover fermion action, which exhibits additive chiral symmetry breaking, additional considerations and discussions regarding its impact on renormalization are necessary in this section.Following a similar strategy employed by Ref. [39], the values of Z S,P,T incorporate two sources of uncertainty.The first one encompasses ensemble-independent statistical and systematic uncertainties, including lattice spacing, finite volume effects, and a 2 µ 2 fit range.The second source of uncertainty arises from perturbative matching, including uncertainties associated with Λ QCD , truncation in the perturbative matching, and perturbative scale running.These uncertainties are fully correlated across different ensembles.As shown in Table VII for the C32P29 ensemble, the truncation error is the largest source of uncertainty for Z S,P .
In Ref. [39], the truncation error resulting from the 3-loop perturbative matching between the RI/MOM and MS schemes is estimated by introducing a fake 4-loop correction.As an illustration, considering the scalar/pseudoscalar case with the largest truncation error, the 3-loop matching and the fake 4-loop one are expressed as follows: The coefficient of fake α 4 s term 2.722 2 /1.007 = 7.358 is the same as that provided by the Padé approximation, Based on the precise 4-loop calculation, the precise 4loop matching is [50], which is larger than C MS,RI/MOM S,n f =3,fake4−loop but smaller than C MS,RI/MOM S,n f =3,Pade3 .Considering that the existence of the 5loop correction, we anticipate that the Padé approximation will provide a more accurate estimate for C S,n f =3 .
Thus we shall use the 4-loop Padé approximation, 29) as the central value of the matching coefficent, and taking the difference between C S,n f =3,Pade4 and C S,n f =3,Pade3 as systematic uncertainty of the truncation on the perturbative matching.
In this section, we will commence with the vector current normalization and explore different choices of quark field renormalization in Section IV B 1. Section IV B 1 will then present the investigation of the chiral symmetry breaking effect between Z V and Z A .Subsequently, in Section IV B 4, we will examine the analogous investigation of Z P and Z S , with additional discussions on Z m in Section IV B 3. Finally, for completeness, the case of the tensor current will be discussed in Section IV B 5. Unlike its continuum counterpart, the local vector current under lattice regularization is subject to both the O(a 2 ) discretization error and O(α s ) loop corrections [51], and then requires additional normalization.The normalization constant Z V can be determined from the vector current conservation condition, where V ν = ψγ ν ψ and H represents an arbitrary hadronic state.Thus extracting Z V from the pion correlator in the rest frame will be the cheapest choice [51], ⟨Tr[S † C (⃗ x, T /2; 0; m q )S C (⃗ x, T /2; 0; where the additional propagator in the denominator can be obtained using the sequential source technique.The ratio at three lattice spacing and their constant fits in the range of 0 ≪ t ≪ T /2, are shown in Fig. 10.As demonstrated in [51], Z V is affected by α s corrections and cannot be accurately extrapolated to 1 using a simple naive a 2n continuum extrapolation.The quark field RC Z RI ′ q in the regularization independent (RI) scheme can be accessed through either its definition with S(p) = x e −ip•x ⟨ψ(x) ψ(0)⟩, or the vertex correction of the vector current which is equivalent in the continuum, where These two definitions are equivalent in the continuum, but they are subject to different discretization errors.Fig. 11 illustrates the comparison between Z RI ′ q and Z RI ′ ,ver q at three lattice spacings.The two definitions become closer either at smaller µ or smaller a, with Z RI ′ ,ver q exhibiting a significantly smaller a 2 p 2 error.This situation is reminiscent of the overlap fermions case [39], where the a 2 µ 2 errors have opposite signs.
Given that Z RI ′ ,ver q involves an off-diagonal projection that introduces additional discretization errors, it is more convenient to define the quark field renormalization through the standard RI definition [29], The discrepancy between Z ω q with different ω and Z RI ′ q should be eliminated through their respective perturbative matchings, which can be calculated in the continuum.For the SMOM condition with ω = 1, the convergence of the perturbative matching is poorer compared to that of Z q using the q-projection definition [31], but the discretization error will be smaller.12.The ratio ZA/ZV through the RI/MOM scheme (data points) or SMOM scheme (bands) at three lattice spacing as functions of a 2 p 2 .
For the clover fermion, the ratio can deviate from unity due to the additive chiral symmetry breaking present in the action.As depicted in Fig. 12 for the RI/MOM scheme (data points, ω = 0) and SMOM scheme (bands, ω = 1) at three different lattice spacings with m π ∼ 300 MeV, we observe that the breaking diminishes as a 2 µ 2 increases, indicating that the chiral symmetry breaking in the mass term becomes less important.However, the ratio does not approach unity in the continuum limit without any α s corrections.It is worth noting that the breaking using the SMOM scheme is much smaller than that with the RI/MOM scheme, while increases rapidly at small a 2 µ 2 .
0.000 0.002 0.004 0.006 0.008 0.010 0.012 a 2 (fm If we fit the Z A /Z V data with a polynomial form in the range of 9GeV 2 ≤ p 2 ≤ 15/a 2 and extrapolate it to a 2 p 2 = 0, the obtained result will be sensitive to the choice of ω.Therefore, we extrapolate the f π at two coarser lattice spacings to the unitary pion mass at the finest lattice spacing, which is m π = 317 MeV.The renormalized f R π = Z A f π is shown in Fig. 13, using both the extrapolated Z A /Z V obtained through RI/MOM and SMOM.The results suggest that the scheme sensitivity is approximately 1% after a linear a 2 continuum extrapolation, with the discretization error through RI/MOM being 25% smaller.After Z ω P is extracted, it can be further converted to the MS scheme at given scale likes µ 0 = 2 GeV through the perturbative matching, and shall be independent of both µ and ω.In Fig. 16, Z MS(2 GeV) P obtained through the RI/MOM scheme at the coarsest lattice spacing a=0.105 fm can be consistent with that using the SMOM scheme, with the 3-4% systematic uncertainties from the perturabative matching.The situation improves at smaller lattice spacing, and systematic uncertainties are also smaller there.In Fig. 17, we show the renormalized quark mass at MS(2 GeV) with m π =317 MeV at three lattice spacing, using Z MS(2 GeV) P through either RI/MOM or SMOM scheme.We can see that the lattice spacing dependence using the RI/MOM scheme is consistent with zero, while SMOM shows a non-vanishing dependence and make the continuum extrapolated value to be 3.1(1.5)%lower.The scalar current RC can be defined similarly with a slightly different parameterization, where S = ψψ, and the A ω S term can be dropped for ω = 1.However, as shown in Fig. 18 through the RI/MOM (ω = 0, upper panel) and SMOM (ω = 1, lower panel) schemes, at three lattice spacing and mπ,ss ∼ 300 MeV.
After converting to the MS scheme at a given scale using the same matching procedure as in the pseudoscalar case, we obtain Z MS(2 GeV) S as shown in Fig. 19, through either the RI/MOM scheme (upper panel) or the SMOM scheme (lower panel) with the corresponding a 2 µ 2 errors.We observe that the scheme dependence of Z S is much stronger than that of Z P , but the difference diminishes as the lattice spacing decreases.
Another comparison we can make is the ratio Z S /Z P obtained through either the RI/MOM or SMOM scheme.As depicted in Fig. 20, the chiral symmetry breaking effect is significantly smaller when using the SMOM scheme (bands) compared to the RI/MOM scheme (data points).Moreover, both schemes exhibit further suppression at smaller lattice spacings.Similar to the Z A /Z V case, the inclusion of the O(α s ) term is necessary to restore chiral symmetry in the continuum.For completeness, we also show the tensor current renormalization using both the RI/MOM and SMOM schemes, where T µν = ψσ µν ψ. through the RI/MOM (ω = 0, upper panel) and SMOM (ω = 1, lower panel) schemes, at three lattice spacing and mπ,ss ∼ 300 MeV.
performing a linear m 2 π extrapolation to the continuum, as shown in Fig. 21.The systematic uncertainty arising from the SMOM matching is estimated to be 2-3%, which implies that Z MS(2 GeV) T is considered independent of the intermediate scheme within this uncertainty.ever, it is important to note that the scheme dependence is not completely eliminated even after performing the continuum extrapolation.Considering the lattice spacing dependence of the renormalized f π , m q , and g S , it is recommended to use the RI/MOM scheme to suppress the discretization error.
In Table VIII, we also show the renormalization constants using the RI/MOM scheme but with an additional 1/p 2 term [32].This term can have an obvious impact at the coarsest lattice spacing but is consistent with zero at the finest lattice spacing.This is understandable as the fitting range at finer lattice spacing is much larger, reducing the possible influence of a 1/p 2 pole in the inferred region.Thus, it can be considered a discretization effect that will be eliminated in the continuum extrapolation.
As a summary of this section, Table IX summarizes the normalization and renormalization constants for all the ensembles used in this work, obtained through the intermediate RI/MOM scheme.The first uncertainty of RCs is ensemble-independent, while the second one is fully correlated on different ensembles and can be suppressed after the continuum extrapolation.

C. Global fit
In order to process this continuum extrapolation systematically, we calculate the quark propagators with unitary light quark mass and also 2 partially quenched quark masses with the constraint m π L > 3.5, on each of the 11 ensembles.Then we use the following NLO partially quenched χPT form [42] to describe the pion masses and decay constants with different valence and sea quark masses, in addition to extra parameters c m/f,a/l for the finite lattice spacing/volume corrections, where N f = 2 for two light flavors and Λ χ = 4πF is the intrinsic scale of χPT with F being the pion decay constant in the chiral limit.The dimensionless expan- involve the chiral condensate Σ and quark mass m v/s for the valence and sea quark masses, respectively.Additionally, α i represents the NLO low energy constants of χPT at the intrinsic scale Λ χ .These constants can be converted to the usual scale-dependent partially quenched χPT NLO coefficients as [53,54].We also introduce additional corrections terms with the coefficients c m/f,a/L/s to account for the discretization, finite volume and strange quark mismatch effects, with m ηs = 689.63(18) MeV from Ref. [49].
Since the statistics on each ensemble are different, we perform 4000 bootstrap re-samplings on each ensemble and conduct the correlated global fit based on these bootstrap samples.In such a strategy, the correlation between different ensembles vanishes within the statistical uncertainty of the resampling and is suppressed by the number of bootstrap samples.The lattice spacing and renormalization constants are sampled for each bootstrap sample using a Gaussian distribution with their uncertainties as the width of the distribution.These values are then applied to the dimensionless quantities extracted from the joint fit defined in Eq. (17)(18)(19)(20).The lattice spacing uncertainty is sampled with an uni-form seed at a given lattice spacing, since this uncertainty is fully corrected on all the ensembles at this lattice spacing.Similarly, The perturbative matching uncertainty of Z P is sampled with a uniform seed on all the ensembles, but with respective rescale factors on each ensemble.On the other hand, the uncertainty of lattice spacing and the non-perturbative uncertainty of Z P are sampled independently on each ensemble.In order to show the impact of two Z P uncertainties in the global fit, we preform the following four cases of fit, and show the results of Σ, F , α 4,5,6,8 and also c m/f,a/L/s in Table X: (1) Fitting with the statistical uncertainty δ Õ from the dimensionless observable Õ with O = m π , m PC q and f π , and also lattice spacing uncertainty δa; (2) Fitting with δ Õ, δa and non-perturbative uncertainty δ np Z P ; (3) Fitting with δ Õ, δa and perturbative uncertainty δ p Z P ; (4) Fitting with δ Õ, δa, δ np Z P and also δ p Z P .All four cases provide reasonable χ 2 /d.o.f.values, and similar values of F which are irrelevant to Z P .However, the uncertainty of Σ is highly sensitive to δZ P , as expected.By imposing the conditions y s = y v , M π,vv = M π,phys = 134.98MeV, a → 0, and L → ∞, we can extract the light quark mass m MS(2 GeV) l in the continuum and infinite volume limits from the global fit.The m l obtained from different cases of fits are also listed in Table X.
As shown in the table, the extrapolated quark mass m l = 3.60(3) obtained in Case (1) is consistent with the value 3.64(8) (11) MeV at a = 0.105 fm, as we argued before.Additionally, it has a smaller uncertainty due to the constraints from the other ensembles and the exclusion of δZ P .However, the uncertainty is enlarged to 2.7% in Case (4) where both the non-perturbative and perturbative uncertainties of Z P are included, while the central value remains almost unchanged.Treating δ p Z P as correlated across all the ensembles significantly suppresses its impact from 3% (at a = 0.105 fm) to 0.5% after the continuum extrapolation, as shown in Case (3); However, the Case (2) suggests that the δ np Z P is enlarged from 1% (at a = 0.105 fm) to 2.5% simultaneously since it is independent for different ensembles.
It is worth mentioning that the relative uncertainty of Σ and m l in Case (2-4) is almost the same and will be cancelled when we consider the renormalization independent combination Σm l .
Currently, 90% of the uncertainty of m l = 3.60(11) comes from the independent statistical uncertainty of Z P using the RI/MOM scheme on different ensembles, which can be suppressed if we adopt a more aggressive treatment on the renormalization constants.However, the current 3% uncertainty of m l is similar to the difference between the continuum extrapolation of m l at m π = 317 MeV using either the RI/MOM or SMOM scheme.Therefore, we reserve this opinion for future study, where more lattice spacings can be utilized to gain a better understanding of the chiral symmetry breaking effect in the renormalization constants.
With the m l extracted above, we can also determine the physical f π in the continuum and infinite volume limits to be 130.7(9) MeV, which is consistent with the experimental value of 130.4(2) MeV [41].Based on F extracted above, we predict F π /F = 1/2f π /F = 1.0675 (19).
As shown in Table X and inspired by the previous sections, the discretization error in quark mass is consistent with zero, while this error in f π is much more significant.On the other hand, the finite volume effect appears to be compariable with the discretization effect, with the correction on the F32P21 ensemble with the smallest m π L being about 4.6(1.1)%.
It is worth mentioning that the dependence of the pion mass on the strange quark mass is almost consistent with zero, as expected based on the direct calculation of the strange content in the pseudo-scalar meson and the Feynman-Hellman theorem [55].However, the dependence of the pion decay constant on the strange quark mass is more significant, which is partially attributed to the smaller uncertainty.Currently, the strange quark mass on the C24P34, C48P14, and H48P32 ensembles is higher than the physical value, while on all the other ensembles it is lower.Therefore, this dependence may arise from other systematic effects, and we cannot exclude this possibility with the present ensembles.Further studies with tuned quark masses would be helpful in verifying this.
To illustrate the lattice spacing dependence and the unitary quark mass dependence, we subtracted the partially quenching effect using bootstrap samples of the fit parameters from the original data points m data π,vv and f data π,vv .We define the corrected m uni π and f uni The m π and f π with unitary valence and sea quark masses y = y v = y s have another widely used parameter-

U⟩ 1 / 4
µ (x) = P exp ig 0 is the tadpole improvement factor, Ṽ = L3 × T is the dimensionless 4-D volume of the lattice, and we use Õ for the dimensionless value of any quantity O.

FIG. 4 .
FIG. 4. The corrected kaon mass m crK and the decay constant f cr K with the physical light quark mass m phys l , varies with the strange quark mass at three lattice spacing (colored data points and bands) and also continuum (gray band).

FIG. 7 .
FIG. 7. The ratio of the PCAC quark mass determined from the quarkoinum (pion and ηs), and also kaon.Two definitions are consistent upto O(a 2 ) correction.
after a naive O(a) + O(a 2 ) extrapolation to the continuum, with a value of -0.0865.Therefore, it is crucial to include the O (α b s ) 2 term, where α b s =

FIG. 8 .
FIG. 8.The lattice spacings determined by w0 with different gauge action improvement factor c1.

FIG. 10 .
FIG. 10.ZV from the vector current conservation of pion correlator at three lattice spacing.The plateau values are summarized in Table IX as the ZV of the ensembles C24P29, F32P30 and H48P32.

FIG. 11 .
FIG. 11.Comparison of Z RI ′ q (data points) and Z RI ′ ,ver q (bands) at three lattice spacing with different scale µ.The discrepancy observed in the figure is primarily attributed to the discretization error.

2 .
FIG.12.The ratio ZA/ZV through the RI/MOM scheme (data points) or SMOM scheme (bands) at three lattice spacing as functions of a 2 p 2 .

4 .
Chiral symmetry breaking between ZP and ZS

1 FIG. 20 .
FIG.20.The ratio ZS/ZP through the RI/MOM scheme at three lattice spacing as functions of a 2 p 2 .

TABLE II .
Summary of our determination on quark masses at MS(2GeV) and the other quantities , through the intermediate RI/MOM or SMOM schemes, with comparison with FLAG

TABLE III .
The bare coupling α b s , critical quark mass mcrti and slope km at three lattice spacing a and mπ ∼ 300 MeV.

TABLE IV .
aw 0 on different ensembles and the fitted values through the functional form in Eq. (24).

TABLE V .
The tadpole improvement factors u I 0 and v I 0 used in the gauge and fermion actions, u0 and v0 measured from the realistic configurations generated using u I 0 and v I 0 , and their averages ū0 and v0.

TABLE VI .
The joint fit results with the same mb

TABLE VII .
Summary of uncertainties of RCs in percentage on the C32P29 ensemble through the intermediate RI/MOM scheme.

TABLE VIII .
Normalization and renormalization constant at MS(2 GeV) at three lattice spacing with mπ,ss ∼ 300 MeV, using RI/MOM (ω = 0) or SMOM (ω = 1) as intermediate scheme.)In the RI/MOM case, the first line represents the result fit without 1/p 2 term, while the second line represents the fit with that term.Based on the values of Z A,S,P,T obtained in this work and collected in Table VIII, it can be observed that the dependence on the intermediate scheme (RI/MOM or SMOM) is reduced as the lattice spacing decreases.How-

TABLE IX .
Normalization and renormalization constants at MS(2 GeV) on all the ensembles.The RCs include two uncertainties where the former one is the ensemble independent statistical and systematic uncertainties, and the latter one is the systematic uncertainties from the perturbative matching which are fully correlated on different ensembles.