Present address:] School of Physics and Astronomy, University of Glasgow,
Glasgow G12 8QQ, United Kingdom.

TUMQCD Collaboration

# Relations between Heavy-light Meson and Quark Masses

N. Brambilla  Physik-Department, Technische Universität München, James-Franck-
Straße 1, 85748 Garching, Germany  Institute for Advanced Study, Technische
Universität München, Lichtenbergstraße 2a, 85748 Garching, Germany  J.
Komijani  [  Physik-Department, Technische Universität München, James-Franck-
Straße 1, 85748 Garching, Germany  Institute for Advanced Study, Technische
Universität München, Lichtenbergstraße 2a, 85748 Garching, Germany  A.S.
Kronfeld  Institute for Advanced Study, Technische Universität München,
Lichtenbergstraße 2a, 85748 Garching, Germany  Theoretical Physics Department,
Fermi National Accelerator Laboratory,
Batavia, IL 60510, USA  A. Vairo  Physik-Department, Technische Universität
München, James-Franck-Straße 1, 85748 Garching, Germany

###### Abstract

The study of heavy-light meson masses should provide a way to determine
renormalized quark masses and other properties of heavy-light mesons. In the
context of lattice QCD, for example, it is possible to calculate hadronic
quantities for arbitrary values of the quark masses. In this paper, we address
two aspects relating heavy-light meson masses to the quark masses. First, we
introduce a definition of the renormalized quark mass that is free of both
scale dependence and renormalon ambiguities, and discuss its relation to more
familiar definitions of the quark mass. We then show how this definition
enters a merger of the descriptions of heavy-light masses in heavy-quark
effective theory and in chiral perturbation theory ($\chi$PT). For practical
implementations of this merger, we extend the one-loop $\chi$PT corrections to
lattice gauge theory with heavy-light mesons composed of staggered fermions
for both quarks. Putting everything together, we obtain a practical formula to
describe all-staggered heavy-light meson masses in terms of quark masses as
well as some lattice artifacts related to staggered fermions. In a companion
paper, we use this function to analyze lattice-QCD data and extract quark
masses and some matrix elements defined in heavy-quark effective theory.

††preprint:  FERMILAB-PUB-17/490-T TUM-EFT 105/17

##  I Introduction

Six of the fundamental parameters of the Standard Model are quark masses, but
because of confinement they can be determined only indirectly via their
influence on hadrons. In order to develop precise predictions of the Standard
Model, it is necessary to define the quark mass in a theoretically sound way,
based on quantum field theory, and establish relations between such
definitions and (nonperturbative) quantities that can be measured in
experiments or computed in lattice gauge theory. In this and a companion paper
Bazavov  _et al._ (2017a), we focus on the connection to lattice QCD and, in
particular, on the relation between the mass of a heavy-light hadron and the
mass, $m_{Q}$, of the heavy quark bound inside. Here, “heavy” refers to masses
much larger than the scale of QCD: $m_{Q}\gg\Lambda_{\text{QCD}}$.

Properties of heavy-light hadrons can be studied with two approximate
symmetries: chiral symmetry for the light quarks and spin-flavor symmetry for
the heavy quarks. Chiral symmetry is a symmetry of the QCD Lagrangian in the
limit of massless light quarks, while spin-flavor symmetry is a symmetry of
the limit of infinitely massive quarks. These symmetries are approximate but
very pertinent in suitable kinematic regions. The small effects of symmetry
breaking can be derived with effective field theories: chiral perturbation
theory ($\chi$PT) for chiral symmetry and heavy-quark effective theory (HQET)
for heavy-quark spin-flavor symmetry.

$\chi$PT and HQET can be combined into heavy-meson chiral perturbation theory
(HM$\chi$PT) Burdman and Donoghue (1992); Yan  _et al._ (1992); Wise (1992),
which describes the low-momentum interactions of pseudo-Goldstone bosons
($\pi$, $K$ and $\eta$) and hadrons containing a heavy quark. Let us use
$B_{x}$ to denote a generic heavy-light pseudoscalar meson composed of a light
quark $x$ and a heavy antiquark $\bar{Q}$, and $B^{*}_{x}$ to denote the
corresponding vector meson. The masses of the $B_{x}$ and $B^{*}_{x}$ mesons,
to first order in $1/m_{Q}$ and light-quark masses, are Wise (1992); Boyd and
Grinstein (1995)

$M_{B^{(*)}_{x}}=M_{0}+2\lambda_{1}B_{0}m_{x}+2\lambda^{\prime}_{1}B_{0}(m_{u}+m_{d}+m_{s})-d_{B^{(*)}}\frac{3\lambda_{2}}{2m_{Q}},$
(1)

where $m_{x}$ is the light valence mass; $m_{u}$, $m_{d}$, and $m_{s}$ are the
light (sea) quark masses; and $m_{Q}$ is the heavy antiquark mass.111Note that
all quark masses can vary in lattice-QCD simulations. The physical meaning of
the terms in Eq. (1) is as follows. The mass $M_{0}$ is the spin-averaged
heavy-light meson mass in the chiral limit. The next two terms parameterize
the light-quark-mass dependence; here $B_{0}$ enters the leading-order
relation between the pion mass and the quark masses,
$m_{\pi}^{2}=B_{0}(m_{u}+m_{d})$. The last term describes the leading spin-
dependent effect; here $d_{B}=1$ for pseudoscalar mesons and
$d_{B^{*}}=-{\textstyle\frac{1}{3}}$ for vector mesons. In the language of
HM$\chi$PT, $\lambda_{1}$, $\lambda^{\prime}_{1}$, and $\lambda_{2}$ are low
energy constants; $\lambda_{1}$ and $\lambda^{\prime}_{1}$ are independent of
the heavy-quark mass, while $\lambda_{2}$ has a logarithmic dependence on
$m_{Q}$ via $\alpha_{s}(m_{Q})$ Wise (1992), as elaborated below. These
quantities can be determined by matching HM$\chi$PT to the underlying
theories, namely HQET and ultimately QCD.

In HQET, the mass of the $B_{x}$ (or $B^{*}_{x}$) meson is Falk and Neubert
(1993)

$M_{B_{x}^{(*)}}=m_{Q}+{\overline{\Lambda}}+\frac{\mu_{\pi}^{2}}{2m_{Q}}-d_{B^{(*)}}\frac{\mu_{G}^{2}(m_{Q})}{2m_{Q}},$
(2)

to first order in $1/m_{Q}$. Each term here has a simple physical
interpretation: $\bar{\Lambda}$ is the energy of the light quarks and gluons,
$\mu_{\pi}^{2}/2m_{Q}$ is the kinetic energy of the heavy quark in the meson’s
rest frame, and $d_{B^{(*)}}\mu_{G}^{2}(m_{Q})/2m_{Q}$ corresponds to the
hyperfine energy of the interaction between the heavy quark’s spin with the
chromomagnetic field inside the meson. Although each quantity corresponds to a
term in the effective Lagrangian of HQET Kronfeld (2000), for Eq. (2) to be
useful in practice, $m_{Q}$ and ${\overline{\Lambda}}$ require a careful
definition. The principal new result of this paper is to propose a new
definition with a natural interpretation in HQET, which in turn hinges on a
recent, improved understanding of the overall normalization of the renormalon
Komijani (2017). This definition is renormalon-free and scale-independent; we
call it the minimal renormalon-subtracted (MRS) mass.

Note that the kinetic term is free of renormalization, as a consequence of
reparametrization invariance Luke and Manohar (1992); Brambilla  _et al._
(2003); Heinonen  _et al._ (2012). On the other hand,

$\mu_{G}^{2}(m_{Q})=C_{\text{cm}}(m_{Q})\tilde{\mu}_{G}^{2},$ (3)

where $C_{\text{cm}}(m_{Q})$, which is known to three loops Grozin  _et al._
(2008), is the Wilson coefficient of the chromomagnetic operator obtained by
matching HQET to QCD at scale $m_{Q}$, and $\tilde{\mu}_{G}^{2}$ is the matrix
element of the chromomagnetic interaction operator in HQET.

The HQET quantities on the right hand side of Eq. (2) depend on masses of
light quarks. Moreover, the heavy-quark mass $m_{Q}$ depends on the sea-quark
masses via two- and higher-loop diagrams. In this work, we exploit $\chi$PT to
incorporate contributions from nonvanishing masses of light quarks. If we set
all light quark masses to zero in Eq. (2) and compare to Eq. (1), then we find

$\displaystyle M_{0}$
$\displaystyle=\left[m_{Q}+{\overline{\Lambda}}+\frac{\mu_{\pi}^{2}}{2m_{Q}}+\text{O}\left(m_{Q}^{-2}\right)\right]_{\text{chiral
limit}},$ (4) $\displaystyle\lambda_{2}$
$\displaystyle=\left[\frac{1}{3}\mu_{G}^{2}(m_{Q})\right]_{\text{chiral
limit}}.$ (5)

The logarithmic dependence of $\lambda_{2}$ on $m_{Q}$ is now manifest via Eq.
(3).

For applications in lattice QCD, it is helpful to modify $\chi$PT to include
lattice discretization errors, especially since lattice-QCD computations often
involve an extrapolation in light-quark masses to the physical up and down
masses, and always require a continuum extrapolation in lattice spacing. With
staggered fermions, the appropriate extension of HM$\chi$PT is known as heavy-
meson, rooted, all-staggered chiral perturbation theory (HMrAS$\chi$PT)
Bernard and Komijani (2013). In this paper, we calculate one-loop corrections
to heavy-light meson masses within the framework of HMrAS$\chi$PT. Together
with the MRS mass, these corrections are needed to analyze lattice-QCD data
for masses of heavy-light mesons composed of staggered quarks, as in Ref.
Bazavov  _et al._ (2017a).

This paper is organized as follows. In Sec. II, we introduce the MRS mass.
Then, in Sec. III, we present the HMrAS$\chi$PT corrections to the heavy-light
meson masses. Section IV combines the MRS mass with these expressions to
provide a practical way to analyze lattice-QCD data for heavy-light meson
masses, which is implemented in a companion paper Bazavov  _et al._ (2017a) to
extract the quark masses and the HQET matrix elements ${\overline{\Lambda}}$
and $\mu_{\pi}^{2}$. Finally, in Sec. V, we summarize the results and offer
some outlook, both to the removal of higher renormalons and to applications in
quarkonium.

##  II Minimal-renormalon-subtracted mass

Let us return to the HQET description of heavy-light meson masses:

$M_{B^{(*)}}=m_{Q}+{\overline{\Lambda}}+\frac{\mu_{\pi}^{2}}{2m_{Q}}-d_{B^{(*)}}\frac{\mu_{G}^{2}(m_{Q})}{2m_{Q}}+\cdots.$
(6)

Clearly, the meaning of ${\overline{\Lambda}}$ is connected to the choice of
the renormalized mass $m_{Q}$, and similar connections persist to the higher
orders in the $1/m_{Q}$ expansion. In the context of HQET, the most natural
choice for $m_{Q}$ is the pole mass, because the pole mass does not depend on
the ultraviolet regulator and therefore has the same physical interpretation
in QCD and HQET. In QCD and other non-Abelian gauge theories, the pole mass is
infrared finite Kronfeld (1998) and gauge independent Kronfeld (1998);
Breckenridge  _et al._ (1995) at every order in perturbation theory, but the
perturbative series relating it to a short-distance mass diverges, owing to
renormalon effects Bigi  _et al._ (1994); Beneke and Braun (1994). Indeed, the
leading infrared renormalon implies an intrinsic ambiguity of order
${\Lambda_{\text{QCD}}}$ in the definition of the pole mass. The meson mass is
unambiguous, so this ambiguity must be canceled by ${\overline{\Lambda}}$,
which itself is of order $\Lambda_{\text{QCD}}$. Similar ambiguities from
subleading infrared renormalons link the pole mass in Eq. (6) to the higher-
order terms Bigi  _et al._ (1994); Beneke and Braun (1994), such as
$\mu_{\pi}^{2}$.

Short-distance masses are free of these long-distance effects and are also
candidates for $m_{Q}$. In interpreting Eq. (6), they all leave something to
be desired. Consider choosing a mass-independent renormalization scheme, such
as $\overline{\rm MS}$, for $m_{Q}$. Then, a term of order $\alpha_{s}m_{Q}$
must be absorbed into $\overline{\Lambda}$, which spoils the $1/m_{Q}$ power
counting of HQET Luke  _et al._ (1995). Threshold masses have been introduced
to alleviate this problem. Uraltsev’s kinetic mass Uraltsev (1997, 2002), the
renormalon-subtracted (RS) mass Pineda (2001), and the MSR mass Hoang  _et
al._ (2008) all introduce a new factorization scale. The 1S mass Hoang (2000)
is based on quarkonium and the potential-subtracted mass Beneke (1998) on the
heavy-quark potential; they also introduce (at least) one more scale into the
study of heavy-light hadrons. Moreover, the heavy-quark potential is infrared
divergent starting at three loops Brambilla  _et al._ (1999).

Here, we propose a new threshold mass that avoids introducing an additional
scale. Our proposed mass is closely related to the RS mass, and we call it the
minimal renormalon-subtracted (MRS) mass. The RS mass removes the leading-
renormalon contribution from the pole mass along with a finite piece that
depends on a renormalization point, $\nu_{f}$, such that (ideally)
${\Lambda_{\text{QCD}}}\ll\nu_{f}\ll m_{Q}$. In the MRS mass, the same goal is
achieved without a new $\nu_{f}$. We advocate formulating HQET using the MRS
mass for the heavy quark, because it retains the advantages of the pole mass,
while circumventing its ambiguity. (We shall also comment on variations on the
MRS mass with similar advantages.) In the construction of HQET as an effective
theory of QCD, the shift from the pole to the MRS mass can be understood as an
application of the notion of a residual mass term Falk  _et al._ (1992).
Shifting the (pole) mass by an amount $\delta m$, $\delta m\ll m_{Q}$, does
not disrupt the matching of HQET to QCD Neubert and Sachrajda (1995); Luke
_et al._ (1995).

Because the pole mass of a quark is infrared finite Kronfeld (1998) and gauge
independent Kronfeld (1998); Breckenridge  _et al._ (1995) at each order in
perturbation theory, at any finite $N$ one can write

$m_{\text{pole}}={\overline{m}}\left(1+\sum_{n=0}^{N}r_{n}\alpha_{s}^{n+1}({\overline{m}})+\text{O}(\alpha_{s}^{N+2})\right),$
(7)

where ${\overline{m}}$ denotes the ${\overline{\rm MS}}$-renormalized mass
self-consistently evaluated at the scale $\mu={\overline{m}}$. Here,
$\alpha_{s}$ is the strong running coupling in a generic mass-independent
scheme; the values of the coefficients $r_{n}$ depend on this scheme, which is
implicit but obvious in the following equations.

At large orders, the coefficients $r_{n}$ grow like $n!$, so the series
diverges for all $\alpha_{s}>0$. To cope with the divergence, one can employ
the method of Borel resummation (Kawai and Takei, 2005, e.g.) to interpret a
divergent series and replace it with a better-behaved expression. The Borel
sum of a series involves an integration in the Borel plane from a base point
to infinity. The integration path is usually defined on the real axis from the
origin to $+\infty$. But, in general, it is possible to have different
integration paths that cannot be deformed to each other. Particularly, if
there are singularities on the positive real axis in the Borel plane, one
cannot attribute an unambiguous Borel sum to the divergent series.

In the case of the pole mass, one finds an ambiguity of order
${\Lambda_{\text{QCD}}}$ stemming from an infrared renormalon singularity.
Apart from the overall normalization, the full nature of this singularity,
and, hence, the large-$n$ behavior of $r_{n}$, stems from the running of the
gauge coupling Beneke (1995). For the normalization, expansions based on the
$r_{n}$ in Eq. (7) have been proposed Pineda (2001); Lee (2003); Hoang  _et
al._ (2008); Lee (2015); Komijani (2017); Hoang  _et al._ (2017); Pineda
(2017). Because the leading infrared renormalon in the pole mass is closely
related to an ultraviolet renormalon in the self-energy of a static quark
Beneke and Braun (1994); Beneke (1995) and, similarly, the heavy-quark
potential Pineda (1998); Hoang  _et al._ (1999), the leading renormalon is
independent of $\overline{m}$. Reinterpreting this observation, Ref. Komijani
(2017) shows that, up to an overall multiplicative factor, the large-$n$
behavior of $r_{n}$ can be obtained from the recurrence relation

$a_{n}=2\left[\beta_{0}na_{n-1}+\beta_{1}(n-1)a_{n-2}+\cdots+\beta_{n-1}a_{0}\right],\quad
n\geq 1,$ (8)

where $\beta_{n}$ are the coefficients of the beta function

$\beta(\alpha_{s})=\frac{d\alpha_{s}}{d\ln\mu^{2}}=-\alpha_{s}^{2}\left(\beta_{0}+\beta_{1}\alpha_{s}+\beta_{2}\alpha_{s}^{2}+\cdots\right),$
(9)

capturing the running of the coupling with renormalization point $\mu$.

The recurrence relation can be solved in closed form in the scheme with

$\beta\left(\alpha_{\text{g}}(\mu)\right)=-\frac{\beta_{0}\alpha_{\text{g}}^{2}(\mu)}{1-(\beta_{1}/\beta_{0})\alpha_{\text{g}}(\mu)}.$
(10)

We shall call the scheme implied by Eq. (10) “geometric” because it
corresponds to changing the scheme of $\alpha_{s}$ such that
$\beta_{j}/\beta_{0}=(\beta_{1}/\beta_{0})^{j}$, $j\geq 0$, and summing the
geometric series. We denote the coupling in this scheme $\alpha_{\text{g}}$
and comment on other schemes below. The geometric scheme, which is very
convenient in reducing the algebra, has been employed before in the literature
on renormalons (Brown  _et al._ , 1992; Lee, 1997, e.g.).

In this scheme, the solution of Eq. (8) is (Komijani, 2017, Eq. (2.14)),

$a_{n}=(2\beta_{0})^{n}\frac{\Gamma(n+1+b)}{\Gamma(2+b)}a_{0},\quad n\geq 1,$
(11)

where $b=\beta_{1}/(2\beta_{0}^{2})$, and $a_{0}$ is the initial term in the
sequence.222Note that Eq. (11) is _not_ valid for $n=0$. Without this gap in
the $a_{n}$, the sequence would not solve Eq. (8). In a general scheme, the
solution would include further terms that are suppressed by powers of $1/n$;
see, for instance, Ref. (Komijani, 2017, Eq. (2.5)). In the geometric scheme,
however, these corrections all vanish.

Equation (10) on its own does not suffice to define the geometric scheme. In
addition, one must specify the finite part dropped in the course of one-loop
renormalization from the bare to the renormalized coupling. Equivalently, we
can choose a value for $b_{1}$ in the relation

$\frac{1}{\alpha_{\text{g}}(\mu)}=\frac{1}{\alpha_{\overline{\rm
MS}}(\mu)}+b_{1}+b_{2}\alpha_{\overline{\rm MS}}(\mu)+\cdots.$ (12)

We choose $b_{1}=0$, which imposes a close relation between the geometric and
$\overline{\rm MS}$ schemes, as we shall see in Eq. (26), below. The other
constants $b_{n}$, $n>1$, are chosen so that the scheme-dependent part of the
beta function takes the geometric form, Eq. (10):
$b_{2}=\beta_{2}/\beta_{0}-(\beta_{1}/\beta_{0})^{2}$,
$b_{3}={\textstyle\frac{1}{2}}[\beta_{3}/\beta_{0}-(\beta_{1}/\beta_{0})^{3}]$,
etc.

The value of the overall normalization $a_{0}$ is needed in order to specify
that $r_{n}\sim a_{n}$ as $n\to\infty$. In the literature, it is more common
to work with $R_{0}=a_{0}/(1+b)$. (Elsewhere, $R_{0}$ is denoted $N_{m}$
Pineda (2001) or $N$ Beneke (1995); Komijani (2017).) Thus, we work with the
sequence

$R_{n}\equiv R_{0}(2\beta_{0})^{n}\frac{\Gamma(n+1+b)}{\Gamma(1+b)},\quad
n\geq 0,$ (13)

even though it is Eq. (11)—and not Eq. (13)—that solves the recurrence
relation. Reference (Komijani, 2017, Eq. (2.22)) provides a series solution
for $R_{0}$ in the geometric scheme:

$R_{0}=\sum_{k=0}^{\infty}r^{\prime}_{k}\frac{\Gamma(1+b)}{\Gamma(2+k+b)}\frac{1+k}{(2\beta_{0})^{k}},$
(14)

where

$r^{\prime}_{k}=r_{k}-2\left[\beta_{0}kr_{k-1}+\beta_{1}(k-1)r_{k-2}+\cdots+\beta_{k-1}r_{0}\right]$
(15)

with the $r_{k}$ of Eq. (7). Note that the $r^{\prime}_{k}$ are free of the
leading pole-mass renormalon in any scheme Komijani (2017), but when
evaluating $r^{\prime}_{k}$ for Eq. (14) the geometric scheme must be used for
the $r_{k}$. Because the leading renormalon cancels in the $r^{\prime}_{k}$,
they grow slowly enough such that the series in Eq. (14) converges.

To tame via Borel summation the leading-renormalon divergence on the right-
hand side of Eq. (7), the pole mass can be expressed as

$m_{\text{pole}}={\overline{m}}+{\overline{m}}\int_{0}^{\infty}dt\,e^{-t/\alpha_{\text{g}}(\mu)}\sum_{n=0}^{\infty}\left[r_{n}-R_{n}+R_{n}\right]\frac{t^{n}}{n!},$
(16)

where the $r_{n}$ are in the geometric scheme. In the following, our
manipulations overlook the subleading renormalons and other sources of
divergence (such as instanton contributions), taking the attitude that the
method could be generalized to handle them, and-or their practical
implications are not important. With this proviso in mind, we rewrite the pole
mass as

$\displaystyle m_{\text{pole}}$
$\displaystyle={\overline{m}}\left(1+\sum_{n=0}^{\infty}\left[r_{n}-R_{n}\right]\alpha_{\text{g}}^{n+1}({\overline{m}})\right)+\mathcal{J}({\overline{m}}),$
(17) $\displaystyle\mathcal{J}(\mu)$
$\displaystyle\equiv\mu\int_{0}^{\infty}dt\,e^{-t/\alpha_{\text{g}}(\mu)}\sum_{n=0}^{\infty}\frac{R_{n}}{n!}t^{n}.$
(18)

In Eq. (17), the power series in $\alpha_{\text{g}}$ with $[r_{n}-R_{n}]$ is
better-behaved than that in Eq. (7), because it is free of the leading
infrared renormalon Pineda (2001).

Let us now examine the second term in Eq. (17), $\mathcal{J}(\mu)$, which
contains the leading infrared renormalon. Plugging Eq. (13) into the series
inside the integral, one finds

$\sum_{n=0}^{\infty}\frac{\Gamma(n+1+b)}{\Gamma(1+b)n!}z^{n}=\frac{1}{(1-z)^{1+b}},$
(19)

where $z=2\beta_{0}t$. Here, following the convention in performing a Borel
sum, a power series is summed in its region of validity, $|z|<1$, and then
extended by analytic continuation. Thus,

$\mathcal{J}(\mu)=\frac{R_{0}}{2\beta_{0}}\mu\int_{0}^{\infty}dz\,\frac{e^{-z/(2\beta_{0}\alpha_{\text{g}}(\mu))}}{(1-z)^{1+b}}.$
(20)

A full definition must specify the path of integration in the complex $z$
plane and the branch cut for the multivalued function $(1-z)^{-(1+b)}$.

Our strategy is to divide the integral into two pieces:

$\mathcal{J}(\mu)=\mathcal{J}_{\text{MRS}}(\mu)+\delta m,$ (21)

where

$\displaystyle\mathcal{J}_{\text{MRS}}(\mu)$
$\displaystyle=\frac{R_{0}}{2\beta_{0}}\mu\int_{0}^{1}dz\,\frac{e^{-z/[2\beta_{0}\alpha_{\text{g}}(\mu)]}}{(1-z)^{1+b}},$
(22) $\displaystyle\delta m$
$\displaystyle=\frac{R_{0}}{2\beta_{0}}\mu\int_{1}^{\infty}dz\,\frac{e^{-z/[2\beta_{0}\alpha_{\text{g}}(\mu)]}}{(1-z)^{1+b}}.$
(23)

Because the integration in $\mathcal{J}_{\text{MRS}}(\mu)$ lies within the
unit disk, it does not depend on the contour. In both pieces, the integral is
defined first for $\mathop{\rm Re}b<0$ and by analytic continuation elsewhere.
As discussed below, the first term, $\mathcal{J}_{\text{MRS}}$, is
unambiguous; the full ambiguity of the leading renormalon is contained in
$\delta m$.

In $\delta m$, the explicit $\mu$ dependence and implicit dependence via
$\alpha_{\text{g}}(\mu)$ cancel. This feature can be shown, and the branch cut
moved into a prefactor, by setting $z=1+2\beta_{0}\alpha_{\text{g}}(\mu)x$,
yielding

$\displaystyle\delta m$
$\displaystyle=(-1)^{1+b}\frac{R_{0}}{2\beta_{0}}\mu\frac{e^{-1/[2\beta_{0}\alpha_{\text{g}}(\mu)]}}{[2\beta_{0}\alpha_{\text{g}}(\mu)]^{b}}\int_{0}^{\infty}dx\,x^{-b-1}e^{-x}$
$\displaystyle=-(-1)^{b}\frac{R_{0}}{2^{1+b}\beta_{0}}\Gamma(-b)\Lambda_{\text{g}},$
(24)

where $(-1)^{b}=e^{\pm i\pi b}$ if the contour is set slightly above ($+$) or
below ($-$) the positive real axis. Remarkably, the conventional Bardeen  _et
al._ (1978); Celmaster and Gonsalves (1979); Chetyrkin  _et al._ (2000)
constant of integration introduced when solving Eq. (10) appears:

$\Lambda_{\text{g}}=\mu\frac{e^{-1/[2\beta_{0}\alpha_{\text{g}}(\mu)]}}{[\beta_{0}\alpha_{\text{g}}(\mu)]^{b}}.$
(25)

This constant is independent of $\mu$ to all orders, so, as announced, $\delta
m$ does not depend on $\mu$. Moreover, owing to the choice $b_{1}=0$ in Eq.
(12), it also follows, to all orders Bardeen  _et al._ (1978); Celmaster and
Gonsalves (1979), that

$\Lambda_{\text{g}}=\Lambda_{\overline{\rm MS}}.$ (26)

Finally, the overall normalization of the renormalon, $R_{0}$, is the same in
both schemes Beneke (1999), and indeed in any scheme such that $b_{1}=0$.
Bazavov  _et al._ (2012); *Bazavov:2014soa

At this point, we can present the main result of this section, namely a new
proposal for the mass to be used in setting up HQET and, thus, Eq. (6).
Reference Pineda (2001) introduces the _renormalon-subtracted_ (RS) mass,
which in our notation can be written

$\displaystyle m_{\text{RS}}(\nu_{f})\equiv
m_{\text{pole}}-\mathcal{J}(\nu_{f}),$ (27)

clearly removing the renormalon ambiguity.333For some practical applications
of the RS mass, see Ref. Pineda (2001); Bali and Pineda (2004); Brambilla  _et
al._ (2009, 2010); Bazavov  _et al._ (2012); Pineda and Segovia (2013); Ayala
_et al._ (2014); Berwein  _et al._ (2015). Our proposal is to subtract the
minimal piece of $\mathcal{J}(\mu)$ that removes the renormalon, namely
$\delta m$:

$\displaystyle m_{\text{MRS}}$ $\displaystyle\equiv m_{\text{pole}}-\delta m,$
(28)
$\displaystyle={\overline{m}}\left(1+\sum_{n=0}^{\infty}\left[r_{n}-R_{n}\right]\alpha_{\text{g}}^{n+1}({\overline{m}})\right)+\mathcal{J}_{\text{MRS}}({\overline{m}}).$
(29)

We call $m_{\text{MRS}}$ the _minimal renormalon-subtracted_ (MRS) mass,
because it removes only the part of $\mathcal{J}(\mu)$ that is ambiguous.
Equation (28) should be considered a formal, intuitive definition (both terms
are ambiguous), while Eq. (29) gives a rigorous (apart from subleading
renormalons and similar obstacles, overlooked in this analysis) and, thus,
operational definition. As the notation suggests, and as Eq. (24) shows, the
difference in using $m_{\text{MRS}}$ instead of $m_{\text{pole}}$ amounts to a
shift by an amount $\delta m\sim\Lambda_{\overline{\rm MS}}\ll m_{Q}$, which
does not disrupt the $1/m_{Q}$ expansion Neubert and Sachrajda (1995); Luke
_et al._ (1995).

To evaluate $m_{\text{MRS}}$, one can proceed as follows. First, one can carry
out the integration in Eq. (22), obtaining

$\mathcal{J}_{\text{MRS}}({\overline{m}})=\frac{R_{0}}{2\beta_{0}}{\overline{m}}e^{-1/[2\beta_{0}\alpha_{\text{g}}({\overline{m}})]}\Gamma(-b)\gamma^{\star}\left(-b,-[2\beta_{0}\alpha_{\text{g}}({\overline{m}})]^{-1}\right),$
(30)

where
$\gamma^{\star}(a,w)\equiv[w^{-a}/\Gamma(a)]\int_{0}^{w}dt\,t^{a-1}e^{-t}$ is
the limiting function of the incomplete gamma function Davis (1972). In
practice, the power series

$\Gamma(-b)\gamma^{\star}(-b,-y)=\sum_{n=0}^{\infty}\frac{y^{n}}{n!(n-b)}$
(31)

converges quickly. This series can be derived from Eq. (22) by changing
variables to $x=1-z$, expanding $e^{x/y}$ in powers of $x$, and integrating
term-by-term. Next, the sum in Eq. (29) is carried out to the order available,
taking care to be consistent with the scheme. The $r_{n}$ are used directly in
Eq. (29) and also in the $R_{n}$ via $R_{0}$ [Eq. (14)]. Because the sequence
corresponding to the leading infrared renormalon has been removed, the power
series for the MRS mass is better-behaved than the uninterpreted series, Eq.
(7). As shown in Sec. IV, Eqs. (77)–(79), the better behavior is realized in
practice.

This concludes the derivation and implementation of the MRS scheme for the
quark mass. We end this section with a few interesting and useful remarks.

It is interesting to examine the asymptotic expansions for
$\mathcal{J}_{\text{MRS}}$ and $m_{\text{MRS}}$, because they underscore a key
feature of asymptotic expansions: the same expansion can hold for more than
one quantity. As a function of $\alpha_{\text{g}}$,
$\mathcal{J}_{\text{MRS}}({\overline{m}})$ has a convergent expansion in
powers of $1/\alpha_{\text{g}}({\overline{m}})$, but it is not an analytic
function about $\alpha_{\text{g}}({\overline{m}})=0$. Its power expansion
leads to the asymptotic expansion

$\mathcal{J}_{\text{MRS}}({\overline{m}})\sim{\overline{m}}\sum_{n=0}^{\infty}R_{n}\alpha_{\text{g}}^{n+1}({\overline{m}}),\quad\text{as
}\alpha_{\text{g}}({\overline{m}})\to 0.$ (32)

The ambiguous function $\mathcal{J}({\overline{m}})$, given in Eq. (18), has
the same asymptotic expansion as $\mathcal{J}_{\text{MRS}}({\overline{m}})$,
because its ambiguous component $\delta m$ is exponentially small in
$1/\alpha_{\text{g}}({\overline{m}})$ and, thus, not visible to conventional
asymptotics. (See Ref. Boyd (1999) for a discussion on exponentially small
terms in asymptotic calculations.) Similarly, putting Eqs. (29) and (32)
together for the MRS mass, one finds

$m_{\text{MRS}}\sim{\overline{m}}+{\overline{m}}\sum_{n=0}^{\infty}r_{n}\alpha_{\text{g}}^{n+1}({\overline{m}}),\quad\text{as
}\alpha_{\text{g}}({\overline{m}})\to 0,$ (33)

which is the same series as the asymptotic expansion of the pole mass. Note
that with suitable changes to the $R_{n}$ and $r_{n}$, Eqs. (32) and (33) may
be extended to any renormalization scheme.

Another method to fix the ambiguity of the pole mass is the so-called
principal value (PV) prescription (Beneke, 1995; Lee, 2003; Bali and Pineda,
2004, e.g.), which removes the ambiguous imaginary part of $m_{\text{pole}}$:

$\displaystyle m_{\text{PV}}$ $\displaystyle\equiv\mathop{\rm
Re}m_{\text{pole}}$ (34) $\displaystyle=m_{\text{MRS}}-\cos(\pi
b)\frac{R_{0}}{2^{1+b}\beta_{0}}\Gamma(-b)\Lambda_{\text{g}}.$ (35)

Like $m_{\text{MRS}}$, $m_{\text{PV}}$ also contains no subtraction scale. The
PV mass has the advantage that it remains finite when $b$ tends to a positive
integer or zero. We prefer the MRS scheme, because the PV scheme entails an
arbitrary decision to replace the multivalued factor $(-1)^{b}$ with its real
part. In the MRS mass, on the other hand, all ambiguity is swept into $\delta
m$ and absorbed into the residual mass of HQET, as discussed below.

As just implied, the MRS mass is not ideal when $b$ tends to a non-negative
integer, in which case $\delta m$ is no longer of order $\Lambda_{\text{g}}$.
Suppose, for example, $b\to 0$. Then, we can use a variant MRS′ such that:

$\displaystyle m_{{\text{MRS}}^{\prime}}$
$\displaystyle=m_{\text{pole}}-\delta m^{\prime},$ (36) $\displaystyle\delta
m^{\prime}$ $\displaystyle=\delta
m-\frac{R_{0}}{2\beta_{0}b}\Lambda_{\text{g}},$ (37)
$\displaystyle\mathcal{J}_{{\text{MRS}}^{\prime}}({\overline{m}})$
$\displaystyle=\mathcal{J}_{\text{MRS}}({\overline{m}})+\frac{R_{0}}{2\beta_{0}b}\Lambda_{\text{g}},$
(38)

which introduces a term chosen to cancel the $n=0$ term in Eq. (31) as $b\to
0$. Other variants could remove the singular term in Eq. (31) as $b$ tends to
any positive integer.

With the MRS mass for heavy quarks, we can now revisit the discussion in Sec.
I on resolving the problem of ambiguities in HQET. We reconstruct the HQET
Lagrangian using the MRS mass for the heavy quark and setting the
corresponding residual mass444Our residual mass $\delta m$ is minus the
residual mass introduced in Ref. Falk  _et al._ (1992). to the ambiguous term
$\delta m$. Then, $M_{0}$ in Eq. (4) reads

$\displaystyle M_{0}$ $\displaystyle=\left[m_{\text{MRS}}+\delta
m+{\overline{\Lambda}}+\frac{\mu_{\pi}^{2}}{2m_{\text{MRS}}}+\text{O}\left(m_{\text{MRS}}^{-2}\right)\right]_{\text{chiral-
limit}},$
$\displaystyle=\left[m_{\text{MRS}}+{\overline{\Lambda}}_{\text{MRS}}+\frac{\mu_{\pi}^{2}}{2m_{\text{MRS}}}+\text{O}\left(m_{\text{MRS}}^{-2}\right)\right]_{\text{chiral-
limit}},$ (39)

where ${\overline{\Lambda}}_{\text{MRS}}={\overline{\Lambda}}+\delta m$ is a
well-defined quantity, because it is associated with two well-defined
quantities, $M_{0}$ and $m_{\text{MRS}}$. We do not attribute the MRS scheme
to $\mu_{\pi}^{2}$ in Eq. (39), because, as defined here, the MRS scheme does
not subtract subleading renormalons in the pole mass, such as the one canceled
by $\mu_{\pi}^{2}$. Renormalon-free definitions of $\overline{\Lambda}$ in
other schemes, such as PV and MRS′, follow by analogy with Eq. (39).

Finally, it is worth recording the relations between the RS and MRS
schemes:555Equation (40) provides an alternative way to calculate the RS mass,
similar to that used in Ref. Bali and Pineda (2004).

$\displaystyle m_{\text{RS}}(\nu_{f})$
$\displaystyle=m_{\text{MRS}}-\mathcal{J}_{\text{MRS}}(\nu_{f}),$ (40)
$\displaystyle{\overline{\Lambda}}_{\text{RS}}(\nu_{f})$
$\displaystyle={\overline{\Lambda}}_{\text{MRS}}+\mathcal{J}_{\text{MRS}}(\nu_{f}),$
(41)

which follow immediately from Eqs. (27) and (28), with
$\mathcal{J}_{\text{MRS}}(\nu_{f})$ evaluated as described in and below Eq.
(30).

##  III Heavy-light meson masses in staggered $\chi$PT

In this section, we discuss the corrections to heavy-light meson masses from
the “pion cloud,” using heavy-meson rooted all-staggered chiral perturbation
theory (HMrAS$\chi$PT), which was developed to compute these corrections for
heavy-light decay constants Bernard and Komijani (2013). HMrAS$\chi$PT
incorporates aspects of staggered fermions into standard HM$\chi$PT Burdman
and Donoghue (1992); Yan  _et al._ (1992); Wise (1992); Boyd and Grinstein
(1995), in particular the extra pseudo-Goldstone bosons arising from the four
fermion species, or tastes, that emerge from each staggered fermion field in
the continuum limit. For the sea quarks, HMrAS$\chi$PT also accounts for
taking the fourth root of the staggered determinant, which is done to simulate
one flavor per staggered fermion field.

The light pseudoscalar sector of the leading order HMrAS$\chi$PT Lagrangian
includes terms to describe the breaking of SU(4) taste symmetry at order
$a^{2}$ Lee and Sharpe (1999). For example, a light pseudoscalar meson with
flavors $x$ and $y$ and mesonic taste $\Xi$ has, at the tree level (Aubin and
Bernard, 2003, Eq. (18)),

$m^{2}_{{xy}_{\Xi}}=B_{0}(m_{x}+m_{y})+a^{2}\Delta_{\Xi}$ (42)

with $\Delta_{P}=0$, where $P$ is the so-called pseudoscalar taste. (The full
set of 16 mesonic tastes is given the labels $\\{I,V,T,A,P\\}$, corresponding
to multiplets of $(1,4,6,4,1)$, respectively Kronfeld (2007).) In contrast,
the heavy-light sector of the leading-order HMrAS$\chi$PT Lagrangian—and also
the one-loop heavy-light meson propagator—are taste invariant (Bernard and
Komijani, 2013, Eq. (27) and (151)). Therefore, all taste violations in the
heavy-light meson masses at next-to-leading order (NLO) come from the NLO
terms in the HMrAS$\chi$PT Lagrangian, treated at the tree level; see, for
instance, Eqs. (125) and (127) and Tables I and II in Ref. Bernard and
Komijani (2013). Because these terms can be added straightforwardly, we drop
the taste index of heavy-light mesons, for convenience, and we neglect the
tree-level lattice corrections to their masses, which can be restored at the
end, as required.

The leading corrections to heavy-light meson masses come from two sources: a
one-loop contribution from the self energy of the heavy-light meson, and the
tree-level contributions from the parts of the HMrAS$\chi$PT Lagrangian that
explicitly break the chiral symmetry of the light quarks and the spin-flavor
symmetry of the heavy quarks. The corresponding terms in the effective
Lagrangian are

$\displaystyle\mathcal{L}_{0}$ $\displaystyle=M_{0}\mathop{\rm
Tr}\left(\overline{H}H\right),$ (43) $\displaystyle\mathcal{L}_{1}$
$\displaystyle=-i\mathop{\rm
Tr}\left(\overline{H}Hv\negmedspace\cdot\negmedspace{\raise
6.45831pt\hbox{$\leftarrow$}\kern-10.00002pt}{D}\right)+g_{\pi}\mathop{\rm
Tr}\left(\overline{H}H\gamma^{\mu}\gamma_{5}\mathbb{A}_{\mu}\right),$ (44)
$\displaystyle\mathcal{L}_{m}$ $\displaystyle=2\lambda_{1}B_{0}\mathop{\rm
Tr}\left(\overline{H}H\mathcal{M}^{+}\right)+2\lambda^{\prime}_{1}B_{0}\mathop{\rm
Tr}\left(\overline{H}H\right)\mathop{\rm
Tr}\left(\mathcal{M}^{+}\right)-\frac{\lambda_{2}}{4m_{Q}}\mathop{\rm
Tr}\left(\overline{H}\sigma^{\mu\nu}H\sigma_{\mu\nu}\right),$ (45)

where $M_{0}$ from Eq. (4) is included for convenience.666One can choose the
generators of heavy-quark flavor symmetry to be compatible with a nontrivial
heavy-quark mass matrix $m_{Q}$ (Kronfeld, 2000, Sec. III). Also, the
$\mu_{\pi}^{2}$ term could, if desired, be moved from $\mathcal{L}_{0}$ to
$\mathcal{L}_{m}$. Here $H$ is a heavy-meson field for pseudoscalar and vector
mesons,
$\mathcal{M}^{+}={\textstyle\frac{1}{2}}[\sigma\mathcal{M}\sigma+\sigma^{\dagger}\mathcal{M}\sigma^{\dagger}]$
with $\mathcal{M}=\mathop{\rm diag}(m_{u},m_{d},m_{s})$,
$\sigma=\exp(i\Phi/2f)$, $\sigma_{\mu\nu}=(i/2)[\gamma_{\mu},\gamma_{\nu}]$,
$v^{\mu}$ is the four-velocity of HQET,
$D_{\mu}=\partial_{\mu}-i\mathbb{V}_{\mu}$,
$\mathbb{V}_{\mu}=\frac{i}{2}(\sigma^{\dagger}\partial_{\mu}\sigma+\sigma\partial_{\mu}\sigma^{\dagger})$,
$\mathbb{A}_{\mu}=\frac{i}{2}(\sigma^{\dagger}\partial_{\mu}\sigma-\sigma\partial_{\mu}\sigma^{\dagger})$,
and $g_{\pi}$ is the $H$-$H^{*}$-$\pi$ coupling. The matrix field $\Phi$
contains the $16n^{2}-1$ pseudo-Goldstone bosons for $n$ staggered light
flavors, and $f$ is the decay constant of these bosons (in the chiral limit).
As discussed elsewhere Aubin and Bernard (2003), $\Phi$ also contains the
flavor-taste singlet $\eta^{\prime}_{I}$, which is removed by taking
$m_{\eta^{\prime}_{I}}$ to infinity, relative to the pseudo-Goldstone boson
masses, once these contributions have been isolated. For full details of the
notation, see Refs. Aubin and Bernard (2003, 2006); Bernard and Komijani
(2013). In Eq. (45), we follow the normalization of Ref. Arndt and Lin (2004)
for $\lambda_{1}$ and $\lambda^{\prime}_{1}$ and of Ref. Stewart (1998) for
$\lambda_{2}$, but with opposite sign so that our $\lambda_{2}>0$. Further, in
Eq. (45), we use the power counting adopted in Ref. Stewart (1998); $m_{q}\sim
1/m_{Q}$ with $m_{q}$ being a generic light quark mass.

The heavy-light mesons carry both light flavor and taste indices. As mentioned
above, the taste index is unnecessary at this order, so we write $B_{x}$ and
$B^{*}_{x}$ for pseudoscalar and vector mesons, respectively, with light
valence flavor $x$. Their masses are then

$\displaystyle M_{B_{x}}$
$\displaystyle=M_{0}+2\lambda_{1}B_{0}m_{x}+2\lambda^{\prime}_{1}B_{0}\mathop{\rm
Tr}\mathcal{M}-\frac{3\lambda_{2}}{2m_{Q}}+\delta M_{B_{x}},$ (46a)
$\displaystyle M_{B^{*}_{x}}$
$\displaystyle=M_{0}+2\lambda_{1}B_{0}m_{x}+2\lambda^{\prime}_{1}B_{0}\mathop{\rm
Tr}\mathcal{M}+\frac{\lambda_{2}}{2m_{Q}}+\delta M_{B^{*}_{x}},$ (46b)

where the three terms after $M_{0}$ are the tree-level corrections, which are
immediate from Eq. (45), and $\delta{}M_{B^{(*)}_{x}}$ is the one-loop
correction from the pion cloud.

To obtain $\delta M_{B^{(*)}_{x}}$, one must evaluate the one-loop self energy
on shell. Although the flavor and hyperfine splittings in the heavy-meson
propagator are formally subleading in $m_{q}$ and $1/m_{Q}$, respectively, it
is quantitatively sensible to include them. In practice, both the flavor
splittings (_e.g._ , $M_{D_{s}}-M_{D}\approx 99$ MeV) and hyperfine splittings
(_e.g._ , $M_{D^{*}}-M_{D}\approx 144$ MeV) are not much different than the
pion mass, so they influence the behavior of the nonanalytic terms in the one-
loop contribution.

The self energy in HMrAS$\chi$PT has already been computed in Ref. Bernard and
Komijani (2013), which in turn relies on Ref. Aubin and Bernard (2006). There
the self energy is used to obtain the wave-function renormalization needed to
compute the one-loop corrections to the heavy-light decay constants. (Some
steps are not spelled out in Refs. Bernard and Komijani (2013); Aubin and
Bernard (2006), so we rederived the self energy from scratch.) Let us first
consider the case in which the three light sea quarks have distinct masses,
_i.e._ , $1+1+1$. We find

$\displaystyle\delta M_{B_{x}}\Bigr{|}_{1+1+1}$
$\displaystyle=-\frac{3g_{\pi}^{2}}{16\pi^{2}f^{2}}\Biggl{\\{}\frac{1}{16}\sum_{\mathscr{S},\Xi}K_{1}(m_{\mathscr{S}x_{\Xi}},\Delta^{*}+\delta_{\mathscr{S}x})$
(47)
$\displaystyle\qquad+\frac{1}{3}\sum_{j\in\mathcal{M}_{I}^{(3,x)}}\frac{\partial}{\partial
m^{2}_{X_{I}}}\left[R^{[3,3]}_{j}(\mathcal{M}_{I}^{(3,x)};\mu^{(3)}_{I})K_{1}(m_{j},\Delta^{*})\right]$
$\displaystyle\qquad+\biggl{(}a^{2}\delta^{\prime}_{V}\sum_{j\in\mathcal{M}_{V}^{(4,x)}}\frac{\partial}{\partial
m^{2}_{X_{V}}}\left[R^{[4,3]}_{j}(\mathcal{M}_{V}^{(4,x)};\mu^{(3)}_{V})K_{1}(m_{j},\Delta^{*})\right]+[V\to
A]\biggr{)}\Biggr{\\}},$ $\displaystyle\delta M_{B^{*}_{x}}\Bigr{|}_{1+1+1}$
$\displaystyle=-\frac{g_{\pi}^{2}}{16\pi^{2}f^{2}}\Biggl{\\{}\frac{1}{16}\sum_{\mathscr{S},\Xi}\Bigl{(}K_{1}(m_{\mathscr{S}x_{\Xi}},-\Delta^{*}+\delta_{\mathscr{S}x})+2K_{1}(m_{\mathscr{S}x_{\Xi}},\delta_{\mathscr{S}x})\Bigr{)}$
(48)
$\displaystyle\qquad+\frac{1}{3}\sum_{j\in\mathcal{M}_{I}^{(3,x)}}\frac{\partial}{\partial
m^{2}_{X_{I}}}\left[R^{[3,3]}_{j}(\mathcal{M}_{I}^{(3,x)};\mu^{(3)}_{I})K_{1}^{*}(m_{j},-\Delta^{*})\right]$
$\displaystyle\qquad+\biggl{(}a^{2}\delta^{\prime}_{V}\sum_{j\in\mathcal{M}_{V}^{(4,x)}}\frac{\partial}{\partial
m^{2}_{X_{V}}}\left[R^{[4,3]}_{j}(\mathcal{M}_{V}^{(4,x)};\mu^{(3)}_{V})K_{1}^{*}(m_{j},-\Delta^{*})\right]+[V\to
A]\biggr{)}\Biggr{\\}},$

where $\mathscr{S}$ runs over the three light sea flavors ($u$, $d$, $s$);
$\Xi$ runs over 16 meson tastes; $m_{\mathscr{S}x_{\Xi}}$ and $m_{X_{\Xi}}$
are masses of mesons with taste $\Xi$ and, respectively, flavor
$\bar{\mathscr{S}}x$ and $\bar{x}x$; and $\delta^{\prime}_{V}$
($\delta^{\prime}_{A}$) is the (axial) vector hairpin low-energy constant.

The notation for the loop functions is

$\displaystyle K_{1}(m,\Delta)$ $\displaystyle=\Delta J_{1}(m,\Delta),$ (49)
$\displaystyle K_{1}^{*}(m,\Delta)$
$\displaystyle=K_{1}(m,\Delta)+2K_{1}(m,0)=K_{1}(m,\Delta)+\frac{4\pi}{3}m^{3},$
(50) $\displaystyle J_{1}(m,\Delta)$
$\displaystyle=\left(\frac{2}{3}\Delta^{2}-m^{2}\right)\ln\left(\frac{m^{2}}{\Lambda_{\chi}^{2}}\right)+\frac{4}{3}\left(\Delta^{2}-m^{2}\right)F\left(\frac{m}{\Delta}\right)-\frac{10}{9}\Delta^{2}+\frac{4}{3}m^{2},$
(51) $\displaystyle F(1/x)$
$\displaystyle=\begin{cases}\displaystyle-\frac{\sqrt{1-x^{2}}}{x_{\phantom{g}}}\left(\frac{\pi}{2}-\tan^{-1}\frac{x}{\sqrt{1-x^{2}}}\right),&\text{if
$|x|\leq 1$},\\\
\displaystyle\frac{\sqrt{x^{2}-1}}{x}\ln\left(x+\sqrt{x^{2}-1}\right),&\text{if
$|x|\geq 1$}.\end{cases}$ (52)

In Eq. (51), $\Lambda_{\chi}$ is the renormalization scale of $\chi$PT. The
residue symbols originating from the hairpin diagrams take two sets of masses
as arguments:

$R_{j}^{[n,k]}(\\{m\\};\\{\mu\\})=\frac{\prod_{i=1}^{k}(\mu^{2}_{i}-m^{2}_{j})}{\prod_{r=1|r\neq
j}^{n}(m^{2}_{r}-m^{2}_{j})}.$ (53)

The sets needed in Eqs. (47) and (48) are

$\displaystyle\mu^{(3)}_{\Xi}$
$\displaystyle=\\{m_{U_{\Xi}},m_{D_{\Xi}},m_{S_{\Xi}}\\},\text{ for
}\Xi\in\\{I,V,A\\},$ (54) $\displaystyle\mathcal{M}^{(3,x)}_{I}$
$\displaystyle=\\{m_{X_{I}},m_{\pi^{0}_{I}},m_{\eta_{I}}\\},$ (55)
$\displaystyle\mathcal{M}^{(4,x)}_{\Xi}$
$\displaystyle=\\{m_{X_{\Xi}},m_{\pi^{0}_{\Xi}},m_{\eta_{\Xi}},m_{\eta^{\prime}_{\Xi}}\\},\text{
for }\Xi\in\\{V,A\\}.$ (56)

Here, $U_{\Xi}$, $D_{\Xi}$, and $S_{\Xi}$ denote $\bar{u}u_{\Xi}$,
$\bar{d}d_{\Xi}$, and $\bar{s}s_{\Xi}$ mesons, respectively, while
$\pi^{0}_{\Xi}$, $\eta_{\Xi}$, and $\eta^{\prime}_{\Xi}$ denote the
eigenstates after mixing among neutral mesons is taken into account. Only the
flavor-taste singlet $\eta^{\prime}_{I}$ receives a large mass related to the
axial anomaly and, as mentioned above, has been removed.

Finally, we have the hyperfine and flavor splittings.
$\Delta^{*}=2\lambda_{2}/m_{Q}$ is the lowest-order hyperfine splitting, and
$\delta_{\mathscr{S}x}$ is the lowest-order flavor splitting between heavy-
light mesons with light quarks of flavor $\mathscr{S}$ and $x$.
$\delta_{\mathscr{S}x}$ can be written in terms of pion masses as

$\delta_{\mathscr{S}x}=2\lambda_{1}B_{0}(m_{\mathscr{S}}-m_{x})\approx\lambda_{1}(m^{2}_{\mathscr{SS}_{P}}-m^{2}_{X_{P}}),$
(57)

where $P$ again denotes pseudoscalar taste.

In lattice-QCD simulations, the up and down sea-quark masses are often taken
equal, which is then denoted as a ($2+1$)-flavor sea. Throughout this paper,
we write $m_{l}\equiv\frac{1}{2}(m_{u}+m_{d})$. For the $2+1$ case, we find

$\displaystyle\delta M_{B_{x}}\Bigr{|}_{2+1}$
$\displaystyle=-\frac{3g_{\pi}^{2}}{16\pi^{2}f^{2}}\Biggl{\\{}\frac{1}{16}\sum_{\mathscr{S},\Xi}K_{1}(m_{\mathscr{S}x_{\Xi}},\Delta^{*}+\delta_{\mathscr{S}x})$
(58)
$\displaystyle\qquad+\frac{1}{3}\sum_{j\in\mathcal{M}_{I}^{(2,x)}}\frac{\partial}{\partial
m^{2}_{X_{I}}}\left[R^{[2,2]}_{j}(\mathcal{M}_{I}^{(2,x)};\mu^{(2)}_{I})K_{1}(m_{j},\Delta^{*})\right]$
$\displaystyle\qquad+\biggl{(}a^{2}\delta^{\prime}_{V}\sum_{j\in\hat{\mathcal{M}}_{V}^{(3,x)}}\frac{\partial}{\partial
m^{2}_{X_{V}}}\left[R^{[3,2]}_{j}(\hat{\mathcal{M}}_{V}^{(3,x)};\mu^{(2)}_{V})K_{1}(m_{j},\Delta^{*})\right]+[V\to
A]\biggr{)}\Biggr{\\}},$ $\displaystyle\delta M_{B^{*}_{x}}\Bigr{|}_{2+1}$
$\displaystyle=-\frac{g_{\pi}^{2}}{16\pi^{2}f^{2}}\Biggl{\\{}\frac{1}{16}\sum_{\mathscr{S},\Xi}\Bigl{(}K_{1}(m_{\mathscr{S}x_{\Xi}},-\Delta^{*}+\delta_{\mathscr{S}x})+2K_{1}(m_{\mathscr{S}x_{\Xi}},\delta_{\mathscr{S}x})\Bigr{)}$
(59)
$\displaystyle\qquad+\frac{1}{3}\sum_{j\in\mathcal{M}_{I}^{(2,x)}}\frac{\partial}{\partial
m^{2}_{X_{I}}}\left[R^{[2,2]}_{j}(\mathcal{M}_{I}^{(2,x)};\mu^{(2)}_{I})K_{1}^{*}(m_{j},-\Delta^{*})\right]$
$\displaystyle\qquad+\biggl{(}a^{2}\delta^{\prime}_{V}\sum_{j\in\hat{\mathcal{M}}_{V}^{(3,x)}}\frac{\partial}{\partial
m^{2}_{X_{V}}}\left[R^{[3,2]}_{j}(\hat{\mathcal{M}}_{V}^{(3,x)};\mu^{(2)}_{V})K_{1}^{*}(m_{j},-\Delta^{*})\right]+[V\to
A]\biggr{)}\Biggr{\\}},$

with sets

$\displaystyle\mu^{(2)}_{\Xi}$
$\displaystyle=\\{m_{L_{\Xi}},m_{S_{\Xi}}\\},\text{ for }\Xi\in\\{I,V,A\\},$
(60) $\displaystyle\mathcal{M}^{(2,x)}_{I}$
$\displaystyle=\\{m_{X_{I}},m_{\eta_{I}}\\},$ (61)
$\displaystyle\hat{\mathcal{M}}^{(3,x)}_{\Xi}$
$\displaystyle=\\{m_{X_{\Xi}},m_{\eta_{\Xi}},m_{\eta^{\prime}_{\Xi}}\\},\text{
for }\Xi\in\\{V,A\\},$ (62)

where $L_{\Xi}$, analogously to Eq. (54), denotes an $\bar{l}l_{\Xi}$ meson.
In the continuum limit of full QCD (where the light valence quark is one of
the sea quarks), the nonanalytic terms in these expressions agree with the
corresponding results, Eqs. (A6), (A7), (A8) and (A9), in HM$\chi$PT Mehen and
Springer (2005); there are unimportant differences in the analytic terms,
stemming from a different choice of counterterms.

As they stand, the expressions for $\delta M_{B^{(*)}_{x}}$ are not completely
independent of the chiral scale $\Lambda_{\chi}$, but this dependence can be
absorbed into counterterms. Most of these are the same as in continuum
HM$\chi$PT, and they are formally of higher order in the $\chi$PT expansion.
An interesting example is suppressed by $a^{2}$ and can be handled via

$\delta M_{B^{(*)}_{x}}\to\delta
M_{B^{(*)}_{x}}+\lambda_{a^{2}}a^{2}\frac{3g_{\pi}^{2}}{16\pi^{2}f^{2}}\left[\bar{\Delta}\sum_{\mathscr{S}}\delta_{\mathscr{S}x}+d_{B^{(*)}}\Delta^{*}\left(3\bar{\Delta}-{\textstyle\frac{1}{3}}\Delta_{I}+\delta^{\prime}_{V}+\delta^{\prime}_{A}\right)\right],$
(63)

where $d_{B^{(*)}}=1(-{\textstyle\frac{1}{3}})$, and the nearly777These
splittings are empirically nearly flavor independent (see, _e.g._ , Ref.
(Bazavov  _et al._ , 2013, Fig. 6)), which can be understood because the
flavor dependence vanishes at leading order in staggered $\chi$PT Aubin and
Bernard (2003). $\mathscr{S}x$-independent quantities

$\displaystyle a^{2}\Delta_{I}$
$\displaystyle=m_{\mathscr{S}x_{I}}^{2}-m_{\mathscr{S}x_{P}}^{2},$ (64)
$\displaystyle a^{2}\bar{\Delta}$
$\displaystyle=\frac{1}{16}\sum_{\Xi}\left(m_{\mathscr{S}x_{\Xi}}^{2}-m_{\mathscr{S}x_{P}}^{2}\right)$
(65)

stem from the taste splittings in the logarithmic terms in $K_{1}$; cf. Eqs.
(49) and (51). The remaining $\Lambda_{\chi}$ dependence can be absorbed into
$\lambda_{1}$, $\lambda^{\prime}_{1}$ and higher order counterterms in
continuum HM$\chi$PT.

Finite-volume effects can be incorporated into $\delta
M_{\raisebox{5.0pt}{$\scriptstyle B^{(*)}_{x}$}}$ by substituting
$K_{1}(m,\Delta)$ with its finite-volume counterpart. Equation (92) of Ref.
Aubin and Bernard (2007) gives the finite-volume version of $J_{1}(m,\Delta)$.
With Eq. (49), we then find, for a spatial volume $L^{3}$,

$K_{1}(m,\Delta)\to K_{1,\text{FV}}(m,\Delta)=K_{1}(m,\Delta)+\delta
K_{1}(m,\Delta,L),$ (66)

where

$\delta
K_{1}(m,\Delta,L)=-\frac{m^{2}}{3}\Delta\delta_{1}(mL)-16\pi^{2}\frac{m^{2}-\Delta^{2}}{3}J_{\text{FV}}(m,\Delta,L),$
(67)

with $\delta_{1}(mL)$ and $J_{\text{FV}}(m,\Delta,L)$ given, respectively, in
Eqs. (84) and (86) of Ref. Aubin and Bernard (2007).

##  IV Synthesis

In this section, we combine our results from the previous two sections,
deriving the HQET and $\chi$PT descriptions of heavy-light meson masses, into
a practical form. In addition, we comment on the effect of the nonzero charm
quark mass in the relation between the $\overline{\rm MS}$ and MRS masses.

Let us first fix notation for quark masses associated with currently available
lattice-QCD ensembles with $2+1+1$ flavors of quarks Bazavov  _et al._ (2013).
We use $m^{\prime}_{l}$, $m^{\prime}_{s}$ and $m^{\prime}_{c}$ to denote the
simulation masses of the light, strange and charm quarks, respectively. We
also use $m_{l}=\frac{1}{2}(m_{u}+m_{d})$, $m_{s}$ and $m_{c}$ to denote the
tuned values of the corresponding quarks. As in Sec. III, we use $B_{x}^{(*)}$
to denote a generic all-staggered heavy-light pseudoscalar (vector) meson with
light valence quark $x$ and heavy valence antiquark $\bar{h}$.

Putting Eqs. (4), (5), and (46) together, using the MRS mass for the heavy
quark, and spelling out the dependence on light quark masses and lattice
spacing, we have

$\displaystyle
M_{B_{x}^{(*)}}(m_{x};\\{m^{\prime}_{l},m^{\prime}_{l},m^{\prime}_{s}\\};a)$
$\displaystyle=\left[m_{h,{\text{MRS}}}+{\overline{\Lambda}}_{\text{MRS}}+\frac{\mu_{\pi}^{2}-d_{B^{(*)}}\mu_{G}^{2}(m_{h})}{2m_{h,{\text{MRS}}}}\right]_{\text{chiral-
limit}}$ (68)
$\displaystyle+2\lambda_{1}B_{0}m_{x}+2\lambda^{\prime}_{1}B_{0}(2m^{\prime}_{l}+m^{\prime}_{s})+\delta
M_{B_{x}^{(*)}}(m_{x};\\{m^{\prime}_{l},m^{\prime}_{l},m^{\prime}_{s}\\};a),$

where the last term is given by the suitable expression from Eqs. (47), (48),
(58), or (59). In Eq. (68), the HQET matrix elements
${\overline{\Lambda}}_{\text{MRS}}$, $\mu_{\pi}^{2}$ and $\mu_{G}(m_{h})$
appear in the SU(3) chiral limit, which is not interesting in practice.
Therefore, we rewrite Eq. (68) as

$\displaystyle
M_{B_{x}^{(*)}}(m_{x};\\{m^{\prime}_{l},m^{\prime}_{l},m^{\prime}_{s}\\};a)$
$\displaystyle=m_{h,{\text{MRS}}}+{\overline{\Lambda}}_{\text{MRS}}+\frac{\mu_{\pi}^{2}-d_{B^{(*)}}\mu_{G}^{2}(m_{h})}{2m_{h,{\text{MRS}}}}+2\lambda_{1}B_{0}m_{x}$
(69)
$\displaystyle+2\lambda^{\prime}_{1}B_{0}(2m^{\prime}_{l}+m^{\prime}_{s})+\delta
M_{B_{x}^{(*)}}(m_{x};\\{m^{\prime}_{l},m^{\prime}_{l},m^{\prime}_{s}\\};a)-\mathcal{C}^{(*)},$

where

$\mathcal{C}^{(*)}=2\lambda_{1}B_{0}m_{q}+2\lambda^{\prime}_{1}B_{0}(2m_{l}+m_{s})+\delta
M_{B_{q}^{(*)}}(m_{q};\\{m_{l},m_{l},m_{s}\\};0)$ (70)

is the continuum, physical quark-mass limit of the $\chi$PT corrections. By
construction, then, the formula reduces to

$M_{B_{q}^{(*)}}(m_{q};\\{m_{l},m_{l},m_{s}\\};0)=m_{h,{\text{MRS}}}+{\overline{\Lambda}}_{\text{MRS}}+\frac{\mu_{\pi}^{2}-d_{B^{(*)}}\mu_{G}^{2}(m_{h})}{2m_{h,{\text{MRS}}}}$
(71)

for $a=0$, tuned sea-quark masses, and $m_{q}$ takes one of the values
${\textstyle\frac{1}{2}}(m_{u}+m_{d})$, $m_{u}$, $m_{d}$, or $m_{s}$,
depending on the light flavor of interest. The resulting values of
${\overline{\Lambda}}_{\text{MRS}}$, $\mu_{\pi}^{2}$ and $\mu_{G}^{2}(m_{h})$
are now appropriate for calculations within HQET.

From HQET power counting, it may seem that chiral corrections may be moved
arbitrarily between the HQET matrix elements and the heavy quark mass. It is
crucial, however, that the heavy-quark mass be defined so that it does not
depend on the hadron of which it is a constituent. More formally, the heavy-
quark mass should be a singlet under chiral symmetry and appear in all hadron
masses in a way consistent with HQET. Therefore, it is necessary to use
HM$\chi$PT to separate all low-energy light-quark mass effects from
$M_{B^{(*)}}$ to obtain $M_{0}$, and then to use HQET to separate $M_{0}$ into
$m_{h,{\text{MRS}}}$, ${\overline{\Lambda}}_{\text{MRS}}$, and
$\mu_{\pi}^{2}/2m_{Q}$. In this way, all chiral corrections are absorbed into
${\overline{\Lambda}}_{\text{MRS}}$, $\mu_{\pi}^{2}$, $\mu_{G}^{2}$, etc. in
Eq. (71).

One can use Eq. (3) to express $\mu_{G}^{2}(m_{h})$ in terms of
$C_{\text{cm}}(m_{h})$ and $\tilde{\mu}_{G}^{2}$, which are the Wilson
coefficient and matrix element of the chromomagnetic operator in HQET,
respectively. Equivalently, one can express it as

$\mu_{G}^{2}(m_{h})=\frac{C_{\text{cm}}(m_{h})}{C_{\text{cm}}(m_{b})}\mu_{G}^{2}(m_{b}).$
(72)

When lattice-QCD data are not available for the $B^{*}$ mass (as in Ref.
Bazavov  _et al._ (2017a)), one can use the experimental value of the
hyperfine splitting, $M_{B^{*}}-M_{B}$, as a prior on the value of
$\mu_{G}^{2}(m_{b})$. Because $\mu_{G}^{2}(m_{b})$ provides the first term in
the $1/m_{h}$ expansion for the physical quantity $M_{B^{*}}-M_{B}$, any
renormalon ambiguity must be suppressed by at least another factor of
$1/m_{h}$. This is not the case for $\mu_{\pi}^{2}$, because it is not closely
related to any physical quantity, as we discuss briefly below.

A further practical step is to exploit the identity888The left-hand side of
Eq. (73) holds with any mass-independent renormalization scheme, and the
right-hand side relies on the remnant chiral symmetry of staggered fermions.

$\frac{m_{r,{\overline{\rm MS}}}(\mu)}{m_{h,{\overline{\rm
MS}}}(\mu)}=\frac{am_{r,0}}{am_{h,0}}+\text{O}(a^{2}),$ (73)

where $am_{0}$ is the bare staggered-fermion mass in lattice units, and the
subscript $r$ denotes a reference mass. In Eq. (69), we trivially rewrite

$m_{h,{\text{MRS}}}=\frac{m_{r,{\overline{\rm MS}}}(\mu)}{m_{h,{\overline{\rm
MS}}}(\mu)}\frac{am_{h,0}}{am_{r,0}}m_{h,{\text{MRS}}}=m_{r,{\overline{\rm
MS}}}(\mu)\frac{{\overline{m}}_{h}}{m_{h,{\overline{\rm
MS}}}(\mu)}\frac{m_{h,{\text{MRS}}}}{{\overline{m}}_{h}}\frac{am_{h,0}}{am_{r,0}}.$
(74)

The factors in the last expression are, respectively, a convenient fit
parameter for lattice-QCD applications, a factor depending on $\alpha_{s}$ and
the mass anomalous dimension to run between ${\overline{m}}_{h}$ and a fixed
scale $\mu$, an expression given in Sec. II, and a variable that is an input
to lattice QCD. The third factor depends on ${\overline{m}}_{h}$ only through
$\alpha_{\text{g}}({\overline{m}}_{h})$, as is easily seen in Eqs. (29) and
(30). A different renormalization scheme, besides $\overline{\rm MS}$, would
amount to substituting the appropriate coefficients for the $r_{n}$ in Eq.
(29), and another renormalon subtraction, besides MRS, would amount to
substituting the appropriate prescription for the $\mathcal{J}$ function.

Let us now address an important detail in the perturbative series relating the
$\overline{\rm MS}$ and MRS definitions of the renormalized quark mass. The
relation between the ${\overline{\rm MS}}$ mass and the pole mass (and hence
the MRS mass) is known through order $\alpha_{s}^{4}$, but the
order-$\alpha_{s}^{4}$ term is known only for massless sea quarks. The effect
of nonzero sea-quark masses is available through order $\alpha_{s}^{3}$ Gray
_et al._ (1990); Bekavac  _et al._ (2007). For the up, down, and strange
quarks, the effects of nonzero masses of quarks can be neglected in
perturbative QCD. For the charmed sea quark, however, the nonzero mass does
affect the heavy-quark mass at the current precision. A bit more generally,
let us consider QCD with $n_{l}$ massless quarks, one charm quark of mass
$0<{\overline{m}}_{c}<{\overline{m}}_{h}$, and $n_{h}$ heavy quarks of mass
${\overline{m}}_{h}$. The relation between the pole and ${\overline{\rm MS}}$
mass of the $n_{h}$ heavy quarks is usually first reported with
$n_{l}+1+n_{h}$ active quarks; see, for example, Ref. Marquard  _et al._
(2016). The result, however, is usually expressed and used with $n_{l}+1$
flavors of active quarks, absorbing the heavy-quark loops into the running of
$\alpha_{s}$. Although the gauge coupling is converted to run with $n_{l}+1$
active flavors, the heavy-quark mass customarily remains in the
$n_{l}+1+n_{h}$ scheme.

When studying the relation between the $\overline{\rm MS}$ and RS masses, Ref.
Ayala  _et al._ (2014) found a better-behaved expansion if the charm-quark
contribution is also decoupled from the renormalon-subtracted series and added
back separately. This can be explained by the fact that the charm-quark mass
cuts off the infrared momenta that generate factorially growing contributions
Ball  _et al._ (1995). In this vein, one can rewrite Eq. (29) as

$m_{h,{\text{MRS}}}={\overline{m}}_{h}\left(1+\sum_{n=0}^{\infty}\big{[}r_{n}^{(n_{l})}-R_{n}^{(n_{l})}\big{]}\alpha_{s}^{n+1}({\overline{m}}_{h};n_{l})\right)+\mathcal{J}_{\text{MRS}}^{(n_{l})}({\overline{m}}_{h})+\Delta
m_{(c)},$ (75)

where $\alpha_{s}({\overline{m}}_{h};n_{l})$ is the value of the gauge
coupling with $n_{l}$ active quarks and the charm quark decoupled at scale
$\mu={\overline{m}}_{c}$; the superscript $(n_{l})$ in $r_{n}^{(n_{l})}$,
$R_{n}^{(n_{l})}$ and $\mathcal{J}_{\text{MRS}}^{(n_{l})}({\overline{m}}_{h})$
indicates that they correspond to the theory with $n_{l}$ active quarks; and
$\Delta m_{(c)}$ includes the contribution of the charmed quark as well as
effects that are generated from matching to the $n_{l}$-flavor $\alpha_{s}$.
More details are provided at the end of this section. Note that, even though
the charmed quark is removed from the perturbative series, the $\overline{\rm
MS}$ mass ${\overline{m}}_{h}$ in Eq. (75) remains in the
($n_{l}+1+n_{h}$)-flavor scheme. In the calculation of the MRS mass, we then
set

$R_{0}^{(3)}=0.535\pm 0.010,$ (76)

which is the overall normalization of the leading infrared renormalon in the
pole mass with three massless active quarks Komijani (2017). The value agrees
with other estimates in the literature: $0.563(26)$ Ayala  _et al._ (2014),
$0.537~{}\text{(no~{}error~{}quoted)}$ Beneke  _et al._ (2016), and
$0.526(12)$ Hoang  _et al._ (2017).

To see how these manipulations work in practice, let us now investigate
numerically the relation between the MRS and ${\overline{\rm MS}}$ masses for
$n_{l}=3$, $R_{0}^{(3)}=0.535$, and $n_{h}=0$ (which is appropriate for
2+1+1-flavor lattice QCD). Exploiting four-loop calculations in Ref. Marquard
_et al._ (2016), we have

$\displaystyle r_{n}^{(3)}$ $\displaystyle=(0.4244,\ 1.0351,\ 3.6932,\
17.4358,\ldots),$ (77) $\displaystyle R_{n}^{(3)}$ $\displaystyle=(0.5350,\
1.0691,\ 3.5966,\ 17.4195,\ldots),$ (78)

for $n=0,1,2,3,\ldots$. Their differences

$r_{n}^{(3)}-R_{n}^{(3)}=(-0.1106,\ -0.0340,\ 0.0966,\ 0.0162,\ldots)$ (79)

are much smaller; consequently, the series in powers of $\alpha_{s}$ in Eq.
(75) can be considered to be a well-behaved series through order
$\alpha_{s}^{4}$. This good behavior indicates that the subleading infrared
renormalon is not an immediate problem in computing the MRS mass. In
principle, this renormalon is canceled by a corresponding renormalon in
$\mu_{\pi}^{2}$, so Eq. (79) also suggests that $\mu_{\pi}^{2}$ will be
numerically stable in the MRS scheme. This outcome is not surprising, because
the effect of this renormalon is known to be small in Lorentz-invariant
regularization schemes Neubert (1997).

Let us return to Eq. (75) and give the explicit expression for $\Delta
m_{(c)}$, which is denoted $\delta m_{c}$ in Ref. Ayala  _et al._ (2014). The
key equation is Eq. (2.15) in Ref. Ayala  _et al._ (2014); in our notation,

$\Delta m_{(c)}=\left(\delta m_{(c,+)}^{(1)}+\delta
m_{(c,\text{dec})}^{(1)}\right)\left[\frac{\alpha_{s}({\overline{m}}_{h};n_{l})}{\pi}\right]^{2}+\left(\delta
m_{(c,+)}^{(2)}+\delta
m_{(c,\text{dec})}^{(2)}\right)\left[\frac{\alpha_{s}({\overline{m}}_{h};n_{l})}{\pi}\right]^{3},$
(80)

where the coefficients $\delta m_{(c,+)}^{(1)}$, $\delta
m_{(c,\text{dec})}^{(1)}$, and $\delta m_{(c,\text{dec})}^{(2)}$ are given in
Eqs. (2.6), (2.16), and (2.17) of Ref. Ayala  _et al._ (2014), and the
coefficient $\delta m_{(c,+)}^{(2)}$ is presented in Ref. Bekavac  _et al._
(2007). Because the quarks that are heavier than the charm quark are (still)
quenched in lattice-QCD simulations, we cannot directly use the expressions
provided in Ref. Ayala  _et al._ (2014). As above, we use $n_{h}$ to denote
the number of dynamical heavy quarks with mass ${\overline{m}}_{h}$; then, for
$n_{l}=3$ and arbitrary $n_{h}$, one has

$\displaystyle\delta m_{(c,\text{dec})}^{(2)}=$
$\displaystyle\,{\overline{m}}_{h}\left[-\frac{115981}{11664}-\frac{29\pi^{2}}{24}+\frac{61\pi^{4}}{1944}-\frac{11}{81}\pi^{2}\ln
2+\frac{2}{81}\pi^{2}\ln^{2}2+\frac{\ln^{4}2}{81}\right.$ (81)
$\displaystyle+\frac{8}{27}\,\text{Li}_{4}({\textstyle\frac{1}{2}})-\frac{511}{216}\zeta(3)-\ln
z^{2}\left(\frac{377}{96}+\frac{\pi^{2}}{27}\ln
2-\frac{1}{18}\zeta(3)-\frac{1}{27}\ln z^{2}\right)$
$\displaystyle\left.\vphantom{\frac{\ln^{4}2}{81}}+\frac{n_{h}}{11664}\left(5917-864\zeta(3)-468\pi^{2}-432\pi^{2}\ln
z^{2}+3861\ln z^{2}\right)\right]-\frac{1}{3}\delta m_{(c,+)}^{(1)}\ln z^{2},$

where $z={\overline{m}}^{\prime}_{c}/{\overline{m}}_{h}$, and
${\overline{m}}^{\prime}_{c}$ is the simulation charmed-quark mass. To obtain
this expression from Ref. Ayala  _et al._ (2014), results in Refs. Chetyrkin
_et al._ (1997); Melnikov and Ritbergen (2000) are needed to restore the
$n_{h}$ dependence. Equation (81) reduces to Eq. (2.17) of Ref. Ayala  _et
al._ (2014) for $n_{h}=1$. For analysis of lattice-QCD data with 2+1+1 flavors
of sea quarks, one should set $n_{h}=0$. Further, to avoid the lengthy
expressions for $\delta m_{(c,+)}^{(1)}$ and especially $\delta
m_{(c,+)}^{(2)}$, we use the approximate expressions provided in Eqs. (3.2)
and (3.4) of Ref. Bekavac  _et al._ (2007). We find numerically, for $n_{l}=3$
and $n_{h}=0$,

$\displaystyle\delta m_{(c,+)}^{(1)}$
$\displaystyle={\overline{m}}^{\prime}_{c}\left(1.596-0.6285z+0.1777z^{2}\right),$
(82) $\displaystyle\delta m_{(c,\text{dec})}^{(1)}$
$\displaystyle=-{\overline{m}}_{h}\left[1.0414+0.4444\ln z\right],$ (83)
$\displaystyle\delta m_{(c,+)}^{(2)}$
$\displaystyle={\overline{m}}^{\prime}_{c}\left[19.987+2.824z+1.288z^{2}-\left(13.644+2.788z-0.0343z^{2}\right)\ln
z\right],$ (84) $\displaystyle\delta m_{(c,\text{dec})}^{(2)}$
$\displaystyle=-{\overline{m}}_{h}\left[22.312+\left(8.227+1.064z-0.419z^{2}+0.118z^{3}-0.148\ln
z\right)\ln z\right].$ (85)

The relations for $\delta m_{(c,+)}^{(2)}$ and $\delta
m_{(c,\text{dec})}^{(2)}$ change slightly for $n_{h}=1$:

$\displaystyle\delta m_{(c,+)}^{(2)}$
$\displaystyle={\overline{m}}^{\prime}_{c}\left[19.987+3.028z+1.141z^{2}-\left(13.644+2.788z-0.0343z^{2}\right)\ln
z\right],$ (86) $\displaystyle\delta m_{(c,\text{dec})}^{(2)}$
$\displaystyle=-{\overline{m}}_{h}\left[22.29+\left(8.296+1.064z-0.419z^{2}+0.118z^{3}-0.148\ln
z\right)\ln z\right],$ (87)

which may prove useful at a later date.

##  V Summary and Outlook

As we show in a companion paper Bazavov  _et al._ (2017a), the theoretical
developments in this paper provide a new way to determine the quark masses as
well as HQET matrix elements such as ${\overline{\Lambda}}$ and
$\mu_{\pi}^{2}$.999The data used in Ref. Bazavov  _et al._ (2017a) lack the
vector-meson masses. An ideal analysis would include them, in order to obtain
information on $\mu_{G}^{2}$ as well. With extremely precise lattice-QCD data
for the heavy-light pseudoscalar-meson correlation functions Bazavov  _et al._
(2017b), Ref. Bazavov  _et al._ (2017a) yields some of the most precise
determinations of all quark masses (other than the top quark). We find there
that the fits of the lattice-QCD data are much more convenient with the MRS
mass than with alternatives. The key is that we do not have to choose a
renormalization scale, such as $\nu_{f}$ for the RS scheme, while preserving
the natural structure of HQET. We have tried the HQET fits with the RS mass
and found that we had to introduce three scales in all,
$\nu_{f}<\mu<{\overline{m}}_{h}$, with $\mu$ being used for $\alpha_{s}$
Komijani  _et al._ (2016). In addition, the end results for
${\overline{m}}_{c}$ and ${\overline{m}}_{b}$ depend on $\nu_{f}$ more than
one would like. These experiences prompted us to reexamine the idea of
renormalon subtraction, leading to the MRS scheme. With the available data, we
also found more stable fits to the $m_{h}$ dependence when using a wide range
of data for $m_{x}$ and $m^{\prime}_{l}$. This experience prompted us to
derive the HMrAS$\chi$PT expressions in Sec. III.

The idea to fit the HQET formula, Eq. (2), to lattice-QCD data to obtain
$\bar{\Lambda}$ and $\mu_{\pi}^{2}$ has been tried before. The first attempt
Kronfeld and Simone (2000) used lattice perturbation theory to convert the
bare lattice mass to the pole mass. In addition to suffering from the
renormalon ambiguity, the conversion was (and usually will be) limited to
order $\alpha_{s}$. Closer in spirit is the recent work Gambino  _et al._
(2017) using nonperturbative renormalization of the lattice quark mass,
together with the kinetic scheme Uraltsev (1997, 2002) in Eq. (2). The
conversion from the kinetic to the $\overline{\rm MS}$ scheme is known to
(only) two loops Gambino (2011). With the MRS scheme, one can use the four-
loop relation between the $\overline{\rm MS}$ mass and the pole mass Marquard
_et al._ (2016) and, in the future, any higher-order terms that may be
obtained.

For lattice QCD with staggered sea and valence quarks, we augment the MRS
scheme with $\chi$PT for staggered fermions. Similar procedures can be applied
with other lattice fermion formulations. For example, $\chi$PT for Wilson
fermions addresses chiral-symmetry breaking from the dimension-five Wilson
term in a way consistent with that of the mass term.

One can, in principle, extend the definition of the MRS mass by subtracting
sequences of higher and higher renormalons from the perturbative series and
collecting their unambiguous contributions into a generalization of the
function $\mathcal{J}_{\text{MRS}}({\overline{m}})$. The key is to split the
Borel integral, as in Eqs. (22) and (23), into unambiguous and ambiguous
pieces, and sweep the ambiguity into $\delta m$. With some imagination and
work, it may be possible to extend the MRS idea beyond renormalon
singularities to nonperturbative saddle points, _i.e._ , instantons. One may
then speculate that the remaining series in powers of $\alpha_{s}$ could
acquire a nonzero radius of convergence. Clearly, such a step will require
generalizations of Ref. Komijani (2017), or other approaches such as direct
investigation of divergences in the Borel plane, and, thus, remains a
conjecture at this stage.

How might this work? As explained in Ref. Berry (1991), “an asymptotic series
… is a compact encoding of a function, and its divergence should be regarded
not as a deficiency but as a source of information about the [underlying]
function.” The crucial point here is that, equipped with the overall
normalization Komijani (2017), we could decode the divergence due to the
leading infrared renormalon in the pole mass. Then we could recover the
contribution of this renormalon in the form of the nonperturbative function
$\mathcal{J}_{\text{MRS}}({\overline{m}})$, which has a convergent expansion
in powers of the reciprocal of $\alpha_{s}({\overline{m}})$, rather than in
powers of $\alpha_{s}({\overline{m}})$ itself.

The MRS scheme should also be useful in other phenomenological settings—the
extra scale(s) of alternative short-distance masses require more choices,
making uncertainty estimation difficult. For example, the rate of inclusive
semileptonic $B$ decay is proportional to $m^{5}_{b,{\text{pole}}}$,
multiplied by a perturbative series that cancels the renormalon ambiguities
Neubert and Sachrajda (1995); Beneke  _et al._ (1994). Using instead the MRS
scheme for mass and corrections, one would have the product of two unambiguous
expressions.

In theoretical treatments of quarkonium, the heavy-quark potential contains a
renormalon ambiguity canceled by $2m_{\text{pole}}$ Pineda (1998); Hoang  _et
al._ (1999). An immediate alternative is the MRS scheme, moving the minimal
renormalon from the mass to the potential. A first application could be in the
extraction of $\alpha_{s}$ from the comparison of the static energy computed
in perturbation theory with lattice QCD, as implemented with a different
renormalon subtraction scheme in Ref. Bazavov  _et al._ (2012);
*Bazavov:2014soa. It would be also interesting to use the MRS mass to examine
the quarkonium spectrum, for example to determine the quark mass, either along
the lines of Refs. Brambilla  _et al._ (2001); Kiyo and Sumino (2014); Mateu
and Ortega (2017) or using lattice QCD.

Another application may be the gluelump mass Campbell  _et al._ (1986), which
arises in the short-range treatment of hybrids Brambilla  _et al._ (2000);
Bali and Pineda (2004); Berwein  _et al._ (2015). As an adjoint-representation
analog of ${\overline{\Lambda}}$, the gluelump mass either contains a
renormalon (as in Ref. (Bali and Pineda, 2004, Eq. (6))) or depends on an
artificial scale (as in some other subtraction schemes). The MRS approach
could be extended to provide a physical definition, making possible a direct
lattice-QCD calculation along the lines of our companion determination of
quark masses and ${\overline{\Lambda}}_{\text{MRS}}$ Bazavov  _et al._
(2017a). Similar quantities can be formulated in baryons and tetraquark states
with one or more heavy quarks. Although all examples mentioned here are
closely related to the renormalons in the pole mass and-or heavy-quark
potential, renormalons appear in other contexts in QCD, and the MRS scheme may
be useful in those areas too.

###### Acknowledgements.

J.K. and A.S.K. thank Martin Beneke, Claude Bernard, and Doug Toussaint for
conversations that influenced this paper. This work was supported in part by
the German Excellence Initiative and the European Union Seventh Framework
Program under grant agreement No. 291763 as well as the European Union’s Marie
Curie COFUND program (J.K., A.S.K.). The work of N.B. and A.V. was supported
by the DFG cluster of excellence “Origin and structure of the universe”.
Fermilab is operated by Fermi Research Alliance, LLC, under Contract No. DE-
AC02-07CH11359 with the United States Department of Energy, Office of Science,
Office of High Energy Physics. The United States Government retains and the
publisher, by accepting the article for publication, acknowledges that the
United States Government retains a non-exclusive, paid-up, irrevocable, world-
wide license to publish or reproduce the published form of this manuscript, or
allow others to do so, for United States Government purposes.

## References

  * Bazavov  _et al._ (2017a) A. Bazavov _et al._ (Fermilab Lattice, MILC, TUMQCD), to appear (2017a).
  * Burdman and Donoghue (1992) G. Burdman and J. F. Donoghue, Phys. Lett. B280, 287 (1992).
  * Yan  _et al._ (1992) T.-M. Yan, H.-Y. Cheng, C.-Y. Cheung, G.-L. Lin, Y. C. Lin, and H.-L. Yu, Phys. Rev. D46, 1148 (1992), (E) Phys. Rev. D55, 5851 (1997).
  * Wise (1992) M. B. Wise, Phys. Rev. D45, R2188 (1992).
  * Boyd and Grinstein (1995) C. G. Boyd and B. Grinstein, Nucl. Phys. B442, 205 (1995), arXiv:hep-ph/9402340 .
  * Falk and Neubert (1993) A. F. Falk and M. Neubert, Phys. Rev. D47, 2965 (1993), arXiv:hep-ph/9209268 .
  * Kronfeld (2000) A. S. Kronfeld, Phys. Rev. D62, 014505 (2000), arXiv:hep-lat/0002008 .
  * Komijani (2017) J. Komijani, JHEP 08, 062 (2017), arXiv:1701.00347 [hep-ph] .
  * Luke and Manohar (1992) M. E. Luke and A. V. Manohar, Phys. Lett. B286, 348 (1992), arXiv:hep-ph/9205228 .
  * Brambilla  _et al._ (2003) N. Brambilla, D. Gromes, and A. Vairo, Phys. Lett. B576, 314 (2003), arXiv:hep-ph/0306107 .
  * Heinonen  _et al._ (2012) J. Heinonen, R. J. Hill, and M. P. Solon, Phys. Rev. D86, 094020 (2012), arXiv:1208.0601 [hep-ph] .
  * Grozin  _et al._ (2008) A. G. Grozin, P. Marquard, J. H. Piclum, and M. Steinhauser, Nucl. Phys. B789, 277 (2008), arXiv:0707.1388 [hep-ph] .
  * Bernard and Komijani (2013) C. Bernard and J. Komijani, Phys. Rev. D88, 094017 (2013), arXiv:1309.4533 [hep-lat] .
  * Kronfeld (1998) A. S. Kronfeld, Phys. Rev. D58, 051501 (1998), arXiv:hep-ph/9805215 .
  * Breckenridge  _et al._ (1995) J. C. Breckenridge, M. J. Lavelle, and T. G. Steele, Z. Phys. C65, 155 (1995), arXiv:hep-th/9407028 .
  * Bigi  _et al._ (1994) I. I. Bigi, M. A. Shifman, N. G. Uraltsev, and A. I. Vainshtein, Phys. Rev. D50, 2234 (1994), arXiv:hep-ph/9402360 .
  * Beneke and Braun (1994) M. Beneke and V. M. Braun, Nucl. Phys. B426, 301 (1994), arXiv:hep-ph/9402364 .
  * Luke  _et al._ (1995) M. E. Luke, A. V. Manohar, and M. J. Savage, Phys. Rev. D51, 4924 (1995), arXiv:hep-ph/9407407 .
  * Uraltsev (1997) N. Uraltsev, Nucl. Phys. B491, 303 (1997), arXiv:hep-ph/9610425 .
  * Uraltsev (2002) N. Uraltsev, Phys. Lett. B545, 337 (2002), arXiv:hep-ph/0111166 .
  * Pineda (2001) A. Pineda, JHEP 06, 022 (2001), arXiv:hep-ph/0105008 .
  * Hoang  _et al._ (2008) A. H. Hoang, A. Jain, I. Scimemi, and I. W. Stewart, Phys. Rev. Lett. 101, 151602 (2008), arXiv:0803.4214 [hep-ph] .
  * Hoang (2000) A. H. Hoang, Phys. Rev. D61, 034005 (2000), arXiv:hep-ph/9905550 .
  * Beneke (1998) M. Beneke, Phys. Lett. B434, 115 (1998), arXiv:hep-ph/9804241 .
  * Brambilla  _et al._ (1999) N. Brambilla, A. Pineda, J. Soto, and A. Vairo, Phys. Rev. D60, 091502 (1999), arXiv:hep-ph/9903355 .
  * Falk  _et al._ (1992) A. F. Falk, M. Neubert, and M. E. Luke, Nucl. Phys. B388, 363 (1992), arXiv:hep-ph/9204229 .
  * Neubert and Sachrajda (1995) M. Neubert and C. T. Sachrajda, Nucl. Phys. B438, 235 (1995), arXiv:hep-ph/9407394 .
  * Kawai and Takei (2005) T. Kawai and Y. Takei,  _Algebraic analysis of singular perturbation theory_, Vol. 227 (American Mathematical Soc., 2005).
  * Beneke (1995) M. Beneke, Phys. Lett. B344, 341 (1995), arXiv:hep-ph/9408380 .
  * Lee (2003) T. Lee, JHEP 10, 044 (2003), arXiv:hep-ph/0304185 .
  * Lee (2015) T. Lee, Phys. Lett. B742, 327 (2015), arXiv:1502.02698 [hep-ph] .
  * Hoang  _et al._ (2017) A. H. Hoang, A. Jain, C. Lepenik, V. Mateu, M. Preisser, I. Scimemi, and I. W. Stewart, (2017), arXiv:1704.01580 [hep-ph] .
  * Pineda (2017) A. Pineda, (2017), arXiv:1704.05095 [hep-ph] .
  * Pineda (1998) A. Pineda,  _Heavy Quarkonium and Nonrelativistic Effective Field Theories_ , Ph.D. thesis, Unversitat de Barcelona (1998).
  * Hoang  _et al._ (1999) A. H. Hoang, M. C. Smith, T. Stelzer, and S. Willenbrock, Phys. Rev. D59, 114014 (1999), arXiv:hep-ph/9804227 .
  * Brown  _et al._ (1992) L. S. Brown, L. G. Yaffe, and C.-X. Zhai, Phys. Rev. D46, 4712 (1992), arXiv:hep-ph/9205213 .
  * Lee (1997) T. Lee, Phys. Rev. D56, 1091 (1997), arXiv:hep-th/9611010 .
  * Bardeen  _et al._ (1978) W. A. Bardeen, A. J. Buras, D. W. Duke, and T. Muta, Phys. Rev. D18, 3998 (1978).
  * Celmaster and Gonsalves (1979) W. Celmaster and R. J. Gonsalves, Phys. Rev. D20, 1420 (1979).
  * Chetyrkin  _et al._ (2000) K. G. Chetyrkin, J. H. Kühn, and M. Steinhauser, Comput. Phys. Commun. 133, 43 (2000), arXiv:hep-ph/0004189 .
  * Beneke (1999) M. Beneke, Phys. Rept. 317, 1 (1999), arXiv:hep-ph/9807443 .
  * Brambilla  _et al._ (2009) N. Brambilla, A. Vairo, X. Garcia i Tormo, and J. Soto, Phys. Rev. D80, 034016 (2009), arXiv:0906.1390 [hep-ph] .
  * Brambilla  _et al._ (2010) N. Brambilla, X. Garcia i Tormo, J. Soto, and A. Vairo, Phys. Rev. Lett. 105, 212001 (2010), (E) Phys. Rev. Lett. 108, 269903(2012), arXiv:1006.2066 [hep-ph] .
  * Bazavov  _et al._ (2012) A. Bazavov, N. Brambilla, X. Garcia i Tormo, P. Petreczky, J. Soto, and A. Vairo, Phys. Rev. D86, 114031 (2012), arXiv:1205.6155 [hep-ph] .
  * Bazavov  _et al._ (2014) A. Bazavov, N. Brambilla, X. Garcia i Tormo, P. Petreczky, J. Soto, and A. Vairo, Phys. Rev. D90, 074038 (2014), arXiv:1407.8437 [hep-ph] .
  * Pineda and Segovia (2013) A. Pineda and J. Segovia, Phys. Rev. D87, 074024 (2013), arXiv:1302.3528 [hep-ph] .
  * Ayala  _et al._ (2014) C. Ayala, G. Cvetič, and A. Pineda, JHEP 09, 045 (2014), arXiv:1407.2128 [hep-ph] .
  * Berwein  _et al._ (2015) M. Berwein, N. Brambilla, J. Tarrús Castellá, and A. Vairo, Phys. Rev. D92, 114019 (2015), arXiv:1510.04299 [hep-ph] .
  * Bali and Pineda (2004) G. S. Bali and A. Pineda, Phys. Rev. D69, 094001 (2004), arXiv:hep-ph/0310130 .
  * Davis (1972) P. J. Davis, “Gamma function and related functions,” in  _Handbook of Mathematical Functions_ , edited by M. Abramowitz and I. A. Stegun (Dover, New York, 1972) Chap. 6, pp. 253–294.
  * Boyd (1999) J. P. Boyd, Acta Applicandae Mathematica 56, 1 (1999).
  * Lee and Sharpe (1999) W.-J. Lee and S. R. Sharpe, Phys. Rev. D60, 114503 (1999), arXiv:hep-lat/9905023 .
  * Aubin and Bernard (2003) C. Aubin and C. Bernard, Phys. Rev. D68, 034014 (2003), arXiv:hep-lat/0304014 .
  * Kronfeld (2007) A. S. Kronfeld, PoS LATTICE2007, 016 (2007), arXiv:0711.0699 [hep-lat] .
  * Aubin and Bernard (2006) C. Aubin and C. Bernard, Phys. Rev. D73, 014515 (2006), arXiv:hep-lat/0510088 .
  * Arndt and Lin (2004) D. Arndt and C. J. D. Lin, Phys. Rev. D70, 014503 (2004), arXiv:hep-lat/0403012 .
  * Stewart (1998) I. W. Stewart, Nucl. Phys. B529, 62 (1998), arXiv:hep-ph/9803227 .
  * Mehen and Springer (2005) T. Mehen and R. P. Springer, Phys. Rev. D72, 034006 (2005), arXiv:hep-ph/0503134 .
  * Bazavov  _et al._ (2013) A. Bazavov _et al._ (MILC), Phys. Rev. D87, 054505 (2013), arXiv:1212.4768 [hep-lat] .
  * Aubin and Bernard (2007) C. Aubin and C. Bernard, Phys. Rev. D76, 014002 (2007), arXiv:0704.0795 [hep-lat] .
  * Gray  _et al._ (1990) N. Gray, D. J. Broadhurst, W. Grafe, and K. Schilcher, Z. Phys. C48, 673 (1990).
  * Bekavac  _et al._ (2007) S. Bekavac, A. Grozin, D. Seidel, and M. Steinhauser, JHEP 10, 006 (2007), arXiv:0708.1729 [hep-ph] .
  * Marquard  _et al._ (2016) P. Marquard, A. V. Smirnov, V. A. Smirnov, M. Steinhauser, and D. Wellmann, Phys. Rev. D94, 074025 (2016), arXiv:1606.06754 [hep-ph] .
  * Ball  _et al._ (1995) P. Ball, M. Beneke, and V. M. Braun, Nucl. Phys. B452, 563 (1995), arXiv:hep-ph/9502300 .
  * Beneke  _et al._ (2016) M. Beneke, P. Marquard, P. Nason, and M. Steinhauser, (2016), arXiv:1605.03609 [hep-ph] .
  * Neubert (1997) M. Neubert, Phys. Lett. B393, 110 (1997), arXiv:hep-ph/9610471 .
  * Chetyrkin  _et al._ (1997) K. G. Chetyrkin, B. A. Kniehl, and M. Steinhauser, Phys. Rev. Lett. 79, 2184 (1997), arXiv:hep-ph/9706430 .
  * Melnikov and Ritbergen (2000) K. Melnikov and T. v. Ritbergen, Phys. Lett. B482, 99 (2000), arXiv:hep-ph/9912391 .
  * Bazavov  _et al._ (2017b) A. Bazavov _et al._ (Fermilab Lattice, MILC), to appear (2017b).
  * Komijani  _et al._ (2016) J. Komijani _et al._ (Fermilab Lattice, MILC, TUMQCD), PoS LATTICE2016, 294 (2016), arXiv:1611.07411 [hep-lat] .
  * Kronfeld and Simone (2000) A. S. Kronfeld and J. N. Simone, Phys. Lett. B490, 228 (2000), (E) Phys. Lett. B495, 441 (2000), arXiv:hep-ph/0006345 .
  * Gambino  _et al._ (2017) P. Gambino, A. Melis, and S. Simula, Phys. Rev. D96, 014511 (2017), arXiv:1704.06105 [hep-lat] .
  * Gambino (2011) P. Gambino, JHEP 09, 055 (2011), arXiv:1107.3100 [hep-ph] .
  * Berry (1991) M. Berry, “Asymptotics, superasymptotics, hyperasymptotics…,” in  _Asymptotics beyond All Orders_, edited by H. Segur, S. Tanveer, and H. Levine (Plenum, New York, 1991) pp. 1–14.
  * Beneke  _et al._ (1994) M. Beneke, V. M. Braun, and V. I. Zakharov, Phys. Rev. Lett. 73, 3058 (1994), arXiv:hep-ph/9405304 .
  * Brambilla  _et al._ (2001) N. Brambilla, Y. Sumino, and A. Vairo, Phys. Lett. B513, 381 (2001), arXiv:hep-ph/0101305 .
  * Kiyo and Sumino (2014) Y. Kiyo and Y. Sumino, Phys. Lett. B730, 76 (2014), arXiv:1309.6571 [hep-ph] .
  * Mateu and Ortega (2017) V. Mateu and P. G. Ortega, (2017), arXiv:1711.05755 [hep-ph] .
  * Campbell  _et al._ (1986) N. A. Campbell, I. H. Jorysz, and C. Michael, Phys. Lett. 167B, 91 (1986).
  * Brambilla  _et al._ (2000) N. Brambilla, A. Pineda, J. Soto, and A. Vairo, Nucl. Phys. B566, 275 (2000), arXiv:hep-ph/9907240 .