# Determination of the up/down-quark mass within QCD sum rules in the scalar
channel

Fang-Hui Yin1 Wen-Ya Tian1 Liang Tang1<EMAIL_ADDRESS>Zhi-Hui Guo2
<EMAIL_ADDRESS>1 College of Physics and Hebei Key Laboratory of Photophysics
Research and Application, Hebei Normal University, Shijiazhuang 050024, China
2 School of Physics, Southeast University, Nanjing 211189, China

###### Abstract

In this work, we determine up/down-quark mass $m_{q=u/d}$ in the isoscalar
scalar channel from both the Shifman-Vainshtein-Zakharov (SVZ) and the Monte-
Carlo-based QCD sum rules. The relevant spectral function, including the
contributions from the $f_{0}(500)$, $f_{0}(980)$ and $f_{0}(1370)$
resonances, is determined from a sophisticated $U(3)$ chiral study. Via the
traditional SVZ QCD sum rules, we give the prediction to the average light-
quark mass
$m_{q}(2\,\text{GeV})=\frac{1}{2}(m_{u}(2\,\text{GeV})+m_{d}(2\,\text{GeV}))=(3.46^{+0.16}_{-0.22}\pm
0.33)\,\text{MeV}$. Meanwhile, by considering the uncertainties of the input
QCD parameters and the spectral functions of the isoscalar scalar channel, we
obtain $m_{q}(2\,\text{GeV})=(3.44\pm 0.14\pm 0.32)\,\text{MeV}$ from the
Monte-Carlo-based QCD sum rules. Both results are perfectly consistent with
each other, and nicely agree with the Particle Data Group value within the
uncertainties.

##  I Introduction

The up- and down- quark masses are fundamental parameters in Quantum
Chromodynamics (QCD), and they play important roles in the phenomenological
study of the Standard Model. However, they are not physical observables, due
to the confinement nature of QCD. Usually one needs to adopt some QCD-based
non-perturbative methods, such as lattice QCD, QCD sum rules, effective field
theories, etc., to determine the quark masses. A recent review of the
determinations of the light quark masses from lattice QCD can be found in
Ref.Aoki:2019cca .

The QCD sum rules, developed more than forty years ago by Shifman, Vainshtein
and Zakharov (SVZ) Shifman ; P.Col ; Latorre:1985uy , has some peculiar
advantages in the study of hadron phenomenology. Its starting point in
evaluating the properties of the concerned hadron is to construct the current,
which possesses the important information about the concerned hadron, like
quantum numbers, constituent quarks and gluons. With the currents one can then
construct the two-point correlation function, which matches the QCD and
phenomenological sides of sum rules. Within the QCD sum rules, in recent years
there are remarkable progresses on the determination of the light quark
massesBecchi:1980vz ; Bijnens:1994ci ; Jamin:2001zr ; Steele:2001ga ;
Gupta:2003fn ; Bodenstein:2013paa ; Narison:2014vka ; Dominguez:2018azt ;
Cherry:2001cj ; Narison:2002hk ; Maltman:2002sb . Many of these works rely on
the pseudoscalar condensates or the divergences of the axial-vector currents
to determine the light quark masses from QCD sum rules. The golden cases to
determine the $m_{d}+m_{u}$ and $m_{d}-m_{u}$ would be the divergence of the
isovector axial-vector current and the divergence of isovector vector current,
respectively Dominguez:1986aa ; Dominguez:2018azt . Alternatively the
isoscalar scalar channel can be also used to calculate the $m_{d}+m_{u}$
Reinders:1981ww ; Cherry:2001cj ; Yuan:2017foa . In this situation, on the
phenomenological side one needs the isoscalar scalar spectral functions, which
can contain the contributions from the scalar resonances
$f_{0}(500)$($\sigma$), $f_{0}(980)$, $f_{0}(1370)$ and other heavier ones.
Indeed there are serve debates on the roles of the different resonances played
in the scalar QCD sum rules Cherry:2001cj ; Shakin:2001hc ; Elias:1998bq ;
Shi:1999hm , specially on the low lying states $f_{0}(500)$ and $f_{0}(980)$.
In the present work, we attempt to address this issue by adopting the
phenomenological spectral functions calculated within the unitarized $U(3)$
resonance chiral theory Guo:2012ym ; Guo:2012yt ; Guo:2011pa . The spectral
functions are constructed by the unitarized scalar form factors, which are
pure predictions once the scattering amplitudes are determined through the
fits to scattering observables. In the complex energy plane, the isoscalar
scalar resonance poles for the $f_{0}(500)$, $f_{0}(980)$, $f_{0}(1370)$ are
found, whose contributions are coherently included in the spectrum functions,
together with two-meson continuum effects from the chiral effective field
theory. The prescription of the unitarized chiral spectrum functions improves
the Breit-Wigner description of the broad scalar $\sigma$ resonance
Yuan:2017foa . Furthermore, there is no free parameter in our spectrum
functions, since they can be predicted from the scattering observables
Guo:2012ym ; Guo:2012yt ; Guo:2011pa .

Despite the success of QCD sum rules, there remain several subtle issues
concerning their application, such as the way to determine the continuum
threshold parameter $s_{0}$ and the uncertainties of the results. In
Refs.Leinweber:1995fn ; Benmerrouche:1995qa ; Wang:2016sdt ; Yuan:2017foa , a
Monte-Carlo-based QCD sum rule is introduced, which starts from the
theoretical formula of the SVZ QCD sum rule and uses Monte-Carlo based
uncertainty analysis to proceed the study. The Monte-Carlo-based QCD sum rule
has some obvious advantages: (a) the continuum threshold $s_{0}$ can be
automatically determined via the fitting procedure; (b) the effects from the
uncertainties of the nonperturbative OPE parameters and the phenomenological
spectral functions can be more reliably taken into account. Consequently the
Monte-Carlo-based QCD sum rule is considered to be able to quantitatively
improve the uncertainty analysis, comparing with the conventional SVZ sum
rule, which relies on some empirical ways to obtain the platform of the
$s_{0}$ and Borel transformation parameters.

In this work, we aim at using the both the conventional SVZ and Monte-Carlo-
based QCD sum rules in the isoscalar scalar channel to determine the average
of the up- and down-quark masses $(m_{u}+m_{d})/2$, which will be shortly
denoted as $m_{q}$ in later discussions. Comparisons from both methods will be
analyzed in detail. Another novelty of this work is to implement the
sophisticated chiral spectrum functions in the phenomenological side, where
one does not need to manually separate the background and the resonance
contributions. Furthermore the unitarized chiral spectrum functions are
considered to be more suitable for the broad scalar resonances than the Breit-
Wigner forms Yuan:2017foa .

The paper is organized as follows. We elaborate the essential formulas of the
QCD Sum Rules and the chiral spectrum functions in Sec.II. The numerical
analysis and results are given in Sec.III. The last part is left for a short
summary and conclusions.

##  II Theoretical formalism

Within the QCD sum rule approach, the connection between quarks/gluons and
hadrons is through the two-point correlation function, constructed by two
currents, i.e.

$\displaystyle\Pi(q^{2})=i\int d^{4}xe^{iq\cdot
x}\langle\Omega|T\\{j_{S}(x),j_{S}^{\dagger}(0)\\}|\Omega\rangle\;,$ (1)

where $j_{S}$ is the relevant hadronic current, the subscript $S$ represents
the scalar current, and $\Omega$ denotes the physical vacuum. In this work, we
focus on the current in the isoscalar scalar channel

$\displaystyle
j_{S}(x)=m_{q}\frac{1}{\sqrt{2}}\bigg{(}\bar{u}(x)u(x)+\bar{d}(x)d(x)\bigg{)}\;,$
(2)

where $m_{q}=\frac{1}{2}(m_{u}+m_{d})$.

It is stressed that our primary aim in this work is to determine the light-
quark mass by using the tool of the QCD sum rule, instead of scrutinizing the
internal structures of the scalar resonances. For the latter purpose, it can
be crucial to delicately choose the proper currents in the sum-rule analyses.
For example, the $\sigma$ resonance has been extensively studied in the QCD
sum rule approach by using the mixed currents of the $\bar{q}q$ type and the
gluons Narison:1988ts ; Narison:2000dh , and also the various tetraquark
currents Jaffe:1976ig ; Jaffe:2004ph ; Achasov:1980gu ; Latorre:1985uy ;
Narison:2005wc ; Brito:2004tv ; Wang:2005xy ; Lee:2006vk ; Chen:2007zzg ;
Chen:2007xr . However, for the purpose to extract the light-quark mass, it is
adequate to analyze the $j_{S}$ current of Eq. (2) within the framework of QCD
sum rule, as long as one can properly calculate the hadronic contributions to
two-point correlation functions of the current, which will be discussed in
detail in the following.

The correlation function has two equivalent representations: the QCD side and
the phenomenological one. The QCD representation of this function can be
calculated by using the Operator Product Expansion (OPE) method. In addition,
it has been noticed that the instanton effects could be important in the
pseudoscalar and scalar channels Shuryak:1982qx ; Cherry:2001cj ; Elias:1998fs
; Elias:1998bq ; Shi:1999hm ; Yuan:2017foa . Therefore we also include the
instanton contributions in the QCD side, which are calculated within the
instanton liquid model Shuryak:1982qx ; Elias:1998fs . The correlation
function from the QCD side then includes both the OPE and instanton parts

$\displaystyle\Pi^{\textrm{(QCD)}}(q^{2})=\Pi^{\textrm{(OPE)}}(q^{2})+\Pi^{\textrm{(Inst)}}(q^{2})\;.$
(3)

After performing the Borel transformation of Eq.(3), one can obtain
Shuryak:1982qx ; Cherry:2001cj ; Elias:1998fs ; Elias:1998bq ; Shi:1999hm ;
Yuan:2017foa ,

$\begin{split}R^{\textrm{(QCD)}}(\tau,\hat{m}_{q})=&R^{\textrm{(OPE)}}(\tau,\hat{m}_{q})+R^{\textrm{(Inst)}}(\tau,\hat{m}_{q})\\\
=&m_{q}^{2}(1/\sqrt{\tau})\left\\{\frac{3}{8\pi^{2}}\frac{1}{\tau^{2}}\left[1+4.821098\frac{\alpha_{s}(1/\tau)}{\pi}+21.97646\left(\frac{\alpha_{s}(1/\tau)}{\pi}\right)^{2}\right.\right.\\\
&+53.14197\left.\left(\frac{\alpha_{s}(1/\tau)}{\pi}\right)^{3}\right]+\frac{\langle\alpha_{s}G^{2}\rangle}{8\pi}\left(1+\frac{11}{2}\frac{\alpha_{s}(1/\tau)}{\pi}\right)\\\
&+3\langle
m_{q}\bar{q}q\rangle\left(1+\frac{13}{3}\frac{\alpha_{s}(1/\tau)}{\pi}\right)-\frac{176}{27}\pi\kappa\alpha_{s}\langle\bar{q}q\rangle^{2}\left[\frac{\alpha_{s}(1/\tau)}{\alpha_{s}(\mu^{2}_{0})}\right]^{1/9}\tau\\\
&\left.+\frac{3}{8\pi^{2}}\frac{e^{\frac{-\rho^{2}}{2\tau}}\rho^{2}}{\tau^{3}}\left[K_{0}\left(\frac{\rho^{2}}{2\tau}\right)+K_{1}\left(\frac{\rho^{2}}{2\tau}\right)\right]\right\\},\end{split}$
(4)

where the last two terms with the modified Bessel functions $K_{0}$ and
$K_{1}$ stand for the contributions from the instanton, and $\rho$ represents
the instanton size in the instanton liquid model Shuryak:1982qx ; Elias:1998fs
. In Eq. (4), $\tau$ is the Borel parameter, $\kappa$ denotes the vacuum
factorization violation factor that parameterizes the deviation of the four-
quark condensate from a product of two-quark condensates, and $\mu_{0}$ is the
renormalization scale for condensates. The running coupling constant and
running mass have been calculated up to four-loop order in the perturbative
corrections Chetyrkin:1996xa ; Shi:1999hm . The four-loop result with three-
flavor quarks for the running coupling constant at scale $1/\sqrt{\tau}$ takes
the form Chetyrkin:1996xa ; Shi:1999hm

$\displaystyle\alpha_{s}(1/\tau)=\pi\left\\{\frac{4}{9}\frac{1}{L}-\frac{256\ln
L}{729L^{2}}+\frac{1}{L^{3}}\left[\frac{16384\ln^{2}L}{59049}-\frac{16384\ln
L}{59049}+\frac{6794}{59049}\right]\right\\}\;,$ (5)

with $L=\ln\left(\frac{1}{\tau\Lambda_{\text{QCD}}^{2}}\right)$. The QCD scale
is taken as $\Lambda_{\text{QCD}}=0.353\,\text{GeV}$ throughout this work
Narison:2009vy . The explicit scale $1/\sqrt{\tau}$ dependence of the light-
quark mass in Eq.(4) up to four-loop order is Chetyrkin:1996xa ; Shi:1999hm

$\displaystyle\begin{split}m_{q}(1/\sqrt{\tau})&=\hat{m}_{q}\frac{1}{\left(\frac{1}{2}L\right)^{4/9}}\left\\{1+\frac{290}{729}\frac{1}{L}-\frac{256}{729}\frac{\ln
L}{L}+\left(\frac{550435}{1062882}-\frac{80\zeta(3)}{729}\right)\frac{1}{L^{2}}\right.\\\
&-\frac{388736}{531441}\frac{\ln
L}{L^{2}}+\frac{106496}{531441}\frac{\ln^{2}L}{L^{2}}+\left(\frac{2121723161}{2324522934}+\frac{8}{6561}\pi^{4}-\frac{119840}{531441}\zeta(3)\right.\\\
&\left.-\frac{8000}{59049}\zeta(5)\right)\frac{1}{L^{3}}+\left(-\frac{611418176}{387420489}+\frac{112640}{531441}\zeta(3)\right)\frac{\ln
L}{L^{3}}\\\
&\left.+\frac{335011840}{387420489}\frac{\ln^{2}L}{L^{3}}-\frac{149946386}{1162261467}\frac{\ln^{3}L}{L^{3}}\right\\}\,,\end{split}$
(6)

where $\hat{m}_{q}$ is the renormalization-group (RG) invariant light-quark
mass and $\zeta$ stands for the Riemann zeta function.

On the phenomenological side, we need to provide the phenomenological spectral
function which is related to the correlation function through the dispersive
integral. It is typical to write the phenomenological spectral function as

${\rm Im}\Pi^{\textrm{(Phen)}}(s)={\rm Im}\Pi^{\textrm{(Res)}}(s)+{\rm
Im}\Pi^{\textrm{(ESC)}}(s)\,,$ (7)

where in the conventional QCD sum rule the first term in the right-hand side
usually includes the resonance contributions via the Breit-Wigner or
$\delta$-function forms, and the second term corresponds to the excited states
and continuum. The ${\rm Im}\Pi^{\textrm{(ESC)}}(s)$ part incorporates the
excited states and higher energy continuum above the energy scale $s_{0}$,
that should be determined via the sum-rule method. In this work we use the
unitarized chiral form factors Guo:2012ym ; Guo:2012yt  to construct the
resonant part ${\rm Im}\Pi^{\textrm{(Res)}}(s)$ of phenomenological spectral
function, which coherently include resonances and the two-meson continuum in
the low energy sector.

In the present convention, the spectral function or the imaginary part of the
two-point correlator is related to the two-meson form factors through

$\displaystyle{\rm
Im}\Pi^{\textrm{(Res)}}(s)=\frac{3m_{\pi}^{4}}{32\pi}\sum_{k=1}^{5}\frac{q_{k}(s)}{\sqrt{s}}\big{|}F_{k}^{\rm{S}}\big{|}^{2}\theta(s-s_{k}^{\rm{threshold}})\,,$
(8)

where the subscript $k$ runs over the relevant two-meson channels:
$\pi\pi,K\bar{K},\eta\eta,\eta\eta^{\prime}$ and $\eta^{\prime}\eta^{\prime}$,
$F_{k}^{\rm S}(s)$ denotes the isoscalar-scalar two-meson form factor of the
$k$th channel for the $j_{S}$ current in Eq. (2), and $q_{k}(s)$ stands for
the center of mass three-momentum of the two particles in the $k$th channel.
We will keep the same conventions of the form factors as those in Refs.
Guo:2012ym ; Guo:2012yt  and accordingly introduce other relevant coefficients
in Eq. (8) to match the normalization of the correlator of $j_{S}$ in Eqs. (1)
and (4). At low energies the intermediate $\pi\pi$ state is the most important
hadronic channel in the isoscalar scalar case. In the meanwhile, we also
include the $K\bar{K},\eta\eta,\eta\eta^{\prime}$ and
$\eta^{\prime}\eta^{\prime}$ channels in this work. The scalar form factors
$F_{PP^{\prime}}^{\rm S}(s)$ are calculated in the the framework of the
unitarized $U(3)$ resonance chiral theory Guo:2012ym ; Guo:2012yt . For the
sake of completeness, we recap the construction of the chiral phenomenological
spectrum function here. The leading order (LO) $U(3)$ chiral Lagrangian is

$\displaystyle\mathcal{L}^{({\rm LO})}=\frac{F^{2}}{4}\langle
u_{\mu}u^{\mu}\rangle+\frac{F^{2}}{4}\langle\chi_{+}\rangle+\frac{F^{2}}{3}M_{0}^{2}\ln^{2}{\det
u}\,,$ (9)

where the first two terms coincide with the conventional $SU(3)$ chiral
perturbation theory Gasser:1984gg  and the last one parameterizes the effect
of the QCD $U_{A}(1)$ anomaly, giving rise to the LO mass $M_{0}$ to the
singlet $\eta_{0}$ meson. The definitions of the basic chiral tensors
$u_{\mu}$ and $\chi_{+}$ can be found in Refs. Guo:2011pa ; Guo:2012ym ;
Guo:2012yt  and references therein. Regarding the higher order effects beyond
the LO, one could include the higher order local operators Gasser:1984gg  or
alternately one can also introduce explicitly the heavier degrees of freedom,
such as the bare resonance states Ecker:1988te . The latter method is employed
to calculate the two-meson scattering amplitudes and the form factors in Refs.
Guo:2011pa ; Guo:2012ym ; Guo:2012yt  within the framework of resonance chiral
theory Ecker:1988te . Three types of resonances, including scalar ($S$),
vector ($V$) and pseudoscalar ($P$) ones, are considered and the relevant
interacting Lagrangians read

$\displaystyle\mathcal{L}_{S}$
$\displaystyle=c_{d}\langle\,S_{8}u_{\mu}u^{\mu}\,\rangle+c_{m}\langle\,S_{8}\chi_{+}\,\rangle+\widetilde{c}_{d}S_{1}\langle\,u_{\mu}u^{\mu}\,\rangle+\widetilde{c}_{m}S_{1}\langle\,\chi_{+}\,\rangle\,,$
(10) $\displaystyle\mathcal{L}_{V}=\frac{iG_{V}}{2\sqrt{2}}\langle
V_{\mu\nu}[u^{\mu},u^{\nu}]\rangle\,,$ (11)
$\displaystyle\mathcal{L}_{P}=id_{m}\langle\,P_{8}\chi_{-}\,\rangle+i\widetilde{d}_{m}P_{1}\langle\,\chi_{-}\,\rangle\,.$
(12)

Two relevant local operators at the next-to-leading order that cannot receive
contributions after integrating the resonances in Eqs. (10), (11) and (12),
are also taken into account

$\displaystyle\mathcal{L}_{\Lambda}$
$\displaystyle=\frac{F^{2}\,\Lambda_{1}}{12}\partial^{\mu}X\partial_{\mu}X-\frac{F^{2}\,\Lambda_{2}}{12}X\langle\,\chi_{-}\,\rangle\,,\quad(X=\ln\det
U)\,.$ (13)

The chiral one-loop calculation with explicit tree-level resonance exchanges
in the $U(3)$ resonance chiral theory is completed in Ref. Guo:2011pa . Later
on the one-loop calculation of the strangeness conserving scalar and
pseudoscalar form factors is also carried out in a similar theoretical
framework in Refs. Guo:2012ym ; Guo:2012yt . An algebraic approximation of the
N/D method is used to perform the unitarization of the scattering amplitudes
and the scalar two-meson form factors. Under this assumption, the unitarized
scalar form factor takes the form

$\displaystyle\mathcal{F}(s)$ $\displaystyle=$
$\displaystyle\big{[}1+N(s)\cdot g(s)\big{]}^{-1}\cdot R(s)\,,$ (14)

with

$\displaystyle N(s)$ $\displaystyle=$ $\displaystyle{T(s)}^{\rm
LO+Res+Loop}+T(s)^{\rm LO}\cdot g(s)\cdot T(s)^{\rm LO}\,,$ (15)
$\displaystyle R(s)$ $\displaystyle=$ $\displaystyle{F(s)}^{\rm
LO+Res+Loop}+N(s)^{\rm LO}\cdot g(s)\cdot F(s)^{\rm LO}\,.$ (16)

Here $T(s)$ and $F(s)$ stand for the $S$-wave two-meson scattering amplitudes
and the scalar form factors from the perturbative calculations, respectively.
The superscripts LO, Res, Loop denote the contributions from the leading
order, resonance exchanges (also including the $\Lambda_{1}$ and $\Lambda_{2}$
terms) and the chiral loops, in order. By construction, the quantities $N(s)$
and $R(s)$ do not contain any right-hand cut. The latter effect is included in
the $g(s)$ function, which is given by

$\displaystyle g(s)$ $\displaystyle=$
$\displaystyle\frac{1}{16\pi^{2}}\bigg{[}a_{SL}(\mu)+\log\frac{m_{2}^{2}}{\mu^{2}}-x_{+}\log\frac{x_{+}-1}{x_{+}}-x_{-}\log\frac{x_{-}-1}{x_{-}}\bigg{]}\,,$
(17) $\displaystyle x_{\pm}$ $\displaystyle=$
$\displaystyle\frac{s+m_{1}^{2}-m_{2}^{2}}{2s}\pm\frac{1}{-2s}\sqrt{-4m_{1}^{2}s+(s+m_{1}^{2}-m_{2}^{2})^{2}}\,.$

In coupled-channel case, it becomes a diagonal matrix spanned in the
scattering channel space and for the $k$th channel one should replace $m_{1}$
and $m_{2}$ in Eq. (17) by the masses of the two particles in that channel.
The subtraction constants $a_{SL}$ in Eq. (17) are usually free parameters
that need to be fixed in the phenomenological fits. In our unitarization
procedure, all the unknown parameters appearing in the unitarized scalar form
factors also enter in the unitarized scattering amplitudes. As a consequence,
the phenomenological fits to the various experimental scattering observables
Guo:2011pa ; Guo:2012ym ; Guo:2012yt  enable one to fix all the unknown
parameters in the scalar form factors. The uncertainties of the resonance
couplings $G_{V}$, $c_{d}$, $c_{m}$, $\tilde{c}_{d}$, $\tilde{c}_{m}$, the
subtraction constants and also the bare masses of the resonances have been
estimated in Refs. Guo:2011pa ; Guo:2012yt  by using the error bars of the
various experimental data, including the phase shifts and inelasticities from
the $\pi\pi\to\pi\pi$, $\pi\pi\to K\bar{K}$, $K\pi\to K\pi$ processes with
different isospin and angular momentum numbers and also the $\pi\eta$ event
distributions. Such uncertainties of the relevant parameters can be also used
to make the errorbar analyses of the form factors and phenomenological
spectrum functions, which will be explicitly shown in the next section. The
spectral function used here clearly improves the one in Ref. Yuan:2017foa ,
with $\sigma$ and $f_{0}(980)$ included in the Breit-Wigner formalism, which
is likely to be improper for the broad resonance. An additional free parameter
is introduced as the weight of the two resonances in their spectral function.
In Ref. Cherry:2001cj , the phenomenological spectral density for the
isoscalar scalar channel is assumed to be solely contributed by the $\pi\pi$
channel. In our study, the heavier-mass channels, such as the
$K\bar{K},\eta\eta$, $\eta\eta^{\prime}$ and $\eta^{\prime}\eta^{\prime}$, are
explicitly included within the $U(3)$ chiral framework.

For the ESC contributions in the phenomenological spectral density in Eq. (7),
the traditional ESC model is utilized in this work Yuan:2017foa ; Shi:1999hm ,

$\begin{split}\frac{1}{\pi}{\rm
Im}\Pi^{\textrm{(ESC)}}(s)&=m_{q}^{2}(1/\sqrt{\tau})\left\\{\frac{3s}{8\pi^{2}}\left[1+\frac{\alpha_{s}}{\pi}\big{(}\frac{17}{3}-2\ln\big{(}s\tau\big{)}\big{)}+\big{(}\frac{\alpha_{s}}{\pi}\big{)}^{2}\big{(}31.8640\right.\right.\\\
&-\frac{95}{3}\ln\big{(}s\tau\big{)}+\frac{17}{4}\ln^{2}\big{(}s\tau\big{)}\big{)}+\big{(}\frac{\alpha_{s}}{\pi}\big{)}^{3}\big{(}89.1564-297.596\ln\big{(}s\tau\big{)}\\\
&\left.\left.+\frac{229}{2}\ln^{2}\big{(}s\tau\big{)}-\frac{221}{24}\ln^{3}\big{(}s\tau\big{)}\right]-\frac{3}{4\pi}sJ_{1}(\sqrt{s}\rho)Y_{1}(\sqrt{s}\rho)\right\\}\theta(s-s_{0})\,,\end{split}$
(18)

where $J_{1}$ and $Y_{1}$ are Bessel functions of the first and second kind,
respectively, and $s_{0}$ is the continuum threshold separating the
contributions from higher excited states and continuum. It should be noted
that we include in Eq. (18) the perturbative contribution from the quark-
hadron duality up to four-loop order Shi:1999hm .

By performing the Borel transformation of Eq.(7), one can obtain the
phenomenological representation for the correlation function

$\displaystyle R^{\textrm{(Phen)}}(\tau,s_{0},\hat{m}_{q})$ $\displaystyle=$
$\displaystyle\frac{1}{\pi}\int_{4m_{\pi}^{2}}^{\infty}e^{-s\tau}{\rm
Im}\Pi^{\textrm{(Phen)}}(s)\,ds$ (19) $\displaystyle=$ $\displaystyle
R^{\textrm{(Res)}}(\tau,s_{0})+R^{\textrm{(ESC)}}(\tau,s_{0},\hat{m}_{q})\,.$

The master equation for QCD sum rules can be established by matching the
results from the QCD calculation in Eq. (4) with its counter part of the
phenomenological one in Eq. (19), $i.e.$,

$R^{\textrm{(QCD)}}(\tau,\hat{m}_{q})=R^{\textrm{(Phen)}}(\tau,s_{0},\hat{m}_{q}).$
(20)

Within the SVZ QCD sum rules, in order to extract the light quark mass
$m_{q}$, we firstly derive the expression of $\hat{m_{q}}$ from Eq.(20). After
an analytical manipulation, we obtain

$\displaystyle\hat{m}_{q}=\frac{m_{\pi}^{2}\sqrt{3}}{8\pi}\sqrt{\frac{\int_{4m_{\pi}^{2}}^{s_{0}}e^{-s\tau}{\rm
Im}\Pi^{\textrm{(phen)}}(s)ds}{f^{2}(\tau)g(\tau)-\int_{s_{0}}^{\infty}e^{-s\tau}f^{2}(\tau)h(\tau,s)ds}}\,,$
(21)

where the functions $f(\tau),\,g(\tau)$, and $h(\tau,s)$ are defined as

$\displaystyle\begin{split}f(\tau)\equiv&\frac{1}{\left(\frac{1}{2}L\right)^{4/9}}\left\\{1+\frac{290}{729}\frac{1}{L}-\frac{256}{729}\frac{\ln
L}{L}+\left(\frac{550435}{1062882}-\frac{80\zeta(3)}{729}\right)\frac{1}{L^{2}}\right.\\\
&-\frac{388736}{531441}\frac{\ln
L}{L^{2}}+\frac{106496}{531441}\frac{\ln^{2}L}{L^{2}}+\left(\frac{2121723161}{2324522934}+\frac{8}{6561}\pi^{4}-\frac{119840}{531441}\zeta(3)\right.\\\
&\left.-\frac{8000}{59049}\zeta(5)\right)\frac{1}{L^{3}}+\left(-\frac{611418176}{387420489}+\frac{112640}{531441}\zeta(3)\right)\frac{\ln
L}{L^{3}}\\\
&\left.+\frac{335011840}{387420489}\frac{\ln^{2}L}{L^{3}}-\frac{149946386}{1162261467}\frac{\ln^{3}L}{L^{3}}\right\\}\,,\end{split}$
(22)
$\displaystyle\begin{split}g(\tau)\equiv&\frac{3}{8\pi^{2}}\frac{1}{\tau^{2}}\left\\{\left[1+4.821098\frac{\alpha_{s}(1/\tau)}{\pi}+21.97646\left(\frac{\alpha_{s}(1/\tau)}{\pi}\right)^{2}+53.14197\left(\frac{\alpha_{s}(1/\tau)}{\pi}\right)^{3}\right]\right.\\\
&+\frac{\langle\alpha_{s}G^{2}\rangle}{8\pi}\left(1+\frac{11}{2}\frac{\alpha_{s}(1/\tau)}{\pi}\right)+3\langle
m_{q}\bar{q}q\rangle\left(1+\frac{13}{3}\frac{\alpha_{s}(1/\tau)}{\pi}\right)\\\
&\left.-\frac{176}{27}\pi\kappa\alpha_{s}\langle\bar{q}q\rangle^{2}\left[\frac{\alpha_{s}(1/\tau)}{\alpha_{s}(\mu^{2}_{0})}\right]^{1/9}\tau+\frac{3}{8\pi^{2}}\frac{e^{\frac{-\rho^{2}}{2\tau}}\rho^{2}}{\tau^{3}}\left[K_{0}\left(\frac{\rho^{2}}{2\tau}\right)+K_{1}\left(\frac{\rho^{2}}{2\tau}\right)\right]\right\\}\,,\end{split}$
(23)
$\displaystyle\begin{split}h(\tau,s)\equiv&\left\\{\frac{3s}{8\pi^{2}}\left[1+\frac{\alpha_{s}(1/\tau)}{\pi}\big{(}\frac{17}{3}-2\ln\big{(}s\tau\big{)}\big{)}+\big{(}\frac{\alpha_{s}(1/\tau)}{\pi}\big{)}^{2}\big{(}31.8640-\frac{95}{3}\ln\big{(}s\tau\big{)}\right.\right.\\\
&+\frac{17}{4}\ln^{2}\big{(}s\tau\big{)}\big{)}+\big{(}\frac{\alpha_{s}(1/\tau)}{\pi}\big{)}^{3}\big{(}89.1564-297.596\ln\big{(}s\tau\big{)}\\\
&\left.\left.+\frac{229}{2}\ln^{2}\big{(}s\tau\big{)}-\frac{221}{24}\ln^{3}\big{(}s\tau\big{)}\big{)}\right]-\frac{3}{4\pi}sJ_{1}(\sqrt{s}\rho)Y_{1}(\sqrt{s}\rho)\right\\}\theta(s-s_{0})\,.\end{split}$
(24)

Then we can obtain the scale-running light-quark mass $m_{q}(1/\sqrt{\tau})$
from Eq.(6).

In this work, we will follow another version of the QCD sum rules, that is the
so-called Monte-Carlo-based QCD sum rules Leinweber:1995fn ; Wang:2011zzx , to
determine the threshold parameter $s_{0}$ and the RG-invariant light quark
mass $\hat{m}_{q}$. To be more specific, the values of $s_{0}$ and
$\hat{m}_{q}$ will be determined by minimizing the following quantity

$\chi^{2}_{k,\xi}=\sum_{j=1}^{n_{B}}\frac{[R^{\textrm{(QCD)}}_{k}(\tau_{j},\hat{m}_{q}^{k,\xi})-R^{\textrm{(phen)}}_{\xi}(\tau_{j},s_{0}^{k,\xi},\hat{m}_{q}^{k,\xi})]^{2}}{\sigma_{\text{QCD}}^{2}(\tau_{j})}\,.$
(25)

In order to quantitatively estimate the error bars of the resulting $s_{0}$
and $\hat{m}_{q}$ caused by the uncertainties from the input QCD parameters,
we generate large amount of $R^{\textrm{(QCD)}}_{k}$ by randomly varying the
OPE condensate parameters and the instanton size $\rho$ within ten percent
uncertainties and randomly sampling the $\kappa$ parameter in the range $2\sim
4$. Meanwhile, we also use the results from Ref. Guo:2012yt  to randomly
generate large samples of phenomenological spectral functions. The index $k$
in Eq. (25) corresponds to the $k$th configuration among the randomly
generated QCD-parameter samples, the index $\xi$ runs over different samples
of the phenomenological spectral functions. The Borel parameter
$\tau_{j}=\tau_{\textrm{min}}+(\tau_{\textrm{max}}-\tau_{\textrm{min}})\times(j-1)/(n_{B}-1)$,
runs over the appropriate window $[\tau_{\text{min}},\tau_{\text{max}}]$,
which range will be explained later in next section. We will take $n_{B}=21$
in Eq. (25) as in Ref. Leinweber:1995fn ; Wang:2011zzx . In fact, we have
tried to increase the values of $n_{B}$ in the calculation, and it turns out
that larger values of $n_{B}$ barely affect the results. The large samples
from the randomly generated $R^{\textrm{(QCD)}}_{k}$ and
$R^{\textrm{(phen)}}_{\xi}$ will then be implemented in Eq. (25) to repeatedly
determine the $s_{0}$ and $\hat{m}_{q}$. The large amount of the fitted
$s_{0}$ and $\hat{m}_{q}$ will be used for uncertainty studies.

The quantity $\sigma_{\text{QCD}}^{2}$ in Eq. (25) corresponds to the variance
of the $R^{\rm{(QCD)}}_{k=1,2,...n_{S}}$. In our study we generate a large
sample of $R^{\rm{(QCD)}}_{k=1,2,...n_{S}}$ with $n_{S}=20000$ at each Borel
parameter $\tau_{j}$. At the $j$th Borel parameter $\tau_{j}$ the variance
$\sigma_{\text{QCD}}^{2}$ takes the standard definition

$\displaystyle\sigma_{\text{QCD}}^{2}(\tau_{j})=\frac{1}{(n_{S}-1)}\sum_{k=1}^{n_{S}}[R_{k}^{\text{(QCD)}}(\tau_{j})-\overline{R^{\text{(QCD)}}}(\tau_{j})]^{2},$
(26)

where $k$ denotes the $k$th set of the randomly generated $n_{S}$ QCD-
parameter samples, and $\overline{R^{\text{(QCD)}}}(\tau_{j})$ corresponds to
the mean value of the huge samples $R^{\rm{(QCD)}}_{k=1,2,...n_{S}}$. It is
pointed out that we estimate the values of $\sigma_{\text{QCD}}^{2}(\tau_{j})$
by taking into account the explicit dependences of $R^{\textrm{(QCD)}}$ on the
OPE condensates, the instanton size $\rho$, and the $\kappa$ in Eq. (4). Since
the quantity $\hat{m}_{q}$ is implicitly dependent on the condensates and the
parameters $\rho$ and $\kappa$, we do not explicitly take the uncertainties of
$\hat{m}_{q}$ when estimating $\sigma_{\text{QCD}}^{2}$ and $\hat{m}_{q}$ will
be fixed at the value obtained by taking the central values of the relevant
input QCD parameters. To be more specific, we firstly calculate the random
configurations of the $R^{\text{(QCD)}}$ by assigning ten-percent
uncertainties to the OPE condensate parameters and the instanton size $\rho$,
and also generating random samples of $\kappa$ between 2 and 4. All the random
samples of the input parameters are generated with Gaussian distributions.
Afterwards, we evaluate the standard variance $\sigma_{\text{QCD}}^{2}$ of the
large data sets of $R^{\text{(QCD)}}$. Next, we use the
$\sigma_{\text{QCD}}^{2}$ in the expression of $\chi^{2}$ to fit the optimal
results $s_{0}$ and $\hat{m}_{q}$. The resulting curves of
$R^{\text{(QCD)}}(\tau)$, together with the uncertainties $\sigma_{\rm
QCD}(\tau_{j=1,2,...,21})$, are explicitly shown in Fig. 1.

Figure 1:  $R^{\text{(QCD)}}$ in Eq. (4), together with its uncertainties
$\sigma_{\rm QCD}(\tau_{j=1,2,...,21})$ in Eq. (26), as functions of the Borel
parameter $\tau$. See the text for details about the estimation of the
uncertainty $\sigma_{\rm QCD}$.

##  III Numerical results

For the numerical analysis, we follow Ref. Yuan:2017foa  to implement ten-
percent uncertainties for the OPE parameters evaluated at
$\mu_{0}=1\;\text{GeV}$, whose explicit values read

$\displaystyle\begin{aligned} &\langle\alpha_{s}G^{2}\rangle=(0.070\pm
0.007)\,\text{GeV}^{4}\;,\;\langle m_{q}\bar{q}q\rangle=-(0.10\pm
0.01)\,\text{GeV}^{4}\;,\\\
&\kappa\alpha_{s}\langle\bar{q}q\rangle^{2}=(1.49\pm 0.15)\times
10^{-4}\kappa\;\text{GeV}^{6}\;,\;\rho=(1.67\pm
0.17)\;\text{GeV}^{-1}\;,\end{aligned}$ (27)

where the factor $\kappa$ is introduced to quantify the violation of
factorization hypothesis in estimating the dimension-6 quark condensates. The
values of $\kappa$ have been determined in the range of $2\sim 4$ in different
channels Narison:2009vy ; Chung:1984gr ; Narison:1995jr . The influences of
the $\kappa$ parameter and the instanton on the OPE expansion, will be
discussed in a moment.

For the input phenomenological spectral functions, we will take the analyses
from Refs. Guo:2012yt ; Guo:2011pa . In the former reference two different
multiplets of bare scalar resonances and one multiplet of pseudoscalar
resonances are included in the study of the spectral functions, while in the
latter reference only one multiplet of bare scalar states are considered. As a
result, somewhat different values of the parameters are also obtained in Refs.
Guo:2011pa ; Guo:2012yt . In the following discussions, the values obtained by
using the phenomenological spectral functions of Ref. Guo:2012yt  will be
considered as the preferred central results. The deviations resulting from the
use of spectral functions of Ref. Guo:2011pa  will be considered as systematic
uncertainties to the central results. In addition, we also randomly generate
large samples of phenomenological spectral functions according to the results
from Ref. Guo:2012yt  (see also the discussions in the last section), which
will be used to estimate the statistical uncertainties of the $s_{0}$ and
$\hat{m}_{q}$. The phenomenological spectral functions after including the
Borel transformation factor $e^{-\tau s}$ are explicitly shown in Fig. 2. It
is obvious that the contributions in the energy region below 1 GeV are
dominant and the effects from the high energy tails beyond the 1.4 GeV are
tiny and negligible.

Figure 2: The phenomenological spectrum functions ${\rm Im}\Pi(s)$ multiplied
by the Borel transformation factor $e^{-\tau s}$ with $\tau=1.0$ GeV-2 as a
function of energy. The red line corresponds to the spectral function by using
the best fit parameters and the corresponding amplitudes given in Ref.
Guo:2012yt  and the surrounding gray areas are randomly generated according to
the uncertainties estimated in the same reference. The blue line stands for
the spectral function by taking the best fit parameters and related amplitudes
from Ref. Guo:2011pa . In both cases, the effects in the energy region beyond
1 GeV are found to be small.

###  III.1  Effects of the parameter $\kappa$ and the instanton in
$R^{\text{(QCD)}}$

The sizes of the dimension-six four-quark condensates have been intensively
discussed in the literature Narison:2009vy ; Chung:1984gr ; Narison:1995jr ;
Gimenez:1989bm ; Kremer:1984bg ; Braaten:1991qm ; Bertlmann:1984ih ;
Launer:1983ib ; Chetyrkin:1998yr ; Cherry:2001cj ; Dominguez:1984yx ;
Wang:2016sdt ; Wang:2021dcb , and there is still no consensus on their precise
values. It is plausible that the magnitudes of the four-quark condensates are
underestimated in the vacuum saturation assumption Shifman , meaning that the
$\kappa$ parameter can be obviously larger than one. Indeed a wide range of
$\kappa$ between 2 and 4 Chung:1984gr ; Narison:2009vy ; Narison:1995jr ;
Gimenez:1989bm ; Cherry:2001cj ; Wang:2016sdt , is determined from various
phenomenological processes, including the $e^{+}e^{-}\to hadrons$, the
$\tau\to hadron$ decays, the baryon spectra, etc. A recent analysis by combing
the electron-positron annihilation and the inclusive $\tau$ decay processes
gives the preferred value of $\kappa=3.0\pm 0.2$ Narison:2009vy , which fits
well in the wide range $2\sim 4$. In this work, we make a very conservative
estimate of the $\kappa$ by taking its value as $3.0\pm 1.0$, and assess its
influence on the determination of the light-quark mass.

In order to clearly show the effects of the $\kappa$, the instanton and other
nonperturbative terms in the QCD spectral function (4), we give in Fig. 3 the
ratios of the various contributions with respect to the perturbative QCD
result by taking $\tau=1.0\,\text{GeV}^{-2}$. It should be noted that the
contribution of the four-quark condensate is negative and we take its absolute
value in Fig. 3 in order to clearly show the sizes of different
nonperturbative contributions.

Figure 3: The ratios of the various nonperturbative contributions with respect
to the perturbative QCD part in $R^{\text{(QCD)}}$. The two vertical dashed
lines represent the conservative range of the parameter $\kappa$ used in our
study.

According to Fig. 3, it is obvious that the magnitude of the four-quark
condensate is larger than the two-quark and the two-gluon condensates, but is
smaller than the instanton contribution for $2<\kappa<4$. It is clear that the
contribution from the perturbative QCD term dominates $R^{\text{(QCD)}}$ in
Eq. (4). Furthermore, we also study the dependence of the various terms in the
OPE with the Borel parameter $\tau$, and the results are given in Fig. 4. To
be definite, we show the curves by taking $\kappa=3.0$ and verify that the
results by taking different values of $\kappa$ are quantitatively similar. As
in Fig. 3, the perturbative QCD term is the dominant part, while the
contributions from the nonperturbative parts, including the instanton, and the
quark and gluon condensates, slightly vary when the Borel parameter $\tau$ is
increased up to 1.2 GeV-2. The sensitivity of the parameter $\kappa$ is
related to the relative size of the four-quark condensate to the perturbative
QCD term, which will give part of the uncertainties of the Borel window
$\tau$, $s_{0}$ and $\hat{m}_{q}$ in our analysis. Detailed studies on the
error-bar analyses will be given in the next section.

Figure 4: The $\tau$ dependence of the ratios of various terms with respect to
the full $R^{\text{(QCD)}}$ in Eq. (4). The results by taking $\kappa=3.0$ are
explicitly given here. The curves by taking different values of $\kappa$
between 2.0 and 4.0 are quantitatively similar, and hence are not explicitly
shown.

###  III.2 Results from the conventional SVZ sum rules

In the conventional SVZ QCD sum rules, the parameters $\tau$ and $s_{0}$ are
fixed according to two well accepted criteria Shifman ; Reinders:1984sr ;
P.Col . The first one is required by the convergence of the OPE. One needs to
find a proper region for $\tau$ to guarantee that the contribution from the
highest dimension operator (HDO) included is less than $10\%$ in the OPE side.
The degree of convergence can be defined as

$\displaystyle\text{ratio}_{\langle\bar{q}q\rangle^{2}}=\frac{R^{\text{(OPE)}}_{\langle\bar{q}q\rangle^{2}}(\tau,\hat{m}_{q})}{R^{\text{(OPE)}}(\tau,\hat{m}_{q})}\,,$
(28)

where $R^{\text{(OPE)}}(\tau,\hat{m}_{q})$ represents the overall OPE
contributions in Eq.(4), and the subscript $\langle\bar{q}q\rangle^{2}$ refers
the HDO condensate here.

The second criterion requires that the portion of resonance contribution (RC),
$i.e.$, the contributions from the energy region below $s_{0}$, should be over
50% P.Col , which can be formulated as

$\displaystyle\text{ratio}^{\text{RC}}=\frac{R^{\text{(QCD)}}(\tau,\hat{m}_{q},s_{0})}{R^{\text{(QCD)}}(\tau,\hat{m}_{q})}\,,$
(29)

being $R^{\text{(QCD)}}(\tau,\hat{m}_{q},s_{0})$ the effects from the QCD side
below $s_{0}$. Under this prerequisite, the contributions of higher excited
and continuum states will be suppressed.

The standard procedure of the SVZ sum rules is to mutually vary the $s_{0}$
and $\tau$ to find their proper values, in order to satisfy the two
aforementioned criteria and find a smooth plateau for the $\hat{m}_{q}$. On
the smooth plateau, the RG-invariant light quark mass $\hat{m}_{q}$ should be
in principle independent of the Borel parameter $\tau$, or at least only shows
weak dependence. During this procedure, the preferred value of $s_{0}$ is
found to be 3.10 GeV2 for $\kappa=3.0$. We will tentatively assign a several-
percent uncertainty for the $s_{0}$ as $3.10\pm 0.20$ GeV2 in later
discussions Qiao:2013raa . The fulfillment of the two criteria and smooth
plateaus are shown in Fig. 5. The dependences of two ratios of Eqs. (28) and
(29) with the Borel parameter $\tau$ at the preferred $s_{0}=3.10$ GeV2 are
given in the left panel, from which one can easily identify the valid region
of $\tau$. In the right panel of Fig. 5, we show the rather weak dependence of
$\hat{m}_{q}$ with the Borel parameter $\tau$. By taking into account of the
error bars of the OPE condensates, the instanton size and the $\kappa$
parameter, we repeatedly take the standard procedure to determine the values
of $s_{0}$ and $\tau$. The results of the $s_{0}$ and $\tau$ with their
uncertainties extracted from the SVZ sum rules are summarized in Table 1.

Figure 5: (a) The OPE convergence $\text{ratio}_{\langle\bar{q}q\rangle^{2}}$ and the resonance contribution $\text{ratio}^{\text{RC}}$ as functions of the Borel parameter $\tau$ with $s_{0}=3.10\;\text{GeV}^{2}$ and $\kappa=3.0$. The blue and red lines represent $\text{ratio}_{\langle\bar{q}q\rangle^{2}}$ and $\text{ratio}^{\text{RC}}$, respectively. (b) The RG-invariant light-quark mass $\hat{m}_{q}$ as a function of the Borel parameter $\tau$ at different values of $s_{0}$ with $\kappa=3.0$. The black dots on each line indicate the lower/upper limits of the valid Borel window.  Table 1: Results from the SVZ QCD sum rules. The uncertainties of the Borel windows $\tau_{\rm min,max}$, $s_{0}$ and $\hat{m}_{q}$ are obtained by taking into account the error bars of the OPE parameters in Eq. (27) and $\kappa=3.0\pm 1.0$.  $\tau_{\text{min}}/\text{GeV}^{-2}$ | $\tau_{\text{max}}/\text{GeV}^{-2}$ | $s_{0}/\text{GeV}^{2}$ | $\hat{m}_{q}/\text{MeV}$ |  $m_{q}$ (2 GeV)/MeV
---|---|---|---|---
$0.37\pm 0.02$ | $1.05^{+0.19}_{-0.13}$ | $3.10\pm 0.20$ | $4.54^{+0.21}_{-0.29}$ | $3.46^{+0.16}_{-0.22}$

A slightly different treatment of the phenomenological spectral function from
Ref. Guo:2011pa , provides another source of uncertainty to $s_{0}$ and
$m_{q}$, compared to the results in Table 1 obtained by using the spectral
functions from Ref. Guo:2012yt . The parameters by taking the phenomenological
spectral functions from Ref. Guo:2011pa  with $\kappa=3.0$ are found to be

$\displaystyle s_{0}=2.90~{}{\rm GeV^{2}}\,,\quad\hat{m}_{q}=4.12~{}{\rm
MeV}\,,\quad m_{q}=3.13~{}{\rm MeV}\,.$ (30)

The deviations between the results in Eq. (30) and the central values in Table
1 are considered as systematic uncertainties caused by the different
treatments of spectral functions. In this way, our determination of the RG-
invariant light-quark mass $\hat{m}_{q}$ via the SVZ QCD sum rules reads

$\displaystyle\hat{m}_{q}=4.54^{+0.21}_{-0.29}\pm 0.42\,\text{MeV},$ (31)

where the first error bar stems from the uncertainties of Borel parameter
$\tau$, threshold parameter $s_{0}$, the condensates, the instanton size
$\rho$ and the factor $\kappa$, and the second one is originated by taking two
different sets of spectral functions Guo:2011pa ; Guo:2012yt . Similarly for
the scale dependent light-quark mass, our determination takes the value

$\displaystyle m_{q}(2\,\text{GeV})=3.46^{+0.16}_{-0.22}\pm
0.33\,\text{MeV}\,.$ (32)

###  III.3 Results from the Monte-Carlo-based sum rules

As an complementary study, we also analyze the RG-invariant light-quark mass
via the Monte-Carlo-based QCD sum rule, which is expected to give a more solid
quantitative estimation of the uncertainty. In this approach, the values of
$s_{0}$ and $\hat{m}_{q}$ are obtained by minimizing the $\chi^{2}$ functions
of Eq. (25) for a given $R^{\rm(QCD)}$ and $R^{\rm(phen)}$. According to the
results from the SVZ approach shown in Fig. 5, we will take a conservative
range for $0.37\,{\rm GeV}^{-2}<\tau<1.05\,{\rm GeV}^{-2}$ in the Monte-Carlo-
based sum rules. We have explicitly verified that the two criteria of the SVZ
sum rules in the conservative ranges of $\tau$ with the solutions of $s_{0}$
from the minimization procedures are generally well satisfied.

Figure 6: $\chi^{2}$ as a function of $s_{0}$ and $\hat{m}_{q}$, where the
central condensate values, the central phenomenological spectrum and
$\kappa=3.0$ are taken. A strong correlation between the $s_{0}$ and
$\hat{m}_{q}$ can be clearly seen.

Figure 7: The histogram distributions of the resulting parameters $s_{0}$ and
$m_{q}$ (obtained from the RG-invariant $\hat{m}_{q}$) from the repeated
$\chi^{2}$ minimization procedures by taking the large amount of random
samples of the OPE condensates, the instanton size $\rho$, the parameter
$\kappa$ and the spectral functions from Ref. Guo:2012yt .

To quantitatively estimate the uncertainties of $s_{0}$ and $\hat{m}_{q}$, we
generate huge amounts of random configurations of the $R^{\rm(QCD)}$ by
assigning ten-percent uncertainties to the OPE condensate parameters and the
instanton size $\rho$. The value of $\kappa$ is allowed to randomly float
between 2.0 and 4.0. In addition, we also generate large random samples of the
spectral functions according to the results from Ref. Guo:2012yt , as
explicitly shown in Fig. 2. In Fig. 6, we show the typical dependence of the
$\chi^{2}$ function with $s_{0}$ and $\hat{m}_{q}$ for a given set of
$R^{\rm(QCD)}$ and $R^{\rm(phen)}$. There exists a minimum valley of the
$\chi^{2}$ function in the $s_{0}$ and $\hat{m}_{q}$ plane, which implies the
clear correlation of the two parameters in a single minimization procedure.
Although the shape of the minimum valley sometimes looks moderately smooth, it
is not difficult to find the minimum point in the valley, which is exactly the
solutions of $s_{0}$ and $\hat{m}_{q}$ that we are looking for. Repeating the
minimization procedures by using the random configurations of $R^{\rm(QCD)}$
and $R^{\rm(phen)}$ gives us a huge amount of solutions for $s_{0}$ and
$\hat{m}_{q}$.

In Fig. 7, we show the histogram distributions of the resulting parameters
$s_{0}$ and $m_{q}$ (obtained from $\hat{m}_{q}$ via Eq. (6)), which seem
somewhat deviating the Gaussian-like distributions. After a scutinized
investigation, we find that the non-Gaussian distributions in Fig. 7 are
mainly caused by the $\rho$ parameter that describes the instanton size,
although the Gaussian uncertainties are assumed for the $\rho$ during the
error-bar analyses. However, by looking at the nontrivial depdence of the
$\rho$ parameter in $R^{\rm(QCD)}$ (4), it is not supprised to see the
somewhat non-Gaussian outpus in Fig. 7. The resulting values of the $s_{0}$,
$\hat{m}_{q}$ and $m_{q}$ with their statistical uncertainties contributed by
the OPE parameters, instanton size, the $\kappa$ parameter and spectral
functions from Ref. Guo:2012yt  obtained in the Monte-Carlo-based sum rules
are summarized in Table 2.

Table 2:  Results from the Monte-Carlo based QCD sum rules: $s_{0}$, $\hat{m}_{q}$ and the scale-dependence light quark masses $m_{q}$. A conservative Borel window $0.37\,{\rm GeV}^{-2}<\tau<1.05\,{\rm GeV}^{-2}$ is taken here. See the text for details. The statistical uncertainties are contributed by the OPE condensates, the instanton size $\rho$, the parameter $\kappa$ and the spectral functions of Ref. Guo:2012yt , $i.e.$, the shaded areas shown in Fig. 2.  $s_{0}/\text{GeV}^{2}$ | $\hat{m}_{q}/\text{MeV}$ |  $m_{q}$ (2 GeV)/MeV
---|---|---
$3.10\pm 0.11$ | $4.52\pm 0.18$ | $3.44\pm 0.14$

To make a conservative estimation of the uncertainties, we also take the
spectral function with the parameters and amplitudes from Ref. Guo:2011pa  to
determine $s_{0}$ and the light-quark masses. The explicit results are found
to be

$\displaystyle s_{0}=2.94~{}{\rm GeV^{2}}\,,\quad\hat{m}_{q}=4.10~{}{\rm
MeV}\,,\quad m_{q}=3.12~{}{\rm MeV}\,.$ (33)

As in the SVZ sum-rule case, the deviations by using the different spectral
functions calculated with the different parameters and amplitudes from Refs.
Guo:2011pa ; Guo:2012yt  are considered as systematic uncertainties, which can
be quadratically added to the statistical ones in Table 2. Therefore our final
determination of the $s_{0}$ is

$\displaystyle s_{0}=3.10\pm 0.11\pm 0.16~{}{\rm GeV^{2}}\,,$ (34)

where the first error bar is the statistical one from Table 2 and the second
error bar, obtained by taking the difference of the central values in Table 2
and Eq. (33), stands for the systematical uncertainty caused by different
spectral functions. Similar rules are also applied in the following
discussions. The combined result for the RG-invariant light-quark mass
$\hat{m}_{q}$ is

$\displaystyle\hat{m}_{q}=4.52\pm 0.18\pm 0.42~{}{\rm MeV}\,,$ (35)

and the scale-dependent light-quark mass obtained at 2 GeV reads

$\displaystyle m_{q}=3.44\pm 0.14\pm 0.32~{}{\rm MeV}\,,$ (36)

which perfectly agrees with the determination from the SVZ sum rule in Eq.
(32) and is also nicely consistent with the PDG value Zyla:2020zbs

$\displaystyle
m_{q}^{\text{PDG}}(2\,\text{GeV})=\frac{1}{2}(m_{u}+m_{d})=3.45^{+0.55}_{-0.15}\,\text{MeV}\,.$
(37)

As shown in Fig. 4, the instanton term gives the largest nonperturbative
contribution in the QCD spectral function. It is interesting to further study
its roles in the determinations of the light-quark masses. If the contribution
of the instanton is turned off, we obtain in the Monte-Carlo QCD sum rule
$s_{0}=2.28\text{GeV}^{2}$, $m_{q}=3.77$ MeV for $\kappa=3.0$, which can be
compared with the values in Eqs. (34) and (36) by taking all the
contributions. In the SVZ QCD sum rule approach, we obtain similar values
$s_{0}=2.25\text{GeV}^{2}$, $m_{q}=3.78$ MeV for $\kappa=3.0$ by simply
turning off the instanton contributions. It is clear that the instanton tends
to reduce the light-quark mass $m_{q}$ to the PDG average value.

Comparing with the results from the Breit-Wigner-type spectral functions used
in Ref. Yuan:2017foa , our predictions with more sophisticated scalar spectral
functions and higher order QCD corrections reduce the magnitudes of $m_{q}$
and turns out to be perfectly consistent with the PDG average value
Zyla:2020zbs , which is mainly based on several lattice determinations. It is
worthy pointing out that our results in Eqs. (32) and (36) obtained from the
sum rules with the scalar currents are compatible with the recent
determinations from the pseudoscalar channels Dominguez:2018azt .

##  IV Conclusions

In this work, we determine the average light-quark mass
$m_{q}=(m_{u}+m_{d})/2$ in the isoscalar scalar channel both from the SVZ and
the Monte-Carlo-based QCD sum rules, where sophisticated phenomenological
spectral functions from the chiral effective theory are used in our analysis.
The phenomenological spectral functions, that are determined from the
unitarized $U(3)$ chiral theory, coherently include the contributions from the
$f_{0}(500),f_{0}(980)$ and $f_{0}(1370)$, together with the two-meson
continuum effects. The present study gives an independent determination of the
light-quark mass from the scalar QCD sum rules, compared to the lattice
calculations and the sum rules in the pseudoscalar channels.

The prediction of the scale-dependent light-quark mass from the conventional
SVZ sum-rule method is $m_{q}(2\,\text{GeV})=(3.46^{+0.16}_{-0.22}\pm
0.33)\,\text{MeV}$, where the first part of the uncertainties are caused by
the OPE condensates, the instanton size, the $\kappa$ parameter, the Borel
parameter $\tau$ and the threshold parameter $s_{0}$, and the second part of
uncertainties are obtained by taking the two different sets of
phenomenological spectral functions. Similarly our prediction from the Monte-
Carlo-based QCD sum rule turns out to be $m_{q}(2\,\text{GeV})=(3.44\pm
0.14\pm 0.32)\,\text{MeV}$, which is nicely consistent with the result from
the SVZ approach and also remarkably agrees with the PDG value
$3.45_{-0.15}^{+0.55}$ MeV. Our work provides a solid proof that a proper way
to include $f_{0}(500)$ and $f_{0}(980)$ in the spectral functions is
important to obtain sensible results in the QCD sum rules. The present
theoretical formalism is also expected to be extended to the strange sector to
determine the strange-quark mass.

Acknowledgments

This work is supported in part by the National Natural Science Foundation of
China(NSFC) under the Grant Nos. 11605039, 11975090, 11575052, Natural Science
Foundation of Hebei Province with contract No. A2017205124, Science Foundation
of Hebei Normal University under Contract No. L2016B08, and the Fundamental
Research Funds for the Central Universities with contract No. 2242021R10099.
Fang-Hui Yin and Wen-Ya Tian are equally contributed to this work.

## References

  * (1) S. Aoki et al. [Flavour Lattice Averaging Group], Eur. Phys. J. C 80, no.2, 113 (2020). [arXiv:1902.08191 [hep-lat]]
  * (2) M.A. Shifman, A.I. Vainshtein and V.I. Zakharov, Nucl. Phys. B147, 385 (1979); ibid, Nucl. Phys. B147, 448 (1979).
  * (3) P. Colangelo and A. Khodjamirian, in At the frontier of particle physics / Handbook of QCD, edited by M. Shifman (World Scientific, Singapore, 2001), arXiv:hep-ph/0010175.
  * (4) J. I. Latorre and P. Pascual, J. Phys. G 11, L231 (1985).
  * (5) C. Becchi, S. Narison, E. de Rafael and F. J. Yndurain, Z. Phys. C 8, 335 (1981).
  * (6) J. Bijnens, J. Prades and E. de Rafael, Phys. Lett. B 348, 226 (1995).
  * (7) M. Jamin, J. A. Oller and A. Pich, Eur. Phys. J. C 24, 237-243 (2002).
  * (8) T. G. Steele, AIP Conf. Proc. 601, no.1, 151-157 (2001).
  * (9) R. Gupta, eConf C0304052, WG503 (2003) [arXiv:hep-ph/0311033 [hep-ph]].
  * (10) S. Bodenstein, C. A. Dominguez and K. Schilcher, JHEP 07, 138 (2013).
  * (11) S. Narison, Phys. Lett. B 738, 346 (2014).
  * (12) C. A. Dominguez, A. Mes and K. Schilcher, JHEP 1902, 057 (2019).
  * (13) S. N. Cherry and M. R. Pennington, hep-ph/0111158.
  * (14) S. Narison, hep-ph/0202200.
  * (15) K. Maltman, Nucl. Phys. B Proc. Suppl. 109, 178-182 (2002).
  * (16) C. A. Dominguez and E. de Rafael, Annals Phys. 174, 372 (1987).
  * (17) L. J. Reinders, S. Yazaki and H. R. Rubinstein, Nucl. Phys. B 196, 125-146 (1982).
  * (18) J. M. Yuan, Z. F. Zhang, T. G. Steele, H. Y. Jin and Z. R. Huang, Phys. Rev. D 96, no. 1, 014034 (2017).
  * (19) C. M. Shakin and H. Wang, Phys. Rev. D 64, 094020 (2001).
  * (20) V. Elias, A. H. Fariborz, F. Shi and T. G. Steele, Nucl. Phys. A 633, 279-311 (1998).
  * (21) F. Shi, T. G. Steele, V. Elias, K. B. Sprague, Y. Xue and A. H. Fariborz, Nucl. Phys. A 671, 416-446 (2000).
  * (22) Z. H. Guo, J. A. Oller and J. Ruiz de Elvira, Phys. Lett. B 712, 407-412 (2012).
  * (23) Z. H. Guo, J. A. Oller and J. Ruiz de Elvira, Phys. Rev. D 86, 054006 (2012).
  * (24) Z. H. Guo and J. A. Oller, Phys. Rev. D 84, 034005 (2011).
  * (25) D. B. Leinweber, Annals Phys. 254, 328 (1997).
  * (26) M. Benmerrouche, G. Orlandini and T. G. Steele, Phys. Lett. B 356, 573 (1995).
  * (27) Q. N. Wang, Z. F. Zhang, T. G. Steele, H. Y. Jin and Z. R. Huang, Chin. Phys. C 41, no. 7, 074107 (2017).
  * (28) S. Narison and G. Veneziano, Int. J. Mod. Phys. A 4, 2751 (1989).
  * (29) S. Narison, Nucl. Phys. B Proc. Suppl. 96, 244-251 (2001).
  * (30) R. L. Jaffe, Phys. Rev. D 15, 267 (1977).
  * (31) R. L. Jaffe, Phys. Rept. 409, 1-45 (2005).
  * (32) N. N. Achasov, S. A. Devyanin and G. N. Shestakov, Sov. J. Nucl. Phys. 32, 566 (1980).
  * (33) S. Narison, Phys. Rev. D 73, 114024 (2006).
  * (34) T. V. Brito, F. S. Navarra, M. Nielsen and M. E. Bracco, Phys. Lett. B 608, 69-76 (2005).
  * (35) Z. G. Wang, W. M. Yang and S. L. Wan, J. Phys. G 31, 971-979 (2005).
  * (36) H. J. Lee and N. I. Kochelev, Phys. Lett. B 642, 358-365 (2006).
  * (37) H. X. Chen, A. Hosaka and S. L. Zhu, Phys. Lett. B 650, 369-372 (2007).
  * (38) H. X. Chen, A. Hosaka and S. L. Zhu, Phys. Rev. D 76, 094025 (2007).
  * (39) E. V. Shuryak, Nucl. Phys. B 214, 237-252 (1983).
  * (40) V. Elias, F. Shi and T. G. Steele, J. Phys. G 24, 267-269 (1998).
  * (41) K. G. Chetyrkin, D. Pirjol and K. Schilcher, Phys. Lett. B 404, 337-344 (1997).
  * (42) S. Narison, Phys. Lett. B 673, 30 (2009).
  * (43) J. Gasser and H. Leutwyler, Nucl. Phys. B 250, 465-516 (1985).
  * (44) G. Ecker, J. Gasser, A. Pich and E. de Rafael, Nucl. Phys. B 321, 311-342 (1989).
  * (45) L. Wang and F. X. Lee, Comput. Phys. Commun. 182, 1721 (2011).
  * (46) Y. Chung, H. G. Dosch, M. Kremer and D. Schall, Z. Phys. C 25, 151 (1984).
  * (47) S. Narison, Phys. Lett. B 361, 121 (1995).
  * (48) Z. G. Wang, [arXiv:2107.09328 [hep-ph]].
  * (49) E. Braaten, S. Narison and A. Pich, Nucl. Phys. B 373, 581-612 (1992).
  * (50) G. Launer, S. Narison and R. Tarrach, Z. Phys. C 26, 433-439 (1984).
  * (51) V. Gimenez, J. A. Penarrocha and J. Bordes, Phys. Lett. B 223, 245-250 (1989).
  * (52) M. Kremer, N. A. Papadopoulos and K. Schilcher, Phys. Lett. B 143, 476-480 (1984).
  * (53) C. A. Dominguez, M. Kremer, N. A. Papadopoulos and K. Schilcher, Z. Phys. C 27, 481 (1985).
  * (54) R. A. Bertlmann, G. Launer and E. de Rafael, Nucl. Phys. B 250, 61-108 (1985).
  * (55) K. G. Chetyrkin, S. Narison and V. I. Zakharov, Nucl. Phys. B 550, 353-374 (1999) [arXiv:hep-ph/9811275 [hep-ph]].
  * (56) L. J. Reinders, H. Rubinstein and S. Yazaki, Phys. Rept. 127, 1 (1985).
  * (57) C. F. Qiao and L. Tang, Eur. Phys. J. C 74, 3122 (2014);ibid, Europhys. Lett. 107, 31001 (2014).
  * (58) P. A. Zyla et al. [Particle Data Group], PTEP 2020, no.8, 083C01 (2020).