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

Zhi-Cheng Hu  Institute of Modern Physics, Chinese Academy of Sciences,
Lanzhou, 730000, China  University of Chinese Academy of Sciences, School of
Physical Sciences, Beijing 100049, China  Bo-Lun Hu  CAS Key Laboratory of
Theoretical Physics, Institute of Theoretical Physics, Chinese Academy of
Sciences, Beijing 100190, China  Ji-Hao Wang  CAS Key Laboratory of
Theoretical Physics, Institute of Theoretical Physics, Chinese Academy of
Sciences, Beijing 100190, China  University of Chinese Academy of Sciences,
School of Physical Sciences, Beijing 100049, China  Ming Gong  Institute of
High Energy Physics, Chinese Academy of Sciences, Beijing 100049, China
Liuming Liu<EMAIL_ADDRESS>Institute of Modern Physics, Chinese Academy
of Sciences, Lanzhou, 730000, China  University of Chinese Academy of
Sciences, Beijing 100049, China.  Peng Sun<EMAIL_ADDRESS>Institute of
Modern Physics, Chinese Academy of Sciences, Lanzhou, 730000, China  Wei Sun
Institute of High Energy Physics, Chinese Academy of Sciences, Beijing 100049,
China  Wei Wang  INPAC, Shanghai Key Laboratory for Particle Physics and
Cosmology, Key Laboratory for Particle Astrophysics and Cosmology (MOE),
School of Physics and Astronomy, Shanghai Jiao Tong University, Shanghai
200240, China  Southern Center for Nuclear-Science Theory (SCNT), Institute of
Modern Physics, Chinese Academy of Sciences, Huizhou 516000, Guangdong
Province, China  Yi-Bo Yang<EMAIL_ADDRESS>University of Chinese Academy of
Sciences, School of Physical Sciences, Beijing 100049, China  CAS Key
Laboratory of Theoretical Physics, Institute of Theoretical Physics, Chinese
Academy of Sciences, Beijing 100190, China  School of Fundamental Physics and
Mathematical Sciences, Hangzhou Institute for Advanced Study, UCAS, Hangzhou
310024, China  International Centre for Theoretical Physics Asia-Pacific,
Beijing/Hangzhou, China  Dian-Jun Zhao  University of Chinese Academy of
Sciences, School of Physical Sciences, Beijing 100049, China  CAS Key
Laboratory of Theoretical Physics, Institute of Theoretical Physics, Chinese
Academy of Sciences, Beijing 100190, China

###### Abstract

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.39(14)(19)$ MeV,
$m_{d}=4.77(16)(09)$ MeV, and $m_{s}=97.8(2.1)(4.0)$ 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.4(3.0)(0.2)$ MeV and the pion decay constant
$F=86.1(0.7)(1.4)$ MeV in the $N_{f}=2$ chiral limit, and the next-to-leading
order low energy constants $\ell_{3}=2.54(46)(5)$ and $\ell_{4}=4.30(7)(10)$.

##  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.

Due to the infamous fermion doubling problem which prevents a straightforward
discretization of the continuum Dirac fermion action, an accurate
determination of the light quark masses is highly non-trivial. Since the
widely used clover fermions suffer the additional chiral symmetry breaking
which induces power divergence with loop corrections, most of the light quark
mass determinations are made with either the staggered fermion (or its
improved versions) Aubin  _et al._ (2004a, b); Mason  _et al._ (2006); Bazavov
_et al._ (2009); Davies  _et al._ (2010); Bazavov  _et al._ (2010); Laiho and
Van de Water (2011); McNeile  _et al._ (2010); Basak  _et al._ (2016); Maezawa
and Petreczky (2016); Bazavov  _et al._ (2018a, b); Basak  _et al._ (2019)
which suffers the mixing between four equivalent “tastes” of a given flavor,
or Ginsburg-Wilson fermion actions like Domain Wall Allton  _et al._ (2008);
Aoki  _et al._ (2011); Blum  _et al._ (2010); Arthur  _et al._ (2013); Blum
_et al._ (2016) or Overlap Yang  _et al._ (2015a) fermion which requires
${\cal O}(10–100)$ times more cost of computational resources than the clover
fermions.

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.
Ishikawa  _et al._ (2008) 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 Sint
(2011), 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_{\pi}L\sim$ 2 gives $m_{u,d}=3.12(24)(08)$ MeV Aoki  _et al._ (2009,
2010a, 2010b, 2012) at ${\overline{\text{MS}}}\rm{(2\ GeV)}$ which is 10%
lower than the present lattice average value $m_{u,d}=3.381(40)$MeV Aoki  _et
al._ (2022) with large uncertainty.

A more systematic study using the SF scheme was conducted by the ALPHA
collaboration with multiple lattice spacing $a\in[0.05,0.086]$ fm but
relatively heavy quark masses $m_{\pi}\geq$ 200 MeV, and their determination
resulted in $m_{u,d}=3.54(12)(9)$MeV Bruno  _et al._ (2020). So far, the most
precise determination of $m_{u,d}=3.469(47)(48)$Durr  _et al._ (2011); Bruno
_et al._ (2020) with clover fermion comes from the BMW collaboration, which
utilized multiple lattice spacings $a\in[0.05,0.012]$ fm with the lightest
pion mass $m_{\pi}=131(2)$ MeV and renormalized the quark mass using the
widely-used RI/MOM scheme Martinelli  _et al._ (1995).

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 higher-order corrections. Thus, the
SMOM scheme Aoki  _et al._ (2008); Sturm  _et al._ (2009) was proposed to
suppress this uncertainty and has been employed in most recent quark mass
determinations using chiral fermions. Nevertheless, a recent study Hasan  _et
al._ (2019) at multiple lattice spacings shows that using either the RI/MOM or
SMOM intermediate scheme can result in the renormalized scalar current under
$\overline{\textrm{MS}}$ scheme differing by 30% at $a\sim 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 Carrasco  _et al._ (2014);
Giusti  _et al._ (2017); Alexandrou  _et al._ (2021) 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 $\sim 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,

$\displaystyle
S_{g}=\frac{1}{N_{c}}\mathrm{Re}\sum_{x,\mu<\nu}\mathrm{Tr}\big{[}1-\hat{\beta}\big{(}\mathcal{P}^{U}_{\mu,\nu}(x)+\frac{c_{1}\mathcal{R}^{U}_{\mu,\nu}(x)}{1-8c_{1}^{0}}\big{)}],$
(1)

where $N_{c}=3$,
$\mathcal{P}^{U}_{\mu,\nu}(x)=U_{\mu}(x)U_{\nu}(x+a\hat{\mu})U^{\dagger}_{\mu}(x+a\hat{\nu})U^{\dagger}_{\nu}(x)$,
$\mathcal{R}^{U}_{\mu,\nu}(x)=U_{\mu}(x)U_{\mu}(x+a\hat{\mu})U_{\nu}(x+2a\hat{\mu})U^{\dagger}_{\mu}(x+a\hat{\mu}+a\hat{\nu})U^{\dagger}_{\mu}(x+a\hat{\nu})U^{\dagger}_{\nu}(x)$,
$U_{\mu}(x)=P[\mathrm{exp}({\rm
i}g_{0}\int_{x+\hat{\mu}a}^{x}\mathrm{d}yA_{\mu}(y))]$,
$\hat{\beta}=(1-8c_{1}^{0})\frac{6}{g_{0}^{2}u_{0}^{4}}\equiv
10/(g_{0}^{2}u_{0}^{4})$ with $c_{1}^{0}=-\frac{1}{12}$,
$c_{1}=\frac{c_{1}^{0}}{u_{0}^{2}}$,
$u_{0}=\langle\frac{\mathrm{Re}\mathrm{Tr}\sum_{x,\mu<\nu}\mathcal{P}^{U}_{\mu\nu}(x)}{6N_{c}\tilde{V}}\rangle^{1/4}$
is the tadpole improvement factor, $\tilde{V}=\tilde{L}^{3}\times\tilde{T}$ is
the dimensionless 4-D volume of the lattice, and we use $\tilde{O}$ for the
dimensionless value of any quantity $O$.

The TITLC fermion action uses 1-step stout smeared link $V$ with smearing
parameter $\rho=0.125$,

$\displaystyle
S_{q}(m)=\sum_{x,\mu=1,...,4,\eta=\pm}\bar{\psi}(x+\eta\hat{\mu}a)\sum\frac{1+\eta\gamma_{\mu}}{2}V_{\eta\mu}(x)\psi(x)$
$\displaystyle\quad+\sum_{x}\psi(x)\big{[}-(4+ma)\delta_{y,x})+c_{\rm
sw}\sigma^{\mu\nu}g_{0}F^{V}_{\mu\nu}\big{]}\psi(x),$ (2)

where $c_{\rm sw}=\frac{1}{v^{3}_{0}}$ with
$v_{0}=\langle\frac{\mathrm{Re}\mathrm{Tr}\sum_{x,\mu<\nu}\mathcal{P}^{V}_{\mu\nu}(x)}{6N_{c}\tilde{V}}\rangle^{1/4}$,
and

$\displaystyle F^{V}_{\mu\nu}$ $\displaystyle=$
$\displaystyle\frac{i}{8a^{2}g_{0}}(\mathcal{P}^{V}_{\mu,\nu}-\mathcal{P}^{V}_{\nu,\mu}+\mathcal{P}^{V}_{\nu,-\mu}-\mathcal{P}^{V}_{-\mu,\nu}$
(3)
$\displaystyle+\mathcal{P}^{V}_{-\mu,-\nu}-\mathcal{P}^{V}_{-\nu,-\mu}+\mathcal{P}^{V}_{-\nu,\mu}-\mathcal{P}^{V}_{\mu,-\nu}).$

Table 1:  Lattice size $\tilde{L}^{3}\times\tilde{T}$, gauge coupling
$\hat{\beta}=10/(g^{2}u_{0}^{4})$, dimensionless bare quark mass parameters
$\tilde{m}^{\rm b}_{l,s}$, renormalized quark masses $m^{R}_{l,s}$ and the
corresponding pseudoscalar mass $m_{\pi,K}$, and the statistics information.

| C24P34 | C24P29 | C32P29 | C32P23 | C48P23 | C48P14 | F32P30 | F48P30 | F32P21 | F48P21 | H48P32
---|---|---|---|---|---|---|---|---|---|---|---
$\tilde{L}^{3}\times\tilde{T}$ | $24^{3}\times 64$ | $24^{3}\times 72$ | $32^{3}\times 64$ | $32^{3}\times 64$ | $48^{3}\times 96$ | $48^{3}\times 96$ | $32^{3}\times 96$ | $48^{3}\times 96$ | $32^{3}\times 64$ | $48^{3}\times 96$ | $48^{3}\times 144$
$\hat{\beta}$ | 6.20 | 6.41 | 6.72
$a$ (fm) | 0.10530(18) | 0.07746(18) | 0.05187(26)
$\tilde{m}^{b}_{l}$ | -0.2770 | -0.2770 | -0.2770 | -0.2790 | -0.2790 | -0.2825 | -0.2295 | -0.2295 | -0.2320 | -0.2320 | -0.1850
$\tilde{m}^{b}_{s}$ | -0.2310 | -0.2400 | -0.2400 | -0.2400 | -0.2400 | -0.2310 | -0.2050 | -0.2050 | -0.2050 | -0.2050 | -0.1700
$m^{R}_{l}$ (MeV) | 22.90(19) | 16.94(12) | 17.35(11) | 10.55(11) | 10.27(10) | 3.638(83) | 18.54(12) | 18.511(92) | 8.58(16) | 8.59(08) | 19.42(05)
$m^{R}_{s}$ (MeV) | 111.41(16) | 87.46(10) | 88.16(10) | 84.48(07) | 84.79(04) | 103.15(05) | 93.23(11) | 93.05(08) | 89.75(10) | 90.43(08) | 95.61(04)
$m_{\pi}$ (MeV) | 341.1(1.8) | 292.7(1.2) | 292.4(1.1) | 228.0(1.2) | 225.6(0.9) | 135.5(1.6) | 303.2(1.3) | 303.4(0.9) | 210.9(2.2) | 207.2(1.1) | 317.2(0.9)
$m_{K}$ (MeV) | 582.7(1.6) | 509.4(1.1) | 509.0(1.1) | 484.1(1.0) | 484.1(1.3) | 510.0(1.0) | 524.6(1.8) | 523.6(1.4) | 492.0(1.7) | 493.0(1.4) | 536.1(3.0)
$n_{\rm cfg}$ | 200 | 476 | 198 | 400 | 62 | 203 | 206 | 99 | 194 | 98 | 176
$n_{\rm src}$ | 32 | 3 | 3 | 3 | 3 | 48 | 3 | 3 | 3 | 12 | 12

The parameters utilized for the simulation, encompassing the lattice size
($\tilde{L}^{3}\times\tilde{T}$), gauge coupling ($\hat{\beta}$), and the
lattice spacing ($a$) determined through the gradient flow Lüscher (2010) with
$w_{0}$ Borsanyi  _et al._ (2012) using the Symanzik action, are outlined in
Table 1. The dimensionless bare degenerated light and strange quark mass
($\tilde{m}^{\rm b}_{l,s}$), renormalized quark masses ($m_{l,s}^{R}$) at
$\overline{\textrm{MS}}$ 2GeV, and the respective pion and kaon masses
($m_{\pi,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 ref . 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_{q}^{R}$ from the bare quark mass parameter $\tilde{m}^{b}_{q}$ can be
subtle since the critical quark mass $\tilde{m}_{\rm crti}$ vanishing the pion
mass is non-zero. A more practical solution defines it through the PCAC
relation Ishikawa  _et al._ (2008):

$Z_{A}\partial_{\mu}A_{\mu}=2m^{R}_{q}Z_{P}P,$ (4)

where $A_{\mu}=\bar{\psi}\gamma_{5}\gamma_{\mu}\psi$ and
$P=\bar{\psi}\gamma_{5}\psi$. The PCAC quark mass $m_{q}^{\rm PC}$ is then
defined through the pion correlation functions:

$\displaystyle m_{q}^{\rm PC}$ $\displaystyle=\frac{m_{\rm
PS}\sum_{\vec{x}}\langle
A_{4}(\vec{x},t)P^{\dagger}(\vec{0},0)\rangle}{2\sum_{\vec{x}}\langle
P(\vec{x},t)P^{\dagger}(\vec{0},0)\rangle}|_{t\rightarrow\infty},$ (5)

where $m_{\mathrm{PS}}$ is the pseudoscalar meson mass. The renormalized quark
mass is subsequently defined as $m^{R}_{q}=Z_{A}/Z_{P}m^{\rm PC}_{q}$.

Figure 1: The dimensionless PCAC quark mass $\tilde{m}_{q}^{\rm PC}=m_{q}^{\rm
PC}a$ v.s. the bare quark mass $\tilde{m}_{q}^{\rm b}=m_{q}^{\rm b}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.

In Fig. 1, we plot the dimensionless PCAC quark mass $\tilde{m}_{q}^{\rm
PC}\equiv m_{q}^{\rm PC}a$ as a function of the dimensionless input bare quark
mass parameters $\tilde{m}_{q}^{b}\equiv m_{q}^{\rm b}a$, at three lattice
spacings with $m_{\pi}\sim 300$ MeV. The figure also includes linear fits
using the following form:

$\displaystyle\tilde{m}^{\rm PC}_{q}=k_{m}(\tilde{m}^{b}_{q}-\tilde{m}_{\rm
crti}),$ (6)

where $\tilde{m}_{\rm crti}$ corresponds to the critical pion mass that makes
the pion mass and $\tilde{m}^{\rm PC}$ vanish. The parameter $k_{m}=1+{\cal
O}(a^{2},\alpha_{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 ${\cal O}(a^{2})$ discretization error and
${\cal O}(\alpha_{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 ${\overline{\text{MS}}}$ scheme, can only be obtained
through regularization-independent (RI) schemes such as RI/MOM Martinelli  _et
al._ (1995) or SMOM Aoki  _et al._ (2008); Sturm  _et al._ (2009). 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 compare to other systematic
uncertainties He  _et al._ (2022).

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_{\pi}$ at three lattice
spacing with $m_{\pi}=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.

Figure 2:  Renormalized scalar matrix element $g_{S,\pi}$ with $m_{\pi}$=317
MeV at three lattice spacing. $g_{S,\pi,{\rm ME}}$ used
$Z^{{\overline{\text{MS}}}\rm{(2\ GeV)}}_{S}$, $g_{S,\pi,{\rm FH}}$ used
$Z^{{\overline{\text{MS}}}\rm{(2\ GeV)}}_{P}$, through either RI/MOM or SMOM
scheme. The extrapolated values deviate by $\sim 7$%.

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,

$\displaystyle g_{S,\pi,{\rm ME}}={\langle\pi|S|\pi\rangle}_{\rm
val}/{\langle\pi|\pi\rangle},$ (7)

where $S=\bar{\psi}\psi$. We show the renormalized $g_{S,\pi,{\rm
ME}}^{{\overline{\text{MS}}}\rm{(2\ GeV)}}=Z^{{\overline{\text{MS}}}\rm{(2\
GeV)}}_{S}g_{S,\pi}$ obtained from RI/MOM and SMOM schemes in Fig. 2, for
$m_{\pi}=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,\pi}$ from the
quark mass dependence of $m_{\pi}$, as

$\displaystyle g_{S,\pi,{\rm FH}}=\frac{1}{2}\frac{\partial
m_{\pi}(m_{q})}{\partial m_{q}}\simeq\frac{m_{\pi}}{4m_{q}}+{\cal
O}(m_{q},a^{2}),$ (8)

where the factor $\frac{1}{2}$ in front of $\frac{\partial
m_{\pi}(m_{q})}{\partial m_{q}}$ is used to average the contribution from two
propagators in the pion correlator. Using the renormalized quark mass
$m_{q}^{{\overline{\text{MS}}}\rm{(2\ GeV)}}$ extracted with RI/MOM or SMOM
scheme, $g_{S,\pi,{\rm FH}}^{{\overline{\text{MS}}}\rm{(2\ GeV)}}$ (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,\pi,{\rm 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,\pi,{\rm
FH}}$ is consistent with the $g_{S,\pi,{\rm ME}}$ using the RI/MOM scheme at
each lattice spacing within two sigma, but have significant difference from
the $g_{S,\pi,{\rm ME}}$ using the SMOM scheme. Thus the deviation between the
$g_{S,\pi,{\rm FH}}$ and $g_{S,\pi,{\rm 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 ref , 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 1, 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_{\pi^{0}}$ of 134.98 MeV within
1% statistical uncertainty. The charged pion mass $m_{\pi^{\pm}}$ 139.57 MeV
receives the QED correction 4.53(6) MeV Feng  _et al._ (2022) and then the
subtracted pure QCD $m_{\pi^{\pm}}$ is consistent with that of $m_{\pi^{0}}$
within the uncertainty.

The corresponding renormalized light quark mass and pion decay constants can
also be determined as:

$\displaystyle m_{l}^{{\overline{\text{MS}}}\rm{(2\
GeV)}}(a=0.105~{}\mathrm{fm})$ $\displaystyle=3.64(8)(11)~{}\mathrm{MeV},$
$\displaystyle f_{\pi}(a=0.105~{}\mathrm{fm})$
$\displaystyle=121.9(5)~{}\mathrm{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){\rm MeV}$ Zyla  _et al._ (2020) 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_{\pi}L>3.5$, on each of the 11
ensembles. Then we use the NLO partially quenched $\chi$PT form Sharpe (1997)
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.

Table 2: Summary of our determination on quark masses at ${\overline{\text{MS}}}{\rm(2GeV)}$ and the other quantities , through the intermediate RI/MOM or SMOM schemes, with comparison with FLAG Aoki  _et al._ (2022) and/or PDG Zyla  _et al._ (2020). The difference between two schemes is considered as systematic uncertainty in the combined determination. |  $m_{l}$ (MeV) |  $m_{u}$ (MeV) |  $m_{d}$ (MeV) |  $m_{s}$ (MeV) |  | $m_{s}/m_{l}$ | $m_{u}/m_{d}$
---|---|---|---|---|---|---|---
RI/MOM | 3.571(94) | 2.39(14) | 4.77(16) | 97.8(2.1) |  | 27.35(36) | 0.501(35)
SMOM | 3.441(46) | 2.20(08) | 4.68(08) | 93.8(1.1) |  | 27.24(23) | 0.471(21)
Combined | 3.57(9)(13) | 2.39(14)(19) | 4.77(16)(09) | 97.8(2.1)(4.0) |  | 27.35(36)(13) | 0.501(35)(30)
FLAG/PDG | 3.381(40) | 2.27(09) | 4.67(09) | 92.2(1.0) |  | 27.42(12) | 0.485(19)
|  $\Sigma^{1/3}$(MeV) |  $F$ (MeV) | $F_{\pi}/F$ |  $f_{\pi}$(MeV) |  $f_{K^{\pm}}$(MeV) | $\ell_{3}$ | $\ell_{4}$
RI/MOM | 268.4(3.0) | 86.1(7) | 1.071(9) | 130.1(1.1) | 153.8(0.7) | 2.54(46) | 4.30(07)
SMOM | 268.6(1.8) | 84.7(6) | 1.089(8) | 128.0(0.9) | 151.1(0.6) | 2.59(20) | 4.20(05)
Combined | 268.4(3.0)(0.2) | 86.1(0.7)(1.4) | 1.071(9)(18) | 130.1(1.1)(2.1) | 153.8(0.7)(2.7) | 2.54(46)(05) | 4.30(07)(10)
FLAG/PDG | 272(5) | 86.8(6) | 1.062(7) | 130.2(0.8) | 155.7(0.7) | 3.07(64) | 4.02(45)

Figure 3: The corrected unitary $(m^{\rm uni}_{\pi})^{2}/m_{q}$ and the decay
constant $f^{\rm uni}_{\pi}$ varies with the quark mass, at three lattice
spacing (colored data points and bands) 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_{\pi})^{2}/m^{{\overline{\text{MS}}}\rm{(2\ GeV)}}_{q}$ (upper panel) and
also $f_{\pi}$ (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_{\pi}$ to be obviously higher, while the impact on the
$(m_{\pi}^{2}/m_{q})$ ratio and, consequently, $m_{q}$ is much weaker.

The $m_{\pi}$ and $f_{\pi}$ with unitary valence and sea quark masses have the
following parameterization,

$\displaystyle m_{\pi}^{2}$
$\displaystyle=\Lambda_{\chi}^{2}2y[1+y(\ln\frac{2y\Lambda^{2}_{\chi}}{m^{2}_{\pi,{\rm
phys}}}-{\ell}_{3})+\mathcal{O}(y^{2})],$ (10) $\displaystyle F_{\pi}$
$\displaystyle=F[1-2y(\ln\frac{2y\Lambda^{2}_{\chi}}{m^{2}_{\pi,{\rm
phys}}}-\ell_{4})+\mathcal{O}(y^{2})],$ (11)

where $\Lambda_{\chi}=4\pi F$, $y=\frac{\Sigma
m_{l}}{F^{2}\Lambda_{\chi}^{2}}$. $\Sigma$, $F$ and $\ell_{3,4}$ are low
energy constants. Our determination of those constants are also collected in
Table 12, consistent with the current $N_{f}=2+1$ FLAG average but have
smaller uncertainties except $F$.

In this work, we use the $m_{K^{\pm}}$ and $m_{K^{0}}$ with the constraint
$m^{\rm phys}_{u}+m^{\rm phys}_{d}=2m^{\rm 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 Giusti  _et al._ (2017). On each
ensemble, we calculate the strange quark propagators with a unitary strange
quark mass $m^{\rm v}_{s}=m^{\rm s}_{s}$, and also 2 partially quenched quark
masses $m^{\rm v}_{s}\sim 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\times 3$ partially quenched Kaon masses and decay
constants on all the ensembles are fitted with the following form proposed in
a recent work Alexandrou  _et al._ (2021),

$\displaystyle m_{K}^{2}(m^{\rm v}_{l},m^{\rm s}_{l},m^{\rm v}_{s},m^{\rm
s}_{s},a)=(b_{s}^{\rm v}m^{\rm v}_{s}+b_{s}^{\rm s}m^{\rm s}_{s}+b_{l}^{\rm
v}m^{\rm v}_{l}+b_{l}^{\rm s}m^{\rm s}_{l})$
$\displaystyle~{}~{}~{}~{}~{}~{}\times\left[1+c^{K}_{l}m^{\rm
v}_{l}+c^{K}_{m}a^{2}+c_{L}^{K}\exp{(-m_{K}L)}\right].$ (12)

Based on the fit of $m_{K}^{2}$, the total strange quark mass dependence
$b_{s}^{\rm v}+b_{s}^{\rm s}=2.40(06)$ is consistent with the leading order
light quark mass dependence $b_{l}^{\rm v}+b_{l}^{\rm s}=2.14(89)$, and the
coefficient of the non-linear quark mass dependence $c_{l}^{K}=2.7(3.4)$ can
not be determined based on current statistics.

Figure 4: The corrected kaon mass $m^{\rm cr}_{K}$ and the decay constant
$f^{\rm cr}_{K}$ with the physical light quark mass $m_{l}^{\rm phys}$, varies
with the strange quark mass at three lattice spacing (colored data points and
bands) and also continuum (gray band).

In Fig. 4, we show the corrected kaon mass $m^{\rm cr}_{K}$ and decay constant
$f^{\rm cr}_{K}$ with the light quark mass $m_{l}$ corrected to its physical
value $m_{l}^{\rm phys}$. 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_{\pi}$ 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 $\chi^{2}$/d.o.f. More
information on the global fit can be found in the supplemental materials ref .

The physical quark masses $m_{u,d,s}$ and also corresponding $f_{\pi,K}$ using
$m_{l}^{\rm phys}$ and intermediate RI/MOM scheme, are collected in Table 12.
In addition, Table 12 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 $\sim$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. With this systematic uncertainty, all
our determinations are consistent with the present lattice averages Aoki  _et
al._ (2022) and/or PDG Zyla  _et al._ (2020) 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
tadpole-improved clover and Symanzik actions. The major results are summarized
in Table 2.

Similar to one of the most precise works Durr  _et al._ (2011) with the clover
fermion, we have skipped the axial current improvement Della Morte  _et al._
(2005) since the improvement coefficient itself strongly depends on the
lattice spacing $a$ and bring the improvement term to be consistent with an
${\cal O}(a^{2})$ correction. As evidence, both $f_{\pi}$ and $m_{q}$ show
good consistency with a simple linear $a^{2}$ lattice spacing dependence. Then
assigning $c_{A}\sim 0.05a/(0.105~{}\mathrm{fm})$ in the improved axial vector
current $A^{\rm imp}_{\mu}=A_{\mu}+c_{A}a\partial_{\mu}P$ Della Morte  _et
al._ (2005) can eliminate the discretization error in $f_{\pi,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^{\rm PC}_{q}$ and
$g_{S}^{R}=Z_{S}g_{S}$ to be smaller. Our final prediction of the quark masses
are $5.6(2.8)$% higher than the 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,\pi}\simeq\frac{m_{\pi}}{4m_{q}}$ to be violated
by 7(3)% after the linear ${\cal 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_{\pi,K}$ to be
2-3% smaller than the physical values after the linear ${\cal O}(a^{2})$
continuum extrapolation. As a result, 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. This indicates 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.

It is worth mentioning that in Ref. Durr  _et al._ (2011), the trace-
subtraction trick $\bar{S}=S-\frac{1}{4}\mathrm{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 generated.

## Conflict of interest

The authors declare that they have no conflict of interest.

## Acknowledgment

We thank Ying Chen, Hengtong Ding, Chuan Liu, Zhaofeng Liu, Qi-An Zhang, and
the other CLQCD members for valuable comments and suggestions. The
calculations were performed using the Chroma software suite Edwards and Joo
(2005) with QUDA Clark  _et al._ (2010); Babich  _et al._ (2011); Clark  _et
al._ (2016) through HIP programming model Bi  _et al._ (2020). The numerical
calculation were carried out on the ORISE Supercomputer, HPC Cluster of ITP-
CAS, the Southern Nuclear Science Computing Center(SNSC), the Siyuan-1 cluster
supported by the Center for High Performance Computing at Shanghai Jiao Tong
University, and Advanced Computing East China Sub-center. This work is
supported in part by the science and education integration young faculty
project of University of Chinese Academy of Sciences, the Strategic Priority
Research Program of Chinese Academy of Sciences, Grant No. XDB34030303 and
XDPB15, NSFC grants No. 12293062, 12293065 and 12047503, and also a NSFC-DFG
joint grant under Grant No. 12061131006 and SCHA 458/22.

## References

  * Aubin  _et al._ (2004a) C. Aubin, C. Bernard, C. E. DeTar, J. Osborn, S. Gottlieb, E. B. Gregory, D. Toussaint, U. M. Heller, J. E. Hetrick, and R. Sugar (MILC), Phys. Rev. D 70, 114501 (2004a), arXiv:hep-lat/0407028 .
  * Aubin  _et al._ (2004b) C. Aubin _et al._ (HPQCD, MILC, UKQCD), Phys. Rev. D 70, 031504 (2004b), arXiv:hep-lat/0405022 .
  * Mason  _et al._ (2006) Q. Mason, H. D. Trottier, R. Horgan, C. T. H. Davies, and G. P. Lepage (HPQCD), Phys. Rev. D 73, 114501 (2006), arXiv:hep-ph/0511160 .
  * Bazavov  _et al._ (2009) A. Bazavov _et al._ (MILC), PoS CD09, 007 (2009), arXiv:0910.2966 [hep-ph] .
  * Davies  _et al._ (2010) C. T. H. Davies, C. McNeile, K. Y. Wong, E. Follana, R. Horgan, K. Hornbostel, G. P. Lepage, J. Shigemitsu, and H. Trottier, Phys. Rev. Lett. 104, 132003 (2010), arXiv:0910.3102 [hep-ph] .
  * Bazavov  _et al._ (2010) A. Bazavov _et al._ (MILC), Rev. Mod. Phys. 82, 1349 (2010), arXiv:0903.3598 [hep-lat] .
  * Laiho and Van de Water (2011) J. Laiho and R. S. Van de Water, PoS LATTICE2011, 293 (2011), arXiv:1112.4861 [hep-lat] .
  * McNeile  _et al._ (2010) C. McNeile, C. T. H. Davies, E. Follana, K. Hornbostel, and G. P. Lepage, Phys. Rev. D 82, 034512 (2010), arXiv:1004.4285 [hep-lat] .
  * Basak  _et al._ (2016) S. Basak _et al._ (MILC), PoS LATTICE2015, 259 (2016), arXiv:1606.01228 [hep-lat] .
  * Maezawa and Petreczky (2016) Y. Maezawa and P. Petreczky, Phys. Rev. D 94, 034507 (2016), arXiv:1606.08798 [hep-lat] .
  * Bazavov  _et al._ (2018a) A. Bazavov _et al._ , Phys. Rev. D 98, 074512 (2018a), arXiv:1712.09262 [hep-lat] .
  * Bazavov  _et al._ (2018b) A. Bazavov _et al._ (Fermilab Lattice, MILC, TUMQCD), Phys. Rev. D 98, 054517 (2018b), arXiv:1802.04248 [hep-lat] .
  * Basak  _et al._ (2019) S. Basak _et al._ (MILC), Phys. Rev. D 99, 034503 (2019), arXiv:1807.05556 [hep-lat] .
  * Allton  _et al._ (2008) C. Allton _et al._ (RBC-UKQCD), Phys. Rev. D 78, 114509 (2008), arXiv:0804.0473 [hep-lat] .
  * Aoki  _et al._ (2011) Y. Aoki _et al._ (RBC, UKQCD), Phys. Rev. D 83, 074508 (2011), arXiv:1011.0892 [hep-lat] .
  * Blum  _et al._ (2010) T. Blum, R. Zhou, T. Doi, M. Hayakawa, T. Izubuchi, S. Uno, and N. Yamada, Phys. Rev. D 82, 094508 (2010), arXiv:1006.1311 [hep-lat] .
  * Arthur  _et al._ (2013) R. Arthur _et al._ (RBC, UKQCD), Phys. Rev. D 87, 094514 (2013), arXiv:1208.4412 [hep-lat] .
  * Blum  _et al._ (2016) T. Blum _et al._ (RBC, UKQCD), Phys. Rev. D 93, 074505 (2016), arXiv:1411.7017 [hep-lat] .
  * Yang  _et al._ (2015a) Y.-B. Yang _et al._ , Phys. Rev. D92, 034517 (2015a), arXiv:1410.3343 [hep-lat] .
  * Ishikawa  _et al._ (2008) T. Ishikawa _et al._ (JLQCD), Phys. Rev. D 78, 011502 (2008), arXiv:0704.1937 [hep-lat] .
  * Sint (2011) S. Sint, Nucl. Phys. B 847, 491 (2011), arXiv:1008.4857 [hep-lat] .
  * Aoki  _et al._ (2009) S. Aoki _et al._ (PACS-CS), Phys. Rev. D 79, 034503 (2009), arXiv:0807.1661 [hep-lat] .
  * Aoki  _et al._ (2010a) S. Aoki _et al._ (PACS-CS), Phys. Rev. D 81, 074503 (2010a), arXiv:0911.2561 [hep-lat] .
  * Aoki  _et al._ (2010b) S. Aoki _et al._ (PACS-CS), JHEP 08, 101 (2010b), arXiv:1006.1164 [hep-lat] .
  * Aoki  _et al._ (2012) S. Aoki _et al._ , Phys. Rev. D 86, 034507 (2012), arXiv:1205.2961 [hep-lat] .
  * Aoki  _et al._ (2022) Y. Aoki _et al._ (Flavour Lattice Averaging Group (FLAG)), Eur. Phys. J. C 82, 869 (2022), arXiv:2111.09849 [hep-lat] .
  * Bruno  _et al._ (2020) M. Bruno, I. Campos, P. Fritzsch, J. Koponen, C. Pena, D. Preti, A. Ramos, and A. Vladikas (ALPHA), Eur. Phys. J. C 80, 169 (2020), arXiv:1911.08025 [hep-lat] .
  * Durr  _et al._ (2011) S. Durr, Z. Fodor, C. Hoelbling, S. D. Katz, S. Krieg, T. Kurth, L. Lellouch, T. Lippert, K. K. Szabo, and G. Vulvert (BMW), JHEP 08, 148 (2011), arXiv:1011.2711 [hep-lat] .
  * Martinelli  _et al._ (1995) G. Martinelli, C. Pittori, C. T. Sachrajda, M. Testa, and A. Vladikas, Nucl. Phys. B445, 81 (1995), arXiv:hep-lat/9411010 [hep-lat] .
  * Aoki  _et al._ (2008) Y. Aoki _et al._ , Phys. Rev. D78, 054510 (2008), arXiv:0712.1061 [hep-lat] .
  * Sturm  _et al._ (2009) C. Sturm, Y. Aoki, N. H. Christ, T. Izubuchi, C. T. C. Sachrajda, and A. Soni, Phys. Rev. D80, 014501 (2009), arXiv:0901.2599 [hep-ph] .
  * Hasan  _et al._ (2019) N. Hasan, J. Green, S. Meinel, M. Engelhardt, S. Krieg, J. Negele, A. Pochinsky, and S. Syritsyn, Phys. Rev. D 99, 114505 (2019), arXiv:1903.06487 [hep-lat] .
  * Carrasco  _et al._ (2014) N. Carrasco _et al._ (European Twisted Mass), Nucl. Phys. B 887, 19 (2014), arXiv:1403.4504 [hep-lat] .
  * Giusti  _et al._ (2017) D. Giusti, V. Lubicz, C. Tarantino, G. Martinelli, F. Sanfilippo, S. Simula, and N. Tantalo, Phys. Rev. D 95, 114504 (2017), arXiv:1704.06561 [hep-lat] .
  * Alexandrou  _et al._ (2021) C. Alexandrou _et al._ (Extended Twisted Mass), Phys. Rev. D 104, 074515 (2021), arXiv:2104.13408 [hep-lat] .
  * Lüscher (2010) M. Lüscher, JHEP 08, 071 (2010), [Erratum: JHEP 03, 092 (2014)], arXiv:1006.4518 [hep-lat] .
  * Borsanyi  _et al._ (2012) S. Borsanyi _et al._ (BMW), JHEP 09, 010 (2012), arXiv:1203.4469 [hep-lat] .
  * (38) Supplementary materials .
  * He  _et al._ (2022) F. He, Y.-J. Bi, T. Draper, K.-F. Liu, Z. Liu, and Y.-B. Yang ($\chi$QCD), Phys. Rev. D 106, 114506 (2022), arXiv:2204.09246 [hep-lat] .
  * Feng  _et al._ (2022) X. Feng, L. Jin, and M. J. Riberdy, Phys. Rev. Lett. 128, 052003 (2022), arXiv:2108.05311 [hep-lat] .
  * Zyla  _et al._ (2020) P. A. Zyla _et al._ (Particle Data Group), PTEP 2020, 083C01 (2020).
  * Sharpe (1997) S. R. Sharpe, Phys. Rev. D 56, 7052 (1997), [Erratum: Phys.Rev.D 62, 099901 (2000)], arXiv:hep-lat/9707018 .
  * Della Morte  _et al._ (2005) M. Della Morte, R. Hoffmann, and R. Sommer, JHEP 03, 029 (2005), arXiv:hep-lat/0503003 .
  * Edwards and Joo (2005) R. G. Edwards and B. Joo (SciDAC, LHPC, UKQCD),  _Lattice field theory. Proceedings, 22nd International Symposium, Lattice 2004, Batavia, USA, June 21-26, 2004_ , Nucl. Phys. Proc. Suppl. 140, 832 (2005), [,832(2004)], arXiv:hep-lat/0409003 [hep-lat] .
  * Clark  _et al._ (2010) M. A. Clark, R. Babich, K. Barros, R. C. Brower, and C. Rebbi, Comput. Phys. Commun. 181, 1517 (2010), arXiv:0911.3191 [hep-lat] .
  * Babich  _et al._ (2011) R. Babich, M. A. Clark, B. Joo, G. Shi, R. C. Brower, and S. Gottlieb, in  _SC11 International Conference for High Performance Computing, Networking, Storage and Analysis Seattle, Washington, November 12-18, 2011_ (2011) arXiv:1109.2935 [hep-lat] .
  * Clark  _et al._ (2016) M. A. Clark, B. Jo, A. Strelchenko, M. Cheng, A. Gambhir, and R. Brower, (2016), arXiv:1612.07873 [hep-lat] .
  * Bi  _et al._ (2020) Y.-J. Bi, Y. Xiao, W.-Y. Guo, M. Gong, P. Sun, S. Xu, and Y.-B. Yang,  _Proceedings, 37th International Symposium on Lattice Field Theory (Lattice 2019): Wuhan, China, June 16-22 2019_ , PoS LATTICE2019, 286 (2020), arXiv:2001.05706 [hep-lat] .
  * Borsanyi  _et al._ (2021) S. Borsanyi _et al._ , Nature 593, 51 (2021), arXiv:2002.12347 [hep-lat] .
  * Gracey (2023) J. A. Gracey, Eur. Phys. J. C 83, 181 (2023), arXiv:2210.12420 [hep-ph] .
  * Zhang  _et al._ (2021) K. Zhang, Y.-Y. Li, Y.-K. Huo, A. Schäfer, P. Sun, and Y.-B. Yang ($\chi$QCD), Phys. Rev. D 104, 074501 (2021), arXiv:2012.05448 [hep-lat] .
  * Chang  _et al._ (2021) L. Chang, Y.-B. Liu, K. Raya, J. Rodríguez-Quintero, and Y.-B. Yang, Phys. Rev. D 104, 094509 (2021), arXiv:2105.06596 [hep-lat] .
  * Yang  _et al._ (2015b) Y.-B. Yang, Y. Chen, T. Draper, M. Gong, K.-F. Liu, Z. Liu, and J.-P. Ma, Phys. Rev. D91, 074516 (2015b), arXiv:1405.4440 [hep-ph] .

## Supplemental Materials

###  IV.1 Simulation details

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,

$\displaystyle
S_{g}=\frac{1}{N_{c}}\mathrm{Re}\sum_{x,\mu<\nu}\mathrm{Tr}\big{[}1-\hat{\beta}\big{(}\mathcal{P}^{U}_{\mu,\nu}(x)+\frac{c_{1}\mathcal{R}^{U}_{\mu,\nu}(x)}{1-8c_{1}^{0}}\big{)}],$
(13)

where $N_{c}=3$,
$\mathcal{P}^{U}_{\mu,\nu}(x)=U_{\mu}(x)U_{\nu}(x+a\hat{\mu})U^{\dagger}_{\mu}(x+a\hat{\nu})U^{\dagger}_{\nu}(x)$,
$\mathcal{R}^{U}_{\mu,\nu}(x)=U_{\mu}(x)U_{\mu}(x+a\hat{\mu})U_{\nu}(x+2a\hat{\mu})U^{\dagger}_{\mu}(x+a\hat{\mu}+a\hat{\nu})U^{\dagger}_{\mu}(x+a\hat{\nu})U^{\dagger}_{\nu}(x)$,
$U_{\mu}(x)=P[\mathrm{exp}({\rm
i}g_{0}\int_{x+\hat{\mu}a}^{x}\mathrm{d}yA_{\mu}(y))]$,
$\hat{\beta}=(1-8c_{1}^{0})\frac{6}{g_{0}^{2}u_{0}^{4}}\equiv
10/(g_{0}^{2}u_{0}^{4})$ with $c_{1}^{0}=-\frac{1}{12}$,
$c_{1}=\frac{c_{1}^{0}}{u_{0}^{2}}$,
$u_{0}=\langle\frac{\mathrm{Re}\mathrm{Tr}\sum_{x,\mu<\nu}\mathcal{P}^{U}_{\mu\nu}(x)}{6N_{c}\tilde{V}}\rangle^{1/4}$
is the tadpole improvement factor, $\tilde{V}=\tilde{L}^{3}\times\tilde{T}$ is
the dimensionless 4-D volume of the lattice, and we use $\tilde{O}$ for the
dimensionless value of any quantity $O$

The TITLC fermion action uses the stout smeared link $V$ with smearing
parameter $\rho=0.125$,

$\displaystyle
S_{q}(m)=\sum_{x,\mu=1,...,4,\eta=\pm}\bar{\psi}(x+\eta\hat{\mu}a)\sum\frac{1+\eta\gamma_{\mu}}{2}V_{\eta\mu}(x)\psi(x)$
$\displaystyle\quad+\sum_{x}\psi(x)\big{[}-(4+ma)\delta_{y,x})+c_{\rm
sw}\sigma^{\mu\nu}g_{0}F^{V}_{\mu\nu}\big{]}\psi(x),$ (14)

where $c_{\rm sw}=\frac{1}{v^{3}_{0}}$ with
$v_{0}=\langle\frac{\mathrm{Re}\mathrm{Tr}\sum_{x,\mu<\nu}\mathcal{P}^{V}_{\mu\nu}(x)}{6N_{c}\tilde{V}}\rangle^{1/4}$,
and

$\displaystyle F^{V}_{\mu\nu}$ $\displaystyle=$
$\displaystyle\frac{i}{8a^{2}g_{0}}(\mathcal{P}^{V}_{\mu,\nu}-\mathcal{P}^{V}_{\nu,\mu}+\mathcal{P}^{V}_{\nu,-\mu}-\mathcal{P}^{V}_{-\mu,\nu}$
(15)
$\displaystyle+\mathcal{P}^{V}_{-\mu,-\nu}-\mathcal{P}^{V}_{-\nu,-\mu}+\mathcal{P}^{V}_{-\nu,\mu}-\mathcal{P}^{V}_{\mu,-\nu}).$

Table 3:  Lattice size $\tilde{L}^{3}\times\tilde{T}$, gauge coupling
$\hat{\beta}$, dimensionless bare quark mass parameters $\tilde{m}^{\rm
b}_{l,s}$, renormalized quark masses $m^{R}_{l,s}$ and the corresponding
pseudoscalar mass $m_{\pi,K,\eta_{s}}$, and the statistics information.

| C24P34 | C24P29 | C32P29 | C32P23 | C48P23 | C48P14 | F32P30 | F48P30 | F32P21 | F48P21 | H48P32
---|---|---|---|---|---|---|---|---|---|---|---
$\tilde{L}^{3}\times\tilde{T}$ | $24^{3}\times 64$ | $24^{3}\times 72$ | $32^{3}\times 64$ | $32^{3}\times 64$ | $48^{3}\times 96$ | $48^{3}\times 96$ | $32^{3}\times 96$ | $48^{3}\times 96$ | $32^{3}\times 64$ | $48^{3}\times 96$ | $48^{3}\times 144$
$\hat{\beta}$ | 6.20 | 6.41 | 6.72
$a$ (fm) | 0.10530(18) | 0.07746(18) | 0.05187(26)
$\tilde{m}^{b}_{l}$ | -0.2770 | -0.2770 | -0.2770 | -0.2790 | -0.2790 | -0.2825 | -0.2295 | -0.2295 | -0.2320 | -0.2320 | -0.1850
$\tilde{m}^{b}_{s}$ | -0.2310 | -0.2400 | -0.2400 | -0.2400 | -0.2400 | -0.2310 | -0.2050 | -0.2050 | -0.2050 | -0.2050 | -0.1700
$m^{R}_{l}$ (MeV) | 22.90(19) | 16.94(12) | 17.35(11) | 10.55(11) | 10.27(10) | 3.638(83) | 18.54(12) | 18.511(92) | 8.58(16) | 8.59(08) | 19.42(05)
$m^{R}_{s}$ (MeV) | 111.41(16) | 87.46(10) | 88.16(10) | 84.48(07) | 84.79(04) | 103.15(05) | 93.23(11) | 93.05(08) | 89.75(10) | 90.43(08) | 95.61(04)
$m_{\pi}$ (MeV) | 340.5(1.7) | 292.7(1.2) | 292.4(1.1) | 228.0(1.2) | 225.6(0.9) | 135.5(1.6) | 303.2(1.3) | 303.4(0.9) | 210.9(2.2) | 207.2(1.1) | 317.2(0.9)
$m_{K}$ (MeV) | 582.7(1.6) | 509.4(1.1) | 509.0(1.1) | 484.1(1.0) | 484.1(1.3) | 525.6(1.3) | 524.6(1.8) | 523.6(1.4) | 492.0(1.7) | 493.0(1.4) | 536.1(3.0)
$m_{\eta_{s}}$ (MeV) | 750.7(0.9) | 659.13(0.6) | 659.7(0.7) | 645.1(0.5) | 645.1(0.3) | 708.1(0.3) | 677.7(0.9) | 676.0(0.6) | 661.8(0.7) | 663.9(0.6) | 700.0(0.5)
$n_{\rm cfg}$ | 200 | 476 | 198 | 400 | 62 | 203 | 206 | 99 | 194 | 98 | 176
$n_{\rm src}$ | 32 | 3 | 3 | 3 | 3 | 48 | 3 | 3 | 3 | 12 | 12

The parameters utilized for the simulation, encompassing the lattice size
($\tilde{L}^{3}\times\tilde{T}$), gauge coupling ($\hat{\beta}$), and the
lattice spacing ($a$), determined through the gradient flow, are outlined in
Table 3. Also included in the table are the dimensionless bare degenerated
light and strange quark mass ($\tilde{m}^{\rm b}_{l,s}$), renormalized quark
masses ($m_{l,s}^{R}$) at $\overline{\textrm{MS}}$ 2GeV, and the respective
pion and $\eta_{s}$ masses ($m_{\pi,K,\eta_{s}}$).

The ensembles are named in the following format, “prefix symbol for given
$\hat{\beta}$ \+ Spacial size + P + pion mass in MeV $\times 0.1$”, e.g.,
C48P14 corresponds to the ensemble at a=0.105 fm (corresponds to the symbol
“C”) with size $48^{3}\times 96$ and pion mass $135.5$ MeV. We use ”C/F/H” for
the ensembles with $\hat{\beta}$= 6.20, 6.41 and 6.72, respectively.

Among the ensembles generated, all of the ensembles except F32P21 satisfy the
finite volume criterion $m_{\pi}L>3.5$, where $L=\tilde{L}a$ denotes the
spatial length in the physical unit. For a more comprehensive control over the
finite volume effect, we include three large volume ensembles (C32P29, C48P23,
and F48P30) with $m_{\pi}L>5$, and a small volume ensemble (F32P21) with
$m_{\pi}L\sim 2.5$. A test ensemble (C24P34) is designed as a similar bare
light quark mass parameter $\tilde{m}^{b}_{l}$ (along with $u_{0}$ and
$v_{0}$) as C24P29, but with a slightly heavier strange quark mass, in order
to control the mistuning effect of the strange quark mass.

The rest part of 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.1.1; Based on the determination of an uniform
lattice spacing at given $\hat{\beta}$ detailed in Sec. IV.1.2, the mistuning
effect of the tadpole improvement factors is not always negligible and will be
addressed in Sec. IV.1.3.

####  IV.1.1 Dimensionless joint fit

We construct two kinds of the two-point functions for the meson states:

$\displaystyle C_{2,wp}^{\Gamma\Gamma^{\prime}}(t;m_{q})$
$\displaystyle=\sum_{\vec{x}}\langle\mathrm{Tr}[S^{\dagger}_{C}(\vec{x},y_{0}+t;y_{0};m_{q})\gamma_{5}\Gamma$
$\displaystyle\quad\quad\quad
S_{C}(\vec{x},y_{0}+t;y_{0};m_{q})\Gamma^{\prime}\gamma_{5}]\rangle/L^{3},$
(16) $\displaystyle C_{2,ww}^{\Gamma\Gamma^{\prime}}(t;m_{q})$
$\displaystyle=\sum_{\vec{x},\vec{z}}\langle\mathrm{Tr}[S^{\dagger}_{C}(\vec{x},y_{0}+t;y_{0};m_{q})\gamma_{5}\Gamma$
$\displaystyle\quad\quad\quad
S_{C}(\vec{z},y_{0}+t;y_{0};m_{q})\Gamma^{\prime}\gamma_{5}]\rangle/L^{3},$
(17)

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

$\displaystyle
S_{C}(x;y_{0};m_{q})=\sum_{\vec{y}}S(x;\vec{y},y_{0};m_{q},U_{C}),$ (18)

$S(x;y;m_{q},U)\equiv\psi(x;m_{q},U)\bar{\psi}(y;m_{q},U)$ is the quark
propagator of the quark field $\psi$ 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 condition
$\mathrm{Im}[\sum_{i=1,2,3}(U_{C}(x)-U_{C}(x-a\hat{n}_{i}))]=0$.

For the Clover fermion action, defining the renormalized quark mass
$m_{q}^{R}$ from the bare quark mass parameter $\tilde{m}^{b}_{q}$ can be
subtle since the critical quark mass $\tilde{m}_{\rm crti}$ vanishing the pion
mass is non-zero. A more practical solution defines it through the PCAC
relation Ishikawa  _et al._ (2008):

$Z_{A}\partial_{\mu}A_{\mu}=2m^{R}_{q}Z_{P}P,$ (19)

where $A_{\mu}=\bar{\psi}\gamma_{5}\gamma_{\mu}\psi$ and
$P=\bar{\psi}\gamma_{5}\psi$. The PCAC quark mass $m_{q}$ is then defined
through the pion correlation functions:

$\displaystyle m_{q}^{\rm PC}$ $\displaystyle=\frac{m_{\rm
PS}\sum_{\vec{x}}\langle
A_{4}(\vec{x},t)P^{\dagger}(\vec{0},0)\rangle}{2\sum_{\vec{x}}\langle
P(\vec{x},t)P^{\dagger}(\vec{0},0)\rangle}|_{t\rightarrow\infty}.$ (20)

The renormalized quark mass is subsequently defined as
$m^{R}_{q}=Z_{A}/Z_{P}m^{\rm PC}_{q}$.

Figure 5: The ratios for $\tilde{m}_{q}$ (left top), $\tilde{f}_{\rm PS}$
(right top), $\tilde{m}_{\rm PS}$ (left bottom) defined in Eq. (21-24), and
the correlator $\tilde{C}_{2,wp}^{{\cal PP}}$ (right bottom) as functions of
$\tilde{t}$, for the physical light quark mass at the coarsest lattice
spacing. The joint fit results are shown on the plots as gray bands. Figure 6:
Similar to Fig. 5 but for the strange quark mass and fewer sources ($n_{\rm
src}=3$).

Through a joint fit ($\tilde{O}$ is the dimensionless value of any quantity
$O$),

$\displaystyle\frac{\tilde{C}_{2,wp}^{{\cal A}_{4}{\cal
P}}(\tilde{t}-1)-\tilde{C}_{2,wp}^{{\cal A}_{4}{\cal
P}}(\tilde{t}+1)}{4\tilde{C}_{2,wp}^{{\cal
PP}}(\tilde{t})}|_{0\ll\tilde{t}\ll\tilde{T}}=\frac{\mathrm{Sinh}(\tilde{m}_{\rm
PS})}{\tilde{m}_{\rm PS}}\tilde{m}_{q}^{\rm PC},$ (21)
$\displaystyle\sqrt{\frac{\tilde{C}_{2,wp}^{{\cal
PP}}(\tilde{t})}{\tilde{C}_{2,ww}^{{\cal
PP}}(\tilde{t})}}|_{0\ll\tilde{t}\ll\tilde{T}}=\frac{\tilde{m}_{\rm
PS}^{2}}{2\tilde{m}_{q}^{\rm PC}\sqrt{Z_{wp}}}\tilde{f}_{\rm PS},$ (22)
$\displaystyle\tilde{C}_{2,wp}^{{\cal
PP}}(\tilde{t})|_{0\ll\tilde{t}\ll\tilde{T}}=\frac{Z_{wp}}{2\tilde{m}_{\rm
PS}}(e^{-\tilde{m}_{\rm PS}\tilde{t}}+e^{-\tilde{m}_{\rm
PS}(\tilde{T}-\tilde{t})}),$ (23)

the PCAC quark mass $\tilde{m}_{q}^{\rm PC}$, pseudo-scalar mass

$\displaystyle\tilde{m}_{\rm
PS}=\mathrm{Cosh}^{-1}\frac{\tilde{C}_{2,wp}^{{\cal
PP}}(\tilde{t}-1)+\tilde{C}_{2,wp}^{{\cal
PP}}(\tilde{t}+1)}{2\tilde{C}_{2,wp}^{{\cal
PP}}(\tilde{t})}|_{0\ll\tilde{t}\ll\tilde{T}},$ (24)

and decay constant $\tilde{f}_{\rm PS}$ 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_{\rm src}$=48 of 96 times slides on $n_{\rm
cfg}$=203 configurations. The values of $n_{\rm src}$ and $n_{\rm cfg}$ of the
other ensembles can be found in Table. 1. The ratios (red data points) defined
in Eq. (21) for $\tilde{m}_{q}$ are consistent with a constant in the region
of $\tilde{t}>10$, as shown in the top left panel. The situation is similar
for $\tilde{f}_{\rm PS}$ (top right panel) and also $\tilde{m}_{\rm PS}$
(bottom left panel). Since the statistics on C48P14 are limited, we performed
a joint uncorrelated fit in the range $\tilde{t}\in[10,40]$ for Eq. (21-23)
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.

Figure 7: The ratio of the PCAC quark mass determined from the quarkoinum
(pion and $\eta_{s}$), and also kaon. Two definitions are consistent upto
${\cal O}(a^{2})$ correction.

Through the above fit, we can extract the light and strange quark masses
through the pion and $\eta_{s}$ correlators, respectively:

$\displaystyle\partial_{\mu}(\bar{l}\gamma_{5}\gamma_{\mu}l)=2m^{\rm
PC}_{l}\bar{l}\gamma_{5}l,\
\partial_{\mu}(\bar{s}\gamma_{5}\gamma_{\mu}s)=2m^{\rm
PC}_{s}\bar{s}\gamma_{5}l.$ (25)

Alternatively, the sum of the light and strange quark masses can also be
extracted from the kaon correlator,

$\displaystyle\partial_{\mu}(\bar{s}\gamma_{5}\gamma_{\mu}l)=(m^{\rm
PC}_{l}+m^{\rm PC}_{s})\bar{s}\gamma_{5}l.$ (26)

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_{\pi}\sim 300$ MeV and
$m_{\eta_{s}}\sim 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.

Figure 8: The dimensionless PCAC quark mass $\tilde{m}_{q}^{\rm PC}=m_{q}^{\rm
PC}a$ v.s. the bare quark mass $\tilde{m}_{q}^{\rm b}=m_{q}^{\rm b}a$ at three
lattice spacings. The slope should approach 1 in the continuum limit.

In Fig. 8, we further compare the PCAC quark mass $\tilde{m}_{q}^{\rm PC}$
with original input quark mass parameters $\tilde{m}_{q}^{b}$, at three
lattice spacings with $m_{\pi}\sim 300$ MeV. The figure also includes linear
fits using the following form:

$\displaystyle\tilde{m}^{\rm PC}=k_{m}(\tilde{m}^{b}-\tilde{m}_{\rm crti}),$
(27)

where $\tilde{m}_{\rm crti}$ corresponds to the critical pion mass that makes
the pion mass and $\tilde{m}^{\rm PC}$ vanish. The parameter $k_{m}=1+{\cal
O}(a^{2},\alpha_{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 ${\cal O}(a^{2})$ discretization error and
${\cal O}(\alpha_{s})$ loop effects at finite lattice spacing.

Table 4:  The bare coupling $\alpha_{s}^{\rm b}$, critical quark mass $\tilde{m}_{\rm crti}$ and slope $k_{m}$ at three lattice spacing $a$ and $m_{\pi}\sim$ 300 MeV. $a$ (fm) | 0.105 | 0.077 | 0.052
---|---|---|---
$\alpha_{s}^{\rm b}(a)$ | 0.2397 | 0.2234 | 0.2035
$\tilde{m}_{\rm crti}$ | -0.28560(4) | -0.23545(3) | -0.18885(1)
$k_{m}$ |  0.881 (1) |  0.953 (1) |  1.009(1)

Based on the numerical results listed in Table 4, we observe that
$\tilde{m}_{\rm crti}$ remains negative even after a naive ${\cal O}(a)+{\cal
O}(a^{2})$ extrapolation to the continuum, with a value of -0.0865. Therefore,
it is crucial to include the ${\cal O}\big{(}(\alpha^{b}_{s})^{2}\big{)}$
term, where
$\alpha_{s}^{b}=\frac{g_{0}^{2}}{4\pi}=\frac{10}{4\pi\hat{\beta}u_{0}^{4}}$,
in the continuum extrapolation to ensure that $\tilde{m}_{\rm crti}\
{}_{\overrightarrow{a\rightarrow 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.

####  IV.1.2 lattice spacing determination

Figure 9: The lattice spacings determined by $w_{0}$ with different gauge
action improvement factor $c_{1}$.

Fig. 9 illustrates the lattice spacing determined by the gradient flow Lüscher
(2010) with $w_{0}$ Borsanyi  _et al._ (2012) at three bare couplings
$\hat{\beta}$ and $m_{\pi}\sim 300$ MeV, using different gauge action
improvement coefficient $c_{1}$ in the flow: $c_{1}=-0.331$ (Iwasaki action),
$-0.2$ and $-1/(12u_{0}^{2})$ 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.

Table 5: $a_{w_{0}}$ on different ensembles and the fitted values through the functional form in Eq. (28).  | $a_{w_{0}}$ | $a^{\rm fit}_{w_{0}}$ | fit parameters
---|---|---|---
C24P34 | $0.11198(17)$ | $0.11202(16)$ | $a(6.20)$ | $0.10530(18)$
C24P29 | $0.10811(15)$ | $0.10834(09)$ | $a(6.41)$ | $0.07746(18)$
C32P29 | $0.10858(11)$ | $0.10836(06)$ | $a(6.72)$ | $0.05187(26)$
C32P23 | $0.10637(13)$ | $0.10637(06)$ | $c_{l}$ | $0.507(18)$
C48P23 | $0.10631(06)$ | $0.10635(05)$ | $c_{s}$ | $0.110(23)$
C48P14 | $0.10583(07)$ | $0.10582(06)$ | $c_{L}$ | $0.001(09)$
F32P30 | $0.08044(10)$ | $0.08044(08)$ | $c_{u_{0}}$ | $-408(125)$
F48P30 | $0.08017(05)$ | $0.08017(05)$ | $c_{v_{0}}$ |  $379(110)$
F32P21 | $0.07758(16)$ | $0.07757(15)$ |  |
F48P21 | $0.07810(07)$ | $0.07809(06)$ |  |
H48P32 | $0.05430(11)$ | $0.05430(11)$ |  |

Using the FLAG average value of the gradient flow scale, $w_{0}=0.1736(9)$ fm
Aoki  _et al._ (2022), the lattice spacing $a_{w_{0}}$ for each ensemble is
determined using the Symanzik flow and summarized in Table 5. It is evident
that $a_{w_{0}}$ primarily depends on the gauge coupling $\hat{\beta}$, while
the quark mass also exhibits significant effects, as indicated by the precise
$a_{w_{0}}$ results with a statistical uncertainty of approximately 0.1%.
Empirically, $a_{w_{0}}$ can be described by the following parameterization,
yielding a $\chi^{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):

$\displaystyle
a_{w_{0}}(\hat{\beta},\tilde{m}_{\pi},\tilde{m}_{\eta_{s}},\delta u_{0},\delta
v_{0})=a(\hat{\beta})\big{[}1+c_{l}(\frac{\tilde{m}_{\pi}^{2}}{a(\hat{\beta})^{2}}-m_{\pi,{\rm
phys}}^{2})$
$\displaystyle\quad\quad\quad+c_{s}(\frac{\tilde{m}_{\eta_{s}}^{2}}{a(\hat{\beta})^{2}}-m_{\eta_{s},{\rm
phys}}^{2})+c_{L}e^{-\tilde{m}_{\pi}\tilde{L}}$
$\displaystyle\quad\quad\quad+c_{u_{0}}(u_{0}-u_{0}^{I})+c_{v_{0}}(v_{0}-v_{0}^{I})\big{]}\;,$
(28)

Here, $m_{\pi,{\rm phys}}=134.98$ MeV Zyla  _et al._ (2020) represents the
physical pion mass without QED correction, and $m_{\eta_{s},{\rm
phys}}=689.63(18)$ MeV Borsanyi  _et al._ (2021) 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(\beta)$
with physical quark masses at each $\beta$, along with the other fitting
parameters $c_{l,s,L,u_{0},v_{0}}$ are also provided in Table 5.

Table 6:  The tadpole improvement factors $u^{\rm I}_{0}$ and $v^{\rm I}_{0}$
used in the gauge and fermion actions, $u_{0}$ and $v_{0}$ measured from the
realistic configurations generated using $u^{\rm I}_{0}$ and $v^{\rm I}_{0}$,
and their averages $\bar{u}_{0}$ and $\bar{v}_{0}$.

| C24P34 | C24P29 | C32P29 | C32P23 | C48P23 | C48P14 | F32P30 | F48P30 | F32P21 | F48P21 | H48P32
---|---|---|---|---|---|---|---|---|---|---|---
$u^{\rm I}_{0}$ | 0.855453 | 0.855453 | 0.855453 | 0.855520 | 0.855520 | 0.855548 | 0.863437 | 0.863473 | 0.863488 | 0.863499 | 0.873378
$u_{0}$ | 0.855255(7) | 0.855439(2) | 0.855429(2) | 0.855528(2) | 0.855523(1) | 0.855530(2) | 0.863460(1) | 0.863459(1) | 0.863519(2) | 0.863515(1) | 0.873372(1)
$\bar{u}_{0}$ | 0.855354(4) | 0.855446(1) | 0.855441(1) | 0.855524(1) | 0.855522(1) | 0.855539(1) | 0.863449(1) | 0.863466(1) | 0.863504(1) | 0.863507(1) | 0.873375(1)
$v^{\rm I}_{0}$ | 0.951479 | 0.951479 | 0.951479 | 0.951545 | 0.951545 | 0.951570 | 0.956942 | 0.956984 | 0.957017 | 0.957006 | 0.963137
$v_{0}$ | 0.951275(6) | 0.951461(2) | 0.951452(2) | 0.951550(2) | 0.951547(1) | 0.951554(2) | 0.956968(1) | 0.956967(1) | 0.957024(1) | 0.957019(1) | 0.963134(1)
$\bar{v}_{0}$ | 0.951377(3) | 0.951470(1) | 0.951466(1) | 0.951548(1) | 0.951546(1) | 0.951562(1) | 0.956955(1) | 0.956976(1) | 0.957021(1) | 0.957013(1) | 0.963136(1)

We observe that the volume dependence is consistent with zero, but there is
still a non-vanishing mismatch effect ($u_{0}^{I}\neq u_{0}$ and
$v_{0}^{I}\neq v_{0}$) in the tadpole improvement factors. It is also
interesting to note that the dependence of $a_{w_{0}}$ 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.

####  IV.1.3 Mismatch effect of the tadpole improvement factors

As shown in Table 5, the mismatch effect of the tadpole improvement factors
can have a non-zero impact on the determination of $a_{w_{0}}$ 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^{\rm I}_{0}$ and $v^{\rm 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 $a{w_{0}}$ is mainly canceled out due to the opposite
signs of $c_{u_{0}}$ and $c_{v_{0}}$.

Figure 10: The combined quark mass $2\tilde{m}_{l}^{\rm PC}+0.55(1)\tilde{m}^{\rm PC}_{s}$ dependence of $\bar{u}_{0}=\frac{u_{0}+u_{0}^{I}}{2}$ (upper panel) and $\bar{v}_{0}=\frac{v_{0}+v_{0}^{I}}{2}$ (lower panel), at $10/g^{2}$=6.2.  Table 7: The joint fit results with the same $\tilde{m}_{l}^{b}$ and either $c_{\rm sw}=1/(v^{I}_{0})^{3}$ or $1/\bar{v}_{0}^{3}$, on the test ensemble C24P34 which has the largest mismatch effect.  | $\tilde{m}_{\rm PS}$ | $\tilde{m}_{l}^{\rm PC}$ | $\tilde{f}_{\rm PS}$ | $Z_{wp}$
---|---|---|---|---
$1/(v^{I}_{0})^{3}$=1.1609 | 0.1832(12) | 0.01191(13) | 0.0768(10) | 10.50(30)
$1/(\bar{v}_{0})^{3}$=1.1613 | 0.1822(12) | 0.01178(13) | 0.0768(11) | 10.50(31)
difference | 0.0011(01) | 0.00013(01) | 0.0000(00) |  0.00(01)

After averaging $u_{0}$ and $u^{\rm I}_{0}$ to obtain an estimate of the self-
consistent tadpole improvement factor $\bar{u}_{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.10, the value of $\bar{u}_{0}$ exhibits a linear
behavior with respect to the combined quark mass $2\tilde{m}_{l}^{\rm
PC}+0.55(1)\tilde{m}^{\rm PC}_{s}\propto
0.5\tilde{m}_{\pi}^{2}+0.14\tilde{m}_{\eta_{s}}^{2}$. This suggests that the
dependence on the strange quark mass exhibits a similar suppression as
observed in $a_{w_{0}}$, as shown in Table 5. 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^{PP}_{2,wp}(t;\tilde{m}_{l}^{b},c_{\rm sw})$ using the same
$\tilde{m}_{l}^{b}=-0.2770$, but with either $c_{\rm sw}=1/(v^{I}_{0})^{3}$ or
$1/\bar{v}_{0}^{3}$. The joint fit parameters defined in Eq. (21-23) are
listed in Table7. We find that the matrix elements $\tilde{f}_{\rm PS}$ and
$Z{wp}$ remain unchanged within 0.05% uncertainties. However, $\tilde{m}_{\rm
PS}$ experiences a 0.6% change, resulting in a 1.1% shift in the quark mass.
This can be understood by noting that increasing $c_{\rm sw}$ leads to an
increase in $\tilde{m}_{\rm crti}$, thereby making the multiplicative
renormalizable quark mass $\tilde{m}^{\rm PC}_{l}$ larger for the same
$\tilde{m}_{l}^{b}$. 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.

###  IV.2 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 ${\overline{\text{MS}}}$ scheme, can only be obtained
through regularization-independent (RI) schemes such as RI/MOM Martinelli  _et
al._ (1995) or SMOM Aoki  _et al._ (2008); Sturm  _et al._ (2009). 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 He  _et al._ (2022). 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.

Table 8: Summary of uncertainties of RCs in percentage on the C32P29 ensemble
through the intermediate RI/MOM scheme.

Source | $Z^{\rm\bar{MS}}_{q}/Z_{V}$ | $Z^{\rm\bar{MS}}_{S}/Z_{V}$ | $Z^{\rm\bar{MS}}_{P}/Z_{V}$ | $Z^{\rm\bar{MS}}_{T}/Z_{V}$
---|---|---|---|---
$\rm{Statistical\ error}$ | $0.13\%$ | $0.37\%$ | $0.92\%$ | $0.06\%$
$\rm{Lattice\ spacing}$ | $<0.01\%$ | $<0.01\%$ | $<0.01\%$ | $<0.01\%$
$\rm{Finite\ volume\ effect}$ | $0.03\%$ | $0.24\%$ | $0.62\%$ | $0.02\%$
$\rm{Fit\ range\ of}\ a^{2}\mu^{2}$ | $0.11\%$ | $0.47\%$ | $1.21\%$ | $0.31\%$
$\Lambda_{QCD}^{\overline{\rm{MS}}}$ | $<0.01\%$ | $1.64\%$ | $1.65\%$ | $<0.01\%$
$\rm{Truncation\ in\ matching}$ | $0.24\%$ | $3.45\%$ | $3.45\%$ | $0.24\%$
$\rm{Perturbative\ running}$ | $0.02\%$ | $0.07\%$ | $0.07\%$ | $0.02\%$
$\rm{Total}$ | $0.32\%$ | $3.74\%$ | $4.26\%$ | $0.31\%$

Following a similar strategy employed by Ref. He  _et al._ (2022), 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}\mu^{2}$ fit
range. The second source of uncertainty arises from perturbative matching,
including uncertainties associated with $\Lambda_{QCD}$, truncation in the
perturbative matching, and perturbative scale running. These uncertainties are
fully correlated across different ensembles. As shown in Table 8 for the
C32P29 ensemble, the truncation error is the largest source of uncertainty for
$Z_{S,P}$.

In Ref. He  _et al._ (2022), the truncation error resulting from the 3-loop
perturbative matching between the RI/MOM and ${\overline{\text{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:

$\displaystyle C_{S,n_{f}=3,{\rm
3-loop}}^{{\overline{\text{MS}}},\text{RI/MOM}}$
$\displaystyle=1+0.4244\alpha_{s}+1.007\alpha_{s}^{2}+2.722\alpha_{s}^{3},$
(29) $\displaystyle C_{S,n_{f}=3,{\rm
fake4-loop}}^{{\overline{\text{MS}}},\text{RI/MOM}}$
$\displaystyle=1+0.4244\alpha_{s}+1.007\alpha_{s}^{2}+2.722\alpha_{s}^{3}$
$\displaystyle\quad\quad+7.358\alpha_{s}^{4}.$ (30)

The coefficient of fake $\alpha_{s}^{4}$ term $2.722^{2}/1.007=7.358$ is the
same as that provided by the Pad$\acute{\rm e}$ approximation,

$\displaystyle C_{S,n_{f}=3,{\rm
Pade}_{3}}^{{\overline{\text{MS}}},\text{RI/MOM}}$
$\displaystyle=\frac{1-2.279\alpha_{s}-0.1402\alpha^{2}_{s}}{1-2.703\alpha_{s}}$
$\displaystyle=1+0.4244\alpha_{s}+1.007\alpha_{s}^{2}+2.722\alpha_{s}^{3}$
$\displaystyle\quad\quad+7.358\alpha_{s}^{4}+\mathcal{O}(\alpha_{s}^{5}).$
(31)

Based on the precise 4-loop calculation, the precise 4-loop matching is Gracey
(2023),

$\displaystyle C_{S,n_{f}=3,{\rm
4-loop}}^{{\overline{\text{MS}}},\text{RI/MOM}}$
$\displaystyle=1+0.4244\alpha_{s}+1.007\alpha_{s}^{2}+2.722\alpha_{s}^{3}$
$\displaystyle\quad\quad+8.263\alpha_{s}^{4},$ (32)

which is larger than $C_{S,n_{f}=3,{\rm
fake4-loop}}^{{\overline{\text{MS}}},\text{RI/MOM}}$ but smaller than
$C_{S,n_{f}=3,{\rm Pade}_{3}}^{{\overline{\text{MS}}},\text{RI/MOM}}$.
Considering that the existence of the 5-loop correction, we anticipate that
the Pad$\acute{\rm e}$ approximation will provide a more accurate estimate for
$C_{S,n_{f}=3}$.

Thus we shall use the 4-loop Pad$\acute{\rm e}$ approximation,

$\displaystyle C_{S,n_{f}=3,{\rm
Pade}_{4}}^{{\overline{\text{MS}}},\text{RI/MOM}}$
$\displaystyle=\frac{1-2.611\alpha_{s}-0.2813\alpha^{2}_{s}-0.3349\alpha^{3}_{s}}{1-3.036\alpha_{s}}$
$\displaystyle=1+0.4244\alpha_{s}+1.007\alpha_{s}^{2}+2.722\alpha_{s}^{3}$
$\displaystyle\quad\quad+8.263\alpha_{s}^{4}+25.084\alpha_{s}^{5}+\mathcal{O}(\alpha_{s}^{6}),$
(33)

as the central value of the matching coefficent, and taking the difference
between $C_{S,n_{f}=3,{\rm Pade}_{4}}$ and $C_{S,n_{f}=3,{\rm Pade}_{3}}$ 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.2.1.
Section IV.2.1 will then present the investigation of the chiral symmetry
breaking effect between $Z_{V}$ and $Z_{A}$. Subsequently, in Section IV.2.4,
we will examine the analogous investigation of $Z_{P}$ and $Z_{S}$, with
additional discussions on $Z_{m}$ in Section IV.2.3. Finally, for
completeness, the case of the tensor current will be discussed in Section
IV.2.5.

####  IV.2.1 Vector normalization and $Z_{q}$

Figure 11: $Z_{V}$ from the vector current conservation of pion correlator at
three lattice spacing.

Unlike its continuum counterpart, the local vector current under lattice
regularization is subject to both the ${\cal O}(a^{2})$ discretization error
and ${\cal O}(\alpha_{s})$ loop corrections Zhang  _et al._ (2021), and then
requires additional normalization. The normalization constant $Z_{V}$ can be
determined from the vector current conservation condition,

$\displaystyle Z_{V}\frac{\langle H|V_{4}|H\rangle}{\langle H|H\rangle}=1,$
(34)

where $V_{\nu}=\bar{\psi}\gamma_{\nu}\psi$ and $H$ represents an arbitrary
hadronic state. Thus extracting $Z_{V}$ from the pion correlator in the rest
frame will be the cheapest choice Zhang  _et al._ (2021),

$\displaystyle Z_{V}=$
$\displaystyle\sum_{\vec{x}}\langle\mathrm{Tr}[S^{\dagger}_{C}(\vec{x},\tilde{T}/2;0;m_{q})S_{C}(\vec{x},\tilde{T}/2;0;m_{q})]\rangle/$
$\displaystyle\big{\\{}2\sum_{\vec{x},\vec{z}}\langle\mathrm{Tr}[S^{\dagger}_{C}(\vec{x},\tilde{T}/2;0;m_{q})S(\vec{x},\tilde{T}/2;\vec{z},t;m_{q})$
$\displaystyle\quad\gamma_{t}S_{C}(\vec{z},t;y_{0};m_{q})]\rangle\big{\\}}|_{\tilde{T}\gg
t\gg 0},$ (35)

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\ll t\ll\tilde{T}/2$, are shown in Fig. 11. As
demonstrated inZhang  _et al._ (2021), $Z_{V}$ is affected by $\alpha_{s}$
corrections and cannot be accurately extrapolated to 1 using a simple naive
$a^{2n}$ continuum extrapolation.

Figure 12: Comparison of $Z_{q}^{\rm RI^{\prime}}$ (data points) and
$Z_{q}^{\rm RI^{\prime},ver}$ (bands) at three lattice spacing with different
scale $\mu$. The discrepancy observed in the figure is primarily attributed to
the discretization error.

The quark field RC $Z_{q}^{\rm RI^{\prime}}$ in the regularization independent
(RI) scheme can be accessed through either its definition

$\displaystyle Z_{q}^{\rm RI^{\prime}}(\mu)$
$\displaystyle=\mathop{\textrm{lim}}\limits_{m\rightarrow
0}\frac{-i}{12p^{2}}\textrm{Tr}\Big{[}S^{-1}(p)\not{p}\Big{]}_{p^{2}=\mu^{2}},$
(36)

with $S(p)=\sum_{x}e^{-ip\cdot x}\langle\psi(x)\bar{\psi}(0)\rangle$, or the
vertex correction of the vector current which is equivalent in the continuum,

$\displaystyle Z_{q}^{\rm RI^{\prime},ver}(\mu)$
$\displaystyle=\mathop{\text{lim}}\limits_{m\rightarrow
0}\frac{Z_{V}}{36}\text{Tr}[\Lambda^{\mu}_{V}(p,p)(\gamma_{\mu}-\frac{\not{p}p_{\mu}}{p^{2}})]_{p^{2}=\mu^{2}},$
(37)

where

$\displaystyle\Lambda_{\mathcal{O}}(p_{1},p_{2})$
$\displaystyle=S^{-1}(p_{1})G_{\mathcal{O}}(p_{1},p_{2})S^{-1}(p_{2}),$
$\displaystyle G_{\mathcal{O}}(p_{1},p_{2})$
$\displaystyle=\sum_{x,y}e^{-i(p_{1}\cdot x-p_{2}\cdot
y)}\langle\psi(x)\mathcal{O}(0)\bar{\psi}(y)\rangle.$ (38)

These two definitions are equivalent in the continuum, but they are subject to
different discretization errors. Fig.12 illustrates the comparison between
$Z_{q}^{\rm RI^{\prime}}$ and $Z_{q}^{\rm RI^{\prime},ver}$ at three lattice
spacings. The two definitions become closer either at smaller $\mu$ or smaller
$a$, with $Z_{q}^{\rm RI^{\prime},ver}$ exhibiting a significantly smaller
$a^{2}p^{2}$ error. This situation is reminiscent of the overlap fermions case
He  _et al._ (2022), where the $a^{2}\mu^{2}$ errors have opposite signs.

Given that $Z_{q}^{\rm RI^{\prime},ver}$ 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
Martinelli  _et al._ (1995),

$\displaystyle Z_{q}^{\omega}(\mu^{2})$
$\displaystyle=\mathop{\text{lim}}\limits_{m\rightarrow
0}\frac{Z_{V}}{48}\text{Tr}[\Lambda^{\mu}_{V}(p_{1},p_{2})\gamma_{\mu}]_{p_{1}^{2}=p_{2}^{2}=\mu^{2},(p_{1}-p_{2})^{2}=\omega\mu^{2}},$
(39)

The discrepancy between $Z_{q}^{\omega}$ with different $\omega$ and
$Z_{q}^{\rm RI^{\prime}}$ should be eliminated through their respective
perturbative matchings, which can be calculated in the continuum. For the SMOM
condition with $\omega=1$, the convergence of the perturbative matching is
poorer compared to that of $Z_{q}$ using the q-projection definition Sturm
_et al._ (2009), but the discretization error will be smaller.

####  IV.2.2 Chiral symmetry breaking between $Z_{A}$ and $Z_{V}$

Figure 13: The ratio $Z_{A}/Z_{V}$ 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

$\displaystyle\frac{Z_{A}}{Z_{V}}(\mu^{2},\omega)$
$\displaystyle=\frac{\text{Tr}[\Lambda^{\mu,{\overline{\text{MS}}}}_{A}(p_{1},p_{2})\gamma_{5}\gamma_{\mu}]}{\text{Tr}[\Lambda^{\mu,{\overline{\text{MS}}}}_{V}(p_{1},p_{2})\gamma_{\mu}]}|_{p_{1}^{2}=p_{2}^{2}=\mu^{2},(p_{1}-p_{2})^{2}=\omega\mu^{2}}$
(40)

can deviate from unity due to the additive chiral symmetry breaking present in
the action. As depicted in Fig. 13 for the RI/MOM scheme (data points,
$\omega=0$) and SMOM scheme (bands, $\omega=1$) at three different lattice
spacings with $m_{\pi}\sim$ 300 MeV, we observe that the breaking diminishes
as $a^{2}\mu^{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 $\alpha_{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}\mu^{2}$.

Figure 14: Renormalized pion decay constant $f_{\pi}$ with $m_{\pi}=317$ MeV
at three lattice spacing, using $Z_{A}=Z_{V}\frac{Z_{A}}{Z_{V}}$ through
either RI/MOM or SMOM scheme. The extrapolated values deviate by 1.3(8)%.

If we fit the $Z_{A}/Z_{V}$ data with a polynomial form in the range of
$9\mathrm{GeV}^{2}\leq p^{2}\leq 15/a^{2}$ and extrapolate it to
$a^{2}p^{2}=0$, the obtained result will be sensitive to the choice of
$\omega$. Therefore, we extrapolate the $f_{\pi}$ at two coarser lattice
spacings to the unitary pion mass at the finest lattice spacing, which is
$m_{\pi}=317$ MeV. The renormalized $f^{R}_{\pi}=Z_{A}f_{\pi}$ is shown in
Fig. 14, 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.

Figure 15: Valence pion mass square $m_{\pi,vv}^{2}$ dependence of
$\hat{Z}^{\omega=0}_{P}$ (RI/MOM, lighter points) and $\hat{Z}^{\omega=1}_{P}$
(SMOM, darker points) with $a^{2}\mu^{2}=4$, at three lattice spacing and
$m_{\pi,ss}\sim 300$ MeV

####  IV.2.3 Pesudoscalar current renormalization $Z_{P}$ and induced quark
mass

The renormalization of $m_{q}^{\rm PC}$ requires both $Z_{A}$ and $Z_{P}$. By
defining $\hat{Z}^{\omega}_{P}$ under either the RI/MOM ($\omega=0$) or SMOM
($\omega=1$) scheme,

$\displaystyle\hat{Z}^{\omega}_{P}(\mu)$
$\displaystyle=\frac{Z^{\omega}_{q}}{\frac{1}{12}\text{Tr}[\Lambda_{P}(p_{1},p_{2})\gamma_{5}]}|_{p_{1}^{2}=p_{2}^{2}=\mu^{2},(p_{1}-p_{2})^{2}=\omega\mu^{2}}$
$\displaystyle=Z_{V}\frac{\frac{1}{48}\text{Tr}[\Lambda^{\mu}_{V}(p_{1},p_{2})\gamma_{\mu}]}{\frac{1}{12}\text{Tr}[\Lambda_{P}(p_{1},p_{2})\gamma_{5}]}|_{p_{1}^{2}=p_{2}^{2}=\mu^{2},(p_{1}-p_{2})^{2}=\omega\mu^{2}},$
(41)

$Z_{P}$ in the chiral limit can be extracted from the following
parameterization,

$\displaystyle\hat{Z}^{\omega}_{P}(\mu;m_{q}^{\rm
PC})=\left[(\frac{A_{P}^{\omega}(\mu)}{m_{q}^{\rm
PC}}+(Z^{\omega}_{P}(\mu))^{-1}+C_{P}^{\omega}(\mu)m_{q}^{\rm
PC}\right]^{-1},$ (42)

where $m_{q}^{\rm PC}$ can be replaced by $m_{\pi}^{2}$ based on the GMOR
relation. As shown in Fig.15, $A_{P}^{\omega=1}$ is negligible in the fit, but
$A_{P}^{\omega=0}$ will be non-zero due to the mass pole of the Goldstone
meson and is related to the dynamical quark mass under the Landau gaugeChang
_et al._ (2021).

Figure 16: Values of $m^{\rm RI}_{q}$ (data points) and $\hat{m}^{\rm RI}_{q}$
(bands) at three lattice spacing with $m_{\pi}\sim 300$ MeV at different
RI/MOM scale $\mu$. $\hat{m}^{\rm RI}_{q}$ (black crosses) using the overlap
fermion are also shown on the figure for comparison.

If one define the RI quark mass at given scale as Chang  _et al._ (2021),

$\displaystyle m^{\rm
RI}_{q}=\frac{\frac{1}{12}\text{Tr}[S^{-1}(p)]|_{p^{2}=\mu^{2}}}{Z^{\rm
RI}_{q}(\mu)},$ (43)

then for the overlap fermion action with exact chiral symmetry, we can further
define $Z^{\rm RI}_{m}=m^{\rm RI}_{q}/m_{q}^{b}$ and have the relation
$Z_{m}\hat{Z}_{p}=1$ holding for arbitrary $m_{q}$ and $\mu$. Thus, the
following definition of the RI quark mass:

$\displaystyle\hat{m}^{\rm RI}_{q}=\frac{Z_{A}m_{q}^{\rm PC}}{\hat{Z}^{\rm
MOM}_{P}(\mu)},$ (44)

ensures that $m^{\rm RI}_{q}=\hat{m}^{\rm RI}_{q}$ given the relation
$Z_{A}m_{q}^{\rm PC}=m_{q}^{b}$ for the overlap fermion. However, the case of
clover fermions can be quite different due to its additive chiral symmetry
breaking.

In Fig.16, we present both $m^{\rm RI}_{q}$ and $\hat{m}^{\rm RI}q$ at three
lattice spacings with $m_{\pi}\sim 300$ MeV, and compare them with the results
obtained from the overlap fermion at $a=0.111$ fm. We observe that
$\hat{m}^{\rm RI}_{q}$ seems to be insensitive to the lattice spacing and is
close to the overlap result $\hat{m}^{\rm RI}_{q}=m^{\rm RI}_{q}$Chang  _et
al._ (2021); but $m^{\rm RI}_{q}$ exhibits a significant $a^{2}\mu^{2}$ error
and some unknown UV effects, resulting in large values at large $\mu$.

Figure 17: $a^{2}\mu^{2}$ dependence of $Z^{{\overline{\text{MS}}}\rm{(2\
GeV)}}_{P}$ through the RI/MOM ($\omega=0$, upper panel) and SMOM ($\omega=1$,
lower panel) schemes, at three lattice spacing and $m_{\pi,ss}\sim 300$ MeV.
Figure 18: Renormalized quark mass with $m_{\pi}$=317 MeV at three lattice
spacing, using $Z^{{\overline{\text{MS}}}\rm{(2\ GeV)}}_{P}$ through either
RI/MOM or SMOM scheme. The extrapolated values deviate by 3.1(1.5)%.

After $Z^{\omega}_{P}$ is extracted, it can be further converted to the
${\overline{\text{MS}}}$ scheme at given scale likes $\mu_{0}$= 2 GeV through
the perturbative matching, and shall be independent of both $\mu$ and
$\omega$. In Fig. 17, $Z^{{\overline{\text{MS}}}\rm{(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. 18, we show the renormalized quark mass at
${\overline{\text{MS}}}\rm{(2\ GeV)}$ with $m_{\pi}$=317 MeV at three lattice
spacing, using $Z^{{\overline{\text{MS}}}\rm{(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.

Figure 19: Valence pion mass square $m_{\pi,vv}^{2}$ dependence of
$\hat{Z}^{\omega=0}_{S}$ (RI/MOM, lighter points) and $\hat{Z}^{\omega=1}_{P}$
(SMOM, darker points) with $a^{2}\mu^{2}=4$, at three lattice spacing and
$m_{\pi,ss}\sim 300$ MeV.

####  IV.2.4 Chiral symmetry breaking between $Z_{P}$ and $Z_{S}$

The scalar current RC can be defined similarly with a slightly different
parameterization,

$\displaystyle\hat{Z}^{\omega}_{S}(\mu)=Z_{V}\frac{\frac{1}{48}\text{Tr}[\Lambda^{\mu}_{V}(p_{1},p_{2})\gamma_{\mu}]}{\frac{1}{12}\text{Tr}[\Lambda_{S}(p_{1},p_{2})]}|_{p_{1}^{2}=p_{2}^{2}=\mu^{2},(p_{1}-p_{2})^{2}=\omega\mu^{2}},$
$\displaystyle\hat{Z}^{\omega}_{S}(\mu;m_{q}^{\rm
PC})=\frac{A_{S}^{\omega}(\mu)}{m_{q}^{\rm
PC}}+(Z^{\omega}_{S}(\mu))+C_{S}^{\omega}(\mu)m_{q}^{\rm PC},$ (45)

where $S=\bar{\psi}\psi$, and the $A_{S}^{\omega}$ term can be dropped for
$\omega=1$. However, as shown in Fig. 19, the effect of the $A_{S}^{\omega}$
term is consistently negligible regardless of the value of $\omega$.

Figure 20: $a^{2}\mu^{2}$ dependence of $Z^{{\overline{\text{MS}}}\rm{(2\
GeV)}}_{S}$ through the RI/MOM ($\omega=0$, upper panel) and SMOM ($\omega=1$,
lower panel) schemes, at three lattice spacing and $m_{\pi,ss}\sim 300$ MeV.
Figure 21:  Renormalized scalar matrix element $g_{S,\pi}$ with $m_{\pi}$=317
MeV at three lattice spacing. $g_{S,\pi,{\rm ME}}$ used
$Z^{{\overline{\text{MS}}}\rm{(2\ GeV)}}_{S}$, $g_{S,\pi,{\rm FH}}$ used
$Z^{{\overline{\text{MS}}}\rm{(2\ GeV)}}_{P}$, through either RI/MOM or SMOM
scheme. The extrapolated values deviate by $\sim 7$%.

After converting to the ${\overline{\text{MS}}}$ scheme at a given scale using
the same matching procedure as in the pseudoscalar case, we obtain
$Z^{{\overline{\text{MS}}}\rm{(2\ GeV)}}_{S}$ through either the RI/MOM scheme
(upper panel) or the SMOM scheme (lower panel) with the corresponding
$a^{2}\mu^{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.

Then we calculate the scalar matrix element with the valence quark
contribution only,

$\displaystyle g_{S,\pi,{\rm ME}}={\langle H|S|H\rangle}_{\rm val}/{\langle
H|H\rangle}$ (46)

using similar propagators as for $Z_{V}$. We show the renormalized
$g_{S,\pi,{\rm ME}}^{{\overline{\text{MS}}}\rm{(2\
GeV)}}=Z^{{\overline{\text{MS}}}\rm{(2\ GeV)}}_{S}g_{S,\pi}$ obtained with
either the RI/MOM or SMOM scheme in Fig. 2, for $m_{\pi}=317$ MeV at three
different lattice spacings. It can be observed that the RI/MOM scheme exhibits
a smaller discretization error compared to the SMOM scheme, and the continuum
extrapolated values differ from each other by approximately 6.6(2.4)%.

Using the Feynman-Hellman theorem, one can also extract $g_{S,\pi}$ from the
quark mass dependence of $m_{\pi}$, as

$\displaystyle g_{S,\pi,{\rm FH}}=\frac{1}{2}\frac{\partial
m_{\pi}(m_{q})}{\partial m_{q}}\simeq\frac{m_{\pi}}{4m_{q}}+{\cal
O}(m_{q},a^{2}),$ (47)

where the factor $\frac{1}{2}$ in front of $\frac{\partial
m_{\pi}(m_{q})}{\partial m_{q}}$ comes from two propagators in the pion
correlator. Using the renormalized quark mass
$m_{q}^{{\overline{\text{MS}}}\rm{(2\ GeV)}}$ extracted from Fig. 18 with
RI/MOM or SMOM scheme, $g_{S,\pi,{\rm FH}}^{{\overline{\text{MS}}}\rm{(2\
GeV)}}$ (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,\pi,{\rm 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,\pi,{\rm
FH}}$ is consistent with the $g_{S,\pi,{\rm ME}}$ using the RI/MOM scheme at
each lattice spacing within two sigma, but have significant difference from
the $g_{S,\pi,{\rm ME}}$ using the SMOM scheme. Thus the deviation between the
$g_{S,\pi,{\rm FH}}$ and $g_{S,\pi,{\rm 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}$.

Figure 22: The ratio $Z_{S}/Z_{P}$ through the RI/MOM scheme at three lattice
spacing as functions of $a^{2}p^{2}$.

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. 22, 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 ${\cal O}(\alpha_{s})$ term is necessary to restore
chiral symmetry in the continuum.

####  IV.2.5 $Z_{T}$ and brief summary on renormalization

Figure 23: Valence pion mass square $m_{\pi,vv}^{2}$ dependence of
$\hat{Z}^{\omega=0}_{S}$ (RI/MOM, lighter points) and $\hat{Z}^{\omega=1}_{T}$
(SMOM, darker points) with $a^{2}\mu^{2}=4$, at three lattice spacing and
$m_{\pi,ss}\sim 300$ MeV.

Figure 24: $a^{2}\mu^{2}$ dependence of $Z^{{\overline{\text{MS}}}\rm{(2\
GeV)}}_{T}$ through the RI/MOM ($\omega=0$, upper panel) and SMOM ($\omega=1$,
lower panel) schemes, at three lattice spacing and $m_{\pi,ss}\sim 300$ MeV.

For completeness, we also show the tensor current renormalization using both
the RI/MOM and SMOM schemes,

$\displaystyle\hat{Z}^{\omega}_{T}(\mu)=Z_{V}\frac{\frac{1}{48}\text{Tr}[\Lambda^{\mu}_{V}(p_{1},p_{2})\gamma_{\mu}]}{\frac{1}{72}\text{Tr}[\Lambda_{T_{\mu\nu}}(p_{1},p_{2})\sigma_{\mu\nu}]}|_{p_{1}^{2}=p_{2}^{2}=\mu^{2},(p_{1}-p_{2})^{2}=\omega\mu^{2}},$
(48)

where $T_{\mu\nu}=\bar{\psi}\sigma_{\mu\nu}\psi$. Then
$Z^{{\overline{\text{MS}}}\rm{(2\ GeV)}}_{T}$ can be obtained through the
corresponding perturbative matching after performing a linear $m_{\pi}^{2}$
extrapolation to the continuum, as shown in Fig. 23. The systematic
uncertainty arising from the SMOM matching is estimated to be 2-3%, which
implies that $Z^{{\overline{\text{MS}}}\rm{(2\ GeV)}}_{T}$ is considered
independent of the intermediate scheme within this uncertainty.

Table 9: Normalization and renormalization constant at
${\overline{\text{MS}}}\rm{(2\ GeV)}$ at three lattice spacing with
$m_{\pi,ss}\sim$ 300 MeV, using RI/MOM ($\omega=0$) or SMOM ($\omega=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.

a(fm) | $\omega$ | $Z_{A}$ | $Z_{S}$ | $Z_{P}$ | $Z_{T}$
---|---|---|---|---|---
0.1053(2) | 0 | 0.8547(13) | 0.957(08)(34) | 0.727(09)(26) | 0.864(01)(05)
| 0.8388(39) | 0.897(20)(44) | 0.672(36)(31) | 0.867(04)(04)
1 | 0.8177(04) | 0.754(02)(09) | 0.698(05)(08) | 0.848(01)(27)
0.0775(2) | 0 | 0.8821(08) | 0.879(05)(18) | 0.701(05)(15) | 0.923(01)(03)
| 0.8792(25) | 0.892(12)(11) | 0.689(22)(09) | 0.919(05)(02)
1 | 0.8533(06) | 0.738(01)(07) | 0.689(02)(06) | 0.906(01)(22)
0.0519(3) | 0 | 0.9011(04) | 0.796(09)(10) | 0.683(05)(08) | 0.978(01)(02)
| 0.8998(07) | 0.822(08)(04) | 0.683(08)(03) | 0.973(01)(01)
1 | 0.8838(02) | 0.720(01)(05) | 0.684(01)(05) | 0.961(00)(19)

Table 10:  Normalization and renormalization constants at
${\overline{\text{MS}}}\rm{(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.

| C24P34 | C24P29 | C32P29 | C32P23 | C48P23 | C48P14 | F32P30 | F48P30 | F32P21 | F48P21 | H48P32
---|---|---|---|---|---|---|---|---|---|---|---
$Z_{V}$ | 0.79676(32) | 0.79814(23) | 0.79810(13) | 0.79957(13) | 0.79954(05) | 0.79957(06) | 0.83548(12) | 0.83511(04) | 0.83579(09) | 0.83567(05) | 0.86855(04)
$Z_{A}$ | 0.85698(89) | 0.85442(85) | 0.8547(13) | 0.85593(87) | 0.85674(96) | 0.85605(94) | 0.88161(64) | 0.88213(77) | 0.88069(60) | 0.88063(98) | 0.90113(36)
$Z_{S}$ | 0.960(05)(33) | 0.953(06)(31) | 0.957(08)(34) | 0.957(06)(33) | 0.966(10)(34) | 0.963(04)(34) | 0.871(09)(22) | 0.879(05)(18) | 0.871(06)(21) | 0.878(04)(17) | 0.796(09)(10)
$Z_{P}$ | 0.732(11)(26) | 0.733(13)(24) | 0.727(09)(26) | 0.733(08)(26) | 0.733(09)(27) | 0.727(11)(26) | 0.699(07)(17) | 0.701(05)(15) | 0.704(05)(17) | 0.701(08)(15) | 0.683(05)(08)
$Z_{T}$ | 0.865(02)(05) | 0.864(01)(05) | 0.864(01)(05) | 0.865(01)(05) | 0.865(01)(05) | 0.866(01)(05) | 0.922(01)(04) | 0.923(01)(03) | 0.920(01)(04) | 0.922(01)(03) | 0.978(01)(02)

Based on the values of $Z_{A,S,P,T}$ obtained in this work and collected in
Table 9, it can be observed that the dependence on the intermediate scheme
(RI/MOM or SMOM) is reduced as the lattice spacing decreases. However, 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_{\pi}$, $m_{q}$, and $g_{S}$, it is
recommended to use the RI/MOM scheme to suppress the discretization error.

In Table 9, we also show the renormalization constants using the RI/MOM scheme
but with an additional $1/p^{2}$ term Hasan  _et al._ (2019). 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 10 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.

###  IV.3 Global fit

Using the lattice spacing we obtain in Sec.II, 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_{\pi^{0}}$ of 134.98
MeV within 1% statistical uncertainty. The charged pion mass $m_{\pi^{\pm}}$
139.57 MeV receives the QED correction 4.53(6) MeV Feng  _et al._ (2022) and
then the subtracted pure QCD $m_{\pi^{\pm}}$ is consistent with that of
$m_{\pi^{0}}$ within the uncertainty.

The corresponding renormalized light quark mass and pion decay constants can
also be determined using the RCs extracted in Sec.IV.2:

$\displaystyle m_{l}^{{\overline{\text{MS}}}\rm{(2\
GeV)}}(a=0.105~{}\mathrm{fm})$ $\displaystyle=3.64(8)(11)~{}\mathrm{MeV},$
$\displaystyle f_{\pi}(a=0.105~{}\mathrm{fm})$
$\displaystyle=121.9(5)~{}\mathrm{MeV},$ (49)

where the second uncertainty of $m_{q}$ comes from that of the renormalization
constant. Based on the continuum extrapolation shown in Fig. 14 with a 317 MeV
pion mass, the pion decay constant can change by approximately 7%, and then
agree with the present PDG value 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_{\pi}L>3.5$, on each of the 11
ensembles. Then we use the following NLO partially quenched $\chi$PT form
Sharpe (1997) 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,

$\displaystyle m^{2}_{\pi,{\rm vv}}$ $\displaystyle=\Lambda_{\chi}^{2}2y_{\rm
v}\big{\\{}1+\frac{2}{N_{f}}[(2y_{\rm v}-y_{\rm s})\mathrm{ln}(2y_{\rm
v})+(y_{\rm v}-y_{\rm s})]$ $\displaystyle\quad+2y_{\rm
v}(2\alpha_{8}-\alpha_{5})+2y_{\rm s}N_{f}(2\alpha_{6}-\alpha_{4})\big{\\}}$
$\displaystyle\quad\times\left[1+c^{\pi}_{L}e^{-m_{\pi}L}+c^{\pi}_{s}(m_{\eta_{s}}^{2}-m_{\eta_{s},{\rm
phys}}^{2})\right]$ $\displaystyle\quad\times(1+d^{\pi}_{a}a^{2}),$ (50)
$\displaystyle F_{\pi,{\rm vv}}$ $\displaystyle=F(1-\frac{N_{f}}{2}(y_{\rm
v}+y_{\rm s})\mathrm{ln}(y_{\rm v}+y_{\rm s})+y_{\rm v}\alpha_{5}+y_{\rm
s}N_{f}\alpha_{4})$
$\displaystyle\quad\times\left[1+d^{\pi}_{L}e^{-m_{\pi}L}+d^{\pi}_{s}(m_{\eta_{s}}^{2}-m_{\eta_{s},{\rm
phys}}^{2})\right]$ $\displaystyle\quad\times(1+d^{\pi}_{a}a^{2})$ (51)

where $N_{f}=2$ for two light flavors and $\Lambda_{\chi}=4\pi F$ is the
intrinsic scale of $\chi$PT with $F$ being the pion decay constant in the
chiral limit. The dimensionless expansion parameters $y_{\rm v/s}=\frac{\Sigma
m^{\rm v/s}_{l}}{F^{2}\Lambda_{\chi}^{2}}$ involve the chiral condensate
$\Sigma$ and quark mass $m_{v/s}$ for the valence and sea quark masses,
respectively. Additionally, $\alpha_{i}$ represents the NLO low energy
constants of $\chi$PT at the intrinsic scale $\Lambda_{\chi}$. These constants
can be converted to the usual scale-dependent NLO coefficients as
$L_{i}(\mu)=128\pi^{2}(\alpha_{i}+c_{i}\mathrm{ln}\frac{\mu}{\Lambda_{\chi}})$
with
$c_{4,5,6,8}=\\{-\frac{1}{2},-\frac{N_{f}}{2},-\frac{2+N_{f}^{2}}{N_{f}^{2}},\frac{N_{f}^{2}-4}{4N_{f}}\\}$,
respectively. 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_{\eta_{s}}=689.63(18)$ MeV from
Ref. Borsanyi  _et al._ (2021).

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. (21-24).

Table 11: Global fits with and without non-perturabtive uncertainty
$\delta_{\rm np}Z_{P}$ and perturabtive uncertainty $\delta_{\rm p}Z_{P}$.

$\delta Z_{P}$ included | N/A | $\delta_{\rm np}$ | $\delta_{\rm p}$ |  $\delta_{\rm np}$+$\delta_{\rm p}$
---|---|---|---|---
$\chi^{2}$/d.o.f. | 1.5
$F({\rm GeV})$ | 0.08615(67) | 0.08615(69) | 0.08615(70) | 0.08614(71)
$\Sigma^{1/3}({\rm GeV})$ | 0.2683(14) | 0.2685(30) | 0.2682(14) | 0.2684(30)
$\alpha_{4}$ | 0.2683(14) | 0.293(95) | 0.294(88) | 0.296(96)
$\alpha_{5}$ | -0.32(17) | -0.32(18) | -0.31(17) | -0.32(18)
$\alpha_{6}$ | 0.019(38) | 0.018(79) | 0.018(38) | 0.020(77)
$\alpha_{8}$ | 0.508(73) | 0.50(14) | 0.510(76) | 0.50(15)
$c^{\pi}_{a}({\rm fm}^{-2})$ | 1.47(65) | 1.4(1.6) | 1.4(2.8) | 1.5(3.1)
$c^{\pi}_{L}$ | 0.54(19) | 0.52(43) | 0.54(19) | 0.54(45)
$c^{\pi}_{s}({\rm GeV}^{2})$ | 0.101(32) | 0.10(13) | 0.101(32) | 0.10(13)
$d^{\pi}_{a}({\rm fm}^{-2})$ | -5.45(43) | -5.45(45) | -5.46(50) | -5.45(52)
$d^{\pi}_{L}$ | -0.69(14) | -0.69(14) | -0.69(14) | -0.69(14)
$d^{\pi}_{s}(\rm{GeV}^{2})$ | 0.224(28) | 0.224(31) | 0.224(28) | 0.224(31)
$m_{l,{\rm phys}}({\rm MeV})$ | 3.573(30) | 3.568(90) | 3.574(38) | 3.571(94)
$f_{\pi,{\rm phys}}({\rm MeV})$ | 130.11(91) | 130.1(1.1) | 130.10(94) | 130.1(1.1)

The lattice spacing uncertainty is sampled with an uniform 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 $\Sigma$, $F$, $\alpha_{4,5,6,8}$ and also
$c_{m/f,a/L/s}$ in Table 11:

(1) Fitting with the statistical uncertainty $\delta\tilde{O}$ from the
dimensionless observable $\tilde{O}$ with $O=m_{\pi},\ m_{q}^{\rm PC}$ and
$f_{\pi}$, and also lattice spacing uncertainty $\delta a$;

(2) Fitting with $\delta\tilde{O}$, $\delta a$ and non-perturbative
uncertainty $\delta_{\rm np}Z_{P}$;

(3) Fitting with $\delta\tilde{O}$, $\delta a$ and perturbative uncertainty
$\delta_{\rm p}Z_{P}$;

(4) Fitting with $\delta\tilde{O}$, $\delta a$, $\delta_{\rm np}Z_{P}$ and
also $\delta_{\rm p}Z_{P}$.

All four cases provide reasonable $\chi^{2}$/d.o.f. values, and similar values
of $F$ which are irrelevant to $Z_{P}$. However, the uncertainty of $\Sigma$
is highly sensitive to $\delta Z_{P}$, as expected. By imposing the conditions
$y_{\rm s}=y_{\rm v}$, $M_{\pi,{\rm vv}}=M_{\pi,{\rm
phys}}=134.98~{}\mathrm{MeV}$, $a\rightarrow 0$, and $L\rightarrow\infty$, we
can extract the light quark mass $m_{l}^{{\overline{\text{MS}}}\rm{(2\ GeV)}}$
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 11.

As shown in the table, the extrapolated quark mass $m_{l}=3.57(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 $\delta 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 $\delta_{\rm 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 $\delta_{\rm 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 $\Sigma$ and $m_{l}$
in Case (2-4) is almost the same and will be cancelled when we consider the
renormalization independent combination $\Sigma m_{l}$.

Currently, 90% of the uncertainty of $m_{l}=3.57(9)$ 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_{\pi}=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_{\pi}$
in the continuum and infinite volume limits to be 130.1(1.1) MeV, which is
consistent with the experimental value of 130.4(2) MeV Zyla  _et al._ (2020).
Based on $F$ extracted above, we predict
$F_{\pi}/F=\sqrt{1/2}f_{\pi}/F={1.071(9)}$.

Our prediction of the light quark mass and the low energy constant are
summarized in Table 2.

Table 12: Summary of our determination on quark masses at ${\overline{\text{MS}}}{\rm(2GeV)}$ and the other quantities , through the intermediate RI/MOM or SMOM schemes, with comparison with FLAG Aoki  _et al._ (2022) and/or PDG Zyla  _et al._ (2020). The difference between two schemes is considered as systematic uncertainty in the combined determination. |  $m_{l}$ (MeV) |  $m_{u}$ (MeV) |  $m_{d}$ (MeV) |  $m_{s}$ (MeV) |  | $m_{s}/m_{l}$ | $m_{u}/m_{d}$
---|---|---|---|---|---|---|---
RI/MOM | 3.571(94) | 2.39(14) | 4.77(16) | 97.8(2.1) |  | 27.35(36) | 0.501(35)
SMOM | 3.441(46) | 2.20(08) | 4.68(08) | 93.8(1.1) |  | 27.24(23) | 0.471(21)
Combined | 3.57(9)(13) | 2.39(14)(19) | 4.77(16)(09) | 97.8(2.1)(4.0) |  | 27.35(36)(13) | 0.501(35)(30)
FLAG/PDG | 3.381(40) | 2.27(09) | 4.67(09) | 92.2(1.0) |  | 27.42(12) | 0.485(19)
|  $\Sigma^{1/3}$(MeV) |  $F$ (MeV) | $F_{\pi}/F$ |  $f_{\pi}$(MeV) |  $f_{K^{\pm}}$(MeV) | $\ell_{3}$ | $\ell_{4}$
RI/MOM | 268.4(3.0) | 86.1(7) | 1.071(9) | 130.1(1.1) | 153.8(0.7) | 2.54(46) | 4.30(07)
SMOM | 268.6(1.8) | 84.7(6) | 1.089(8) | 128.0(0.9) | 151.1(0.6) | 2.59(20) | 4.20(05)
Combined | 268.4(3.0)(0.2) | 86.1(0.7)(1.4) | 1.071(9)(18) | 130.1(1.1)(2.1) | 153.8(0.7)(2.7) | 2.54(46)(05) | 4.30(07)(10)
FLAG/PDG | 272(5) | 86.8(6) | 1.062(7) | 130.4(0.2) | 155.7(0.7) | 3.07(64) | 4.02(45)

Figure 25: The corrected unitary $(m^{\rm uni}_{\pi})^{2}/m_{q}$ and the decay
constant $f^{\rm uni}_{\pi}$ varies with the quark mass, at three lattice
spacing (colored data points and bands) and also continuum (gray band).

As shown in Table 11 and inspired by the previous sections, the discretization
error in quark mass is consistent with zero, while this error in $f_{\pi}$ 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_{\pi}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 Yang  _et al._ (2015b). 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^{\rm
data}_{\pi,{\rm vv}}$ and $f^{\rm data}_{\pi,{\rm vv}}$. We define the
corrected $m^{\rm uni}_{\pi}$ and $f^{\rm uni}_{\pi}$ as follows,

$\displaystyle(m^{\rm uni}_{\pi})^{2}$ $\displaystyle=\frac{(m^{\rm
data}_{\pi,{\rm
vv}})^{2}}{1+c_{m,L}e^{-m_{\pi}L}+c_{m,s}(m_{\eta_{s}}^{2}-m_{\eta_{s},{\rm
phys}}^{2})}$ $\displaystyle\quad\quad-\Lambda_{\chi}^{2}y_{\rm v}(y_{\rm
s}-y_{\rm v})\big{\\{}-\frac{2}{N_{f}}[\mathrm{ln}(2y_{\rm v})+1]$
$\displaystyle\quad\quad\quad+2N_{f}(2\alpha_{6}-\alpha_{4})\big{\\}}(1+c_{m,a}a^{2}),$
$\displaystyle f^{\rm uni}_{\pi}$ $\displaystyle=\frac{f^{\rm data}_{\pi,{\rm
vv}}}{1+c_{f,l}e^{-m_{\pi}L}+c_{f,s}(m_{\eta_{s}}^{2}-m_{\eta_{s},{\rm
phys}}^{2})}$ $\displaystyle\quad\quad-F(y_{\rm s}-y_{\rm
v})\\{-\frac{N_{f}}{2}\mathrm{ln}(y_{\rm v}+y_{\rm s})+N_{f}\alpha_{4}\\}$
$\displaystyle\quad\quad\quad(1+c_{f,a}a^{2}).$ (52)

In Fig. 25, we show the ratio $(m^{\rm
uni}_{\pi})^{2}/m^{{\overline{\text{MS}}}\rm{(2\ GeV)}}_{q}$ (upper panel) and
also $f^{\rm uni}_{\pi}$ (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. We can observe that the
continuum extrapolation actually pushes $f_{\pi}$ to be obviously higher,
while the impact on the $(m_{\pi}^{2}/m_{q})$ ratio and, consequently, $m_{q}$
is much weaker.

The $m_{\pi}$ and $f_{\pi}$ with unitary valence and sea quark masses
$y=y_{\rm v}=y_{\rm s}$ have another widely used parameterization,

$\displaystyle m_{\pi}^{2}$
$\displaystyle=\Lambda_{\chi}^{2}2y[1+y(\ln\frac{2y\Lambda^{2}_{\chi}}{m^{2}_{\pi,{\rm
phys}}}-{\ell}_{3})+\mathcal{O}(y^{2})],$ (53) $\displaystyle F_{\pi}$
$\displaystyle=F[1-2y(\ln\frac{2y\Lambda^{2}_{\chi}}{m^{2}_{\pi,{\rm
phys}}}-\ell_{4})+\mathcal{O}(y^{2})].$ (54)

where $\ell_{3,4}$ is related to $\alpha_{4,5,6,8}$ by

$\displaystyle\ell_{3}=$
$\displaystyle\ln\frac{\Lambda^{2}_{\chi}}{m_{\pi,{\rm
phys}}^{2}}-2[(2\alpha_{8}-\alpha_{5})+2(2\alpha_{6}-\alpha_{4})],$
$\displaystyle\ell_{4}=$
$\displaystyle\ln\frac{\Lambda^{2}_{\chi}}{m_{\pi,{\rm
phys}}^{2}}+\frac{1}{2}(\alpha_{5}+2\alpha_{4}).$ (55)

Our determination of $\ell_{3,4}$ are also collected in Table 12 , consistent
with the current FLAG average but have smaller uncertainties.

In this work, we use the $m_{K^{\pm}}$ and $m_{K^{0}}$ with the constraint
$m^{\rm phys}_{u}+m^{\rm phys}_{d}=2m^{\rm phys}_{l}$, to determine the up,
down and strange quark masses $m_{u,d,s}$. On each ensemble, we calculate the
strange quark propagators with a unitary strange quark mass $m^{\rm
v}_{s}=m^{\rm s}_{s}$, and also 2 partially quenched quark masses $m^{\rm
v}_{s}\sim 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\times 3$ partially quenched Kaon masses and decay constants on all the
ensembles are fitted with the following form proposed in a recent work
Alexandrou  _et al._ (2021),

$\displaystyle m_{K}^{2}(m^{\rm v}_{l},m^{\rm s}_{l},m^{\rm v}_{s},m^{\rm
s}_{s},a)=(b_{s}^{\rm v}m^{\rm v}_{s}+b_{s}^{\rm s}m^{\rm s}_{s}+b_{l}^{\rm
v}m^{\rm v}_{l}+b_{l}^{\rm s}m^{\rm s}_{l})$
$\displaystyle~{}~{}~{}~{}~{}~{}\times\left[1+c^{K}_{l}m^{\rm
v}_{l}+c^{K}_{m}a^{2}+c_{L}^{K}\exp{(-m_{K}L)}\right],$ (56) $\displaystyle
f_{K}(m^{\rm v}_{l},m^{\rm s}_{l},m^{\rm v}_{s},m^{\rm
s}_{s},a)=(d_{f}+d_{s}^{\rm v}m^{\rm v}_{s}+d_{s}^{\rm s}m^{\rm
s}_{s}+d_{l}^{\rm v}m^{\rm v}_{l}+d_{l}^{\rm s}m^{\rm s}_{l})$
$\displaystyle~{}~{}~{}~{}~{}~{}\times\left[1+d^{K}_{a}a^{2}+d_{L}^{K}\exp{(-m_{K}L)}\right].$
(57)

Table 13: Global fit of the Kaon mass and decay constant, with both
$\delta_{\rm np}Z_{P}$ and $\delta_{\rm p}Z_{P}$.

$\chi^{2}$/d.o.f. | $b_{l}^{v}({\rm GeV})$ | $b_{l}^{s}({\rm GeV})$ | $b_{s}^{v}({\rm GeV})$ | $b_{s}^{s}({\rm GeV})$ | $c_{l}^{K}({\rm GeV^{-1}})$ | $c_{a}^{K}({\rm fm^{-2}})$ | $c_{L}^{k}$
---|---|---|---|---|---|---|---
1.1 | 1.86(89) | 0.28(24) | 2.29(11) | 0.11(12) | 2.7(3.4) | 0.6(3.3) | 5.6(9.2)
$\chi^{2}$/d.o.f. | $d_{f}({\rm GeV})$ | $d_{l}^{v}$ | $d_{l}^{s}$ | $d_{s}^{v}$ | $d_{s}^{s}$ | $d_{a}^{K}({\rm fm^{-2}})$ | $d_{L}^{k}$
0.93 | 0.1269(32) | 0.313(31) | 0.528(44) | 0.1442(83) | 0.105(30) | -5.80(72) | -7.3(2.3)

As shown in Table 13, the global fit of the Kaon mass and decay constant
provides reasonable $\chi^{2}$/d.o.f., which suggests a good agreement between
the data points and their fit. Based on the fit of $m_{K}^{2}$, the total
strange quark mass dependence $b_{s}^{\rm v}+b_{s}^{\rm s}=2.40(06)$ is
consistent with the leading order light quark mass dependence $b_{l}^{\rm
v}+b_{l}^{\rm s}=2.14(89)$, and the coefficient of the non-linear quark mass
dependence $c_{l}^{K}=2.7(3.4)$ is also consistent with zero.

Based on the QED correction $\Delta_{\rm QED}m_{K}$ obtained in previous
literature Giusti  _et al._ (2017),

$\displaystyle\Delta_{\rm QED}m_{K^{\pm}}-\Delta_{\rm QED}m_{K^{0}}=2.07(15)\
{\rm MeV},$ $\displaystyle\Delta_{\rm QED}m^{2}_{K^{0}}=0.174(24)\times
10^{-3}\ {\rm GeV}^{2},$ (58)

we can obtain $\Delta_{\rm QED}m_{K^{0}}=0.17(2)$ MeV and $\Delta_{\rm
QED}m_{K^{\pm}}=2.24(15)$ MeV, respectively. Thus, the pure QCD Kaon mass with
physical quark mass will be

$\displaystyle m_{K^{0},{\rm QCD}}=m_{K^{0}}^{\rm phys}-\Delta_{\rm
QED}m_{K^{0}}=497.44(02)\ {\rm MeV},$ $\displaystyle m_{K^{\pm},{\rm
QCD}}=m_{K^{\pm}}^{\rm phys}-\Delta_{\rm QED}m_{K^{\pm}}=491.44(15)\ {\rm
MeV},$

and the physical quark masses $m^{\rm phys}_{u}$, $m^{\rm phys}_{d}$ and
$m^{\rm phys}_{s}$ can be determined using the following three conditions

$\displaystyle m_{K}(m^{\rm phys}_{d},m^{\rm phys}_{l},m^{\rm phys}_{s},m^{\rm
phys}_{s},0)$ $\displaystyle=m_{K^{0},{\rm QCD}},$ $\displaystyle m_{K}(m^{\rm
phys}_{u},m^{\rm phys}_{l},m^{\rm phys}_{s},m^{\rm phys}_{s},0)$
$\displaystyle=m_{K^{\pm},{\rm QCD}},$ $\displaystyle m^{\rm phys}_{u}+m^{\rm
phys}_{d}$ $\displaystyle=2m^{\rm phys}_{l}.$ (60)

Note that we ignored the isospin symmetry breaking effect in the sea quark
masses, since it should be ${\cal O}((m_{d}-m_{u})^{2})$ and negligible given
the current uncertainties.

Figure 26: The corrected kaon mass $m^{\rm cr}_{K}$ and the decay constant
$f^{\rm cr}_{K}$ with the physical light quark mass $m_{l}^{\rm phys}$, varies
with the strange quark mass at three lattice spacing (colored data points and
bands) and also continuum (gray band).

In Fig. 26, we show the corrected kaon mass $m^{\rm cr}_{K}$ and the decay
constant $f^{\rm cr}_{K}$ with the light quark mass $m_{l}$ corrected to its
physical value $m_{l}^{\rm phys}$ using the bootstrap samples of the fitted
$m_{K}$,

$\displaystyle(m^{\rm cr}_{K})^{2}$ $\displaystyle=(m^{\rm
data}_{K})^{2}-m^{2}_{K}(m^{\rm v}_{l},m^{\rm s}_{l},m^{\rm v}_{s},m^{\rm
s}_{s},a,L)$ $\displaystyle\quad+m^{2}_{K}(m^{\rm phys}_{l},m^{\rm
phys}_{l},m^{\rm v}_{s},m^{\rm v}_{s},a,L\rightarrow\infty),$ (61)
$\displaystyle f^{\rm cr}_{K}$ $\displaystyle=f^{\rm data}_{K}-f_{K}(m^{\rm
v}_{l},m^{\rm s}_{l},m^{\rm v}_{s},m^{\rm s}_{s},a,L)$
$\displaystyle\quad+f_{K}(m^{\rm phys}_{l},m^{\rm phys}_{l},m^{\rm
v}_{s},m^{\rm v}_{s},a,L\rightarrow\infty).$ (62)

We can see that $f_{K}$ also exhibits a strong lattice spacing dependence,
similar to the $f_{\pi}$ case, while the Kaon mass is insensitive to the
lattice spacing.

As a consistency check, we also calculate $m_{\eta_{s}}$ on each ensemble and
perform a similar global fit with the following functional form,

$\displaystyle m_{\eta_{s}}^{2}=$ $\displaystyle[h_{l}^{\rm s}m_{l}^{\rm
s}+h_{s}^{\rm v}m_{s}^{\rm v}+h_{s}^{\rm s}m_{s}^{\rm s}+h_{l^{2}}(m_{l}^{\rm
s})^{2}]$ $\displaystyle\quad\quad(1+h_{a}a^{2}+h_{L}e^{m_{\eta_{s}}L}),$ (63)

The $\chi^{2}$/d.o.f. of the fit is 1.4 an all the coefficients except
$h_{s}^{l}$ are consistent with zero within the uncertainty, but dropping more
terms will enlarge the $\chi^{2}$/d.o.f. significantly and then we keep all
the coefficients. Using the physical light and strange quark masses, we
predict $m_{\eta_{s}}=687.4(2.2)$ MeV in the infinite volume and continuum
limit, which is consistent with the BMWc value of 689.63(18) MeV Borsanyi  _et
al._ (2021) we used in the fits for lattice spacing and $f_{\pi}$. Additively,
we extract the dependence of $m_{\eta_{s}}^{2}$ on the strange quark mass,
$\frac{\partial m_{\eta_{s}}^{2}}{\partial m_{s}}=4.9(1)$ GeV, which is
slightly lower than the value of $m_{\pi}^{2}/m_{q}$ shown in Fig. 25.

The physical quark masses $m_{u,d,s}$ and also corresponding $f_{\pi,K}$ using
$m_{l}^{\rm phys}$ and intermediate RI/MOM scheme, are collected in Table 12.
In addition, Table 12 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 $\sim$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. With this systematic uncertainty, all
our determinations are consistent with the present lattice averages Aoki  _et
al._ (2022) and/or PDG Zyla  _et al._ (2020) within 1-2 $\sigma$, and the low
energy constants $\Sigma$ and $\ell_{3,4}$ have smaller uncertainties.