# Chiral effective Lagrangian for heavy-light mesons from QCD: $1/m_{Q}$
correction

Qing-Sen Chen  Center for Theoretical Physics, College of Physics, Jilin
University, Changchun 130012, China  Hui-Feng Fu<EMAIL_ADDRESS>Center
for Theoretical Physics, College of Physics, Jilin University, Changchun
130012, China  Yong-Liang Ma<EMAIL_ADDRESS>School of Fundamental Physics
and Mathematical Sciences, Hangzhou Institute for Advanced Study, UCAS,
Hangzhou, 310024, China  International Center for Theoretical Physics Asia-
Pacific (ICTP-AP) (Beijing/Hangzhou), UCAS, Beijing 100190, China  Qing Wang
<EMAIL_ADDRESS>Department of Physics, Tsinghua University, Beijing
100084, China  Center for High Energy Physics, Tsinghua University, Beijing
100084, China

###### Abstract

As a successive work to [Phys.Rev.D 102 (2020), 034034], we derive the
$1/m_{Q}$ corrections to chiral effective Lagrangian for heavy-light mesons
from QCD under proper approximations. The low energy constants in the
effective Lagrangian are expressed in terms of the light quark self-energy and
heavy quark mass $m_{Q}$. Numerical results of the low energy constants with
$1/m_{Q}$ corrections are given. We find that the results of pion decay
constant and the masses of heavy-light mesons are improved coherently compared
to that obtained in the heavy quark limit.

##  I Introduction

Establishing the analytic relationships between the low energy constants
(LECs) of the chiral effective theory and QCD has great significance in hadron
physics. Such relationships fill up the gap between the effective theories and
the fundamental theory, and allow one to calculate the LECs from the QCD Green
functions. This type of researches has been carried out in Refs. Wang  _et
al._ (2000a); Yang  _et al._ (2002); Jiang  _et al._ (2010); Jiang and Wang
(2010); Jiang  _et al._ (2013, 2015) for the traditional chiral effective
theory and in Refs. Wang and Wang (2000); Wang  _et al._ (2000b); Ren  _et
al._ (2017) for some extensions to the chiral effective theory. Recently, the
chiral effective theory for heavy-light mesons were derived from QCD in the
heavy quark limit $m_{Q}\rightarrow\infty$ Chen  _et al._ (2020).

In Ref. Chen  _et al._ (2020), we focused on the heavy-light meson doublets
with the spin-parity of the light quark cloud $s_{l}^{P}=\frac{1}{2}^{-}$ (or
$J^{P}=(0^{-},1^{-})$ for the spin-parity of the heavy-light mesons) and
$s_{l}^{P}=\frac{1}{2}^{+}$ (or $J^{P}=(0^{+},1^{+})$) which are denoted as
$H$ and $G$ respectively. The LECs in the effective Lagrangian (including the
masses of the mesons and the coupling constants) were expressed in terms of
the light quark self-energy which can be calculated using Dyson-Schwinger
equations or lattice QCD. The resulted numerical values turned out to be
roughly consistent with the experimental data. The mass splitting of $H$ and
$G$ doublets emerges evidently from dynamical chiral symmetry breaking of QCD
which is realized by the solution to the quark gap equation.

Soon after the early works on the chiral effective theory for heavy-light
mesons in the heavy quark limit Wise (1992); Yan  _et al._ (1992); Nowak  _et
al._ (1993); Casalbuoni  _et al._ (1993), the effects due to finite heavy
quark mass began to be explored Cheng  _et al._ (1994); Balk  _et al._ (1994);
Di Bartolomeo  _et al._ (1995) and continue to be an important topic until now
Cheung and Hwang (2016); Alhakami (2020). Since our framework establishes the
relationships between the LECs and the fundamental theory QCD, it provides a
method to calculate the $1/m_{Q}$ corrections from QCD. We study in this work
the $1/m_{Q}$ corrections to the LECs of the chiral Lagrangian for the heavy-
light mesons in order to improve our previous results. It is found that the
$1/m_{Q}$ corrections are indeed helpful in improving the numerical results.
Especially, the tension between the pion decay constant $f_{\pi}$ and the mass
splitting, when calculated up to the $1/m_{Q}$ order, is released. We study
both the charmed mesons and the bottom mesons in this paper, for these two
meson sectors suffer from different values of $1/m_{Q}$ corrections, which
provides more detailed information on the validity of the heavy quark
expansion.

This paper is organized as follows. In Sec. II, we introduce the chiral
effective Lagrangian for the heavy-light mesons. In Sec. III, the LECs of the
chiral Lagrangian up to $1/m_{Q}$ order is derived from QCD, which turn to be
integrals of the dressing functions of the light quark propagator. The
numerical results based on these formula are given in Sec. IV, where the
relevant dressing functions of the quark propagator are obtained from the
quark Dyson-Schwinger equation as well as from lattice QCD. A summary is given
in Section V.

##  II chiral effective Lagrangian for heavy-light mesons

For convenience, we introduce the chiral effective Lagrangian for heavy-light
mesons in this section. The heavy-light meson doublets $H$ and $G$ can be
expressed as

$\displaystyle H$ $\displaystyle=$
$\displaystyle\frac{1+v/}{2}\left(P^{\ast;\mu}\gamma_{\mu}+iP\gamma_{5}\right),$
$\displaystyle G$ $\displaystyle=$
$\displaystyle\frac{1+v/}{2}\left(S^{\ast;\mu}\gamma_{\mu}\gamma_{5}+S\right),$
(1)

where $(P,P^{\ast;\mu})$ refer to $J^{P}=(0^{-},1^{-})$ states, and
$(S,S^{\ast;\mu})$ refer to $J^{P}=(0^{+},1^{+})$ states, respectively.
$v^{\mu}$ is the velocity of an on-shell heavy quark, i.e.
$p^{\mu}=m_{Q}v^{\mu}$ with $v^{2}=1$. As in Ref. Chen  _et al._ (2020), we
only consider two light flavors. Then the chiral effective Lagrangian for
heavy-light mesons (up to first order derivatives) can be written as Manohar
and Wise (2000); Nowak  _et al._ (2004)

$\displaystyle{\cal L}$ $\displaystyle=$ $\displaystyle{\cal L}_{H}+{\cal
L}_{G}+{\cal L}_{HG},$ (2)

where 111Here, since we are not interested in the mass splitting of the
hadrons in a heavy quark doublet, we shall not write done the effective
Lagrangian containing magnetic moment operator.

$\displaystyle{\cal L}_{H}$ $\displaystyle=$ $\displaystyle{}-i{\rm
Tr}\left(\bar{H}v\cdot\nabla H\right)-g_{H}{\rm
Tr}\left(H\gamma^{\mu}\gamma_{5}A_{\mu}\bar{H}\right)$
$\displaystyle{}+m_{H}{\rm Tr}\left(\bar{H}H\right),$ $\displaystyle{\cal
L}_{G}$ $\displaystyle=$ $\displaystyle{}-i{\rm Tr}\left(\bar{G}v\cdot\nabla
G\right)+g_{G}{\rm Tr}\left(G\gamma^{\mu}\gamma_{5}A_{\mu}\bar{G}\right)$
$\displaystyle{}+m_{G}{\rm Tr}\left(\bar{G}G\right),$ $\displaystyle{\cal
L}_{HG}$ $\displaystyle=$ $\displaystyle g_{HG}{\rm
Tr}\left(H\gamma^{\mu}\gamma_{5}A_{\mu}\bar{G}\right)+h.c.,$ (3)

where the covariant derivative $\nabla_{\mu}=\partial_{\mu}-iV_{\mu}$ with
$V_{\mu}=\frac{i}{2}\left(\Omega^{\dagger}\partial_{\mu}\Omega+\Omega\partial_{\mu}\Omega^{\dagger}\right)$,
and
$A_{\mu}=\frac{i}{2}\left(\Omega^{\dagger}\partial_{\mu}\Omega-\Omega\partial_{\mu}\Omega^{\dagger}\right)$.
The $\Omega$ field is related to the chiral field $U(x)=\exp(i\pi(x)/f_{\pi})$
through $U=\Omega^{2}$. The parameters $g_{H},g_{G},g_{HG},m_{H}$, and $m_{G}$
are the LECs of the Lagrangian. They are free ones at the level of effective
theory since they cannot be controlled by symmetry argument which the
effective theory relies on.

The states associated with $H$ and $G$ are called chiral partners Nowak  _et
al._ (1993); Bardeen and Hill (1994), and their mass splitting arises from the
chiral symmetry breaking. In the $D$ meson family, where $H$ is associated to
$(D,D^{*})$ and $G$ is associated to $(D^{*}_{0},D_{1})$, the spin-averaged
masses of the chiral partners are Zyla  _et al._ (2020)

$\displaystyle m_{H}$ $\displaystyle\simeq$ $\displaystyle 1970~{}{\rm
MeV},\;\;\;m_{G}\simeq 2400~{}{\rm MeV},$ (4)

which yields the mass splitting

$\displaystyle\Delta m$ $\displaystyle=$ $\displaystyle m_{G}-m_{H}=430~{}{\rm
MeV}.$ (5)

In the $B$ meson family, where $H$ is associated to $(B,B^{*})$ and $G$ is
associated to $(B^{*}_{0},B_{1})$, the spin-averaged masses of the chiral
partners are Zyla  _et al._ (2020)

$\displaystyle m_{H}$ $\displaystyle\simeq$ $\displaystyle 5302~{}{\rm
MeV},\;\;\;m_{G}\simeq 5726~{}{\rm MeV},$ (6)

and the mass splitting is

$\displaystyle\Delta m$ $\displaystyle=$ $\displaystyle m_{G}-m_{H}=424~{}{\rm
MeV}.$ (7)

##  III  The $1/m_{Q}$ corrections to chiral effective Lagrangian for heavy-
light mesons from QCD

To obtain the $1/m_{Q}$ corrections to the chiral effective Lagrangian for the
heavy-light mesons, we follow our previous analysis but retain only the
$1/m_{Q}$ order contributions in the heavy quark expansion. The QCD generating
functional can be written as

$\displaystyle Z[J]$ $\displaystyle=$
$\displaystyle\int\mathcal{D}q\mathcal{D}\bar{q}\mathcal{D}Q\mathcal{D}\bar{Q}\mathcal{D}G_{\mu}\varDelta_{F}(G_{\mu})\exp\left\\{i\int
d^{4}x\left[\bar{q}(i\not{D}+J)q+\bar{Q}(i\not{D}-m_{Q})Q-\frac{1}{4}G_{\mu\nu}^{i}G^{\mu\nu}_{i}-\dfrac{1}{2\xi}[F^{i}(G_{\mu})]^{2}\right]\right\\},$
(8)

where $q(x)$, $Q(x)$ and $G_{\mu}(x)$ are the light-quark, heavy-quark and
gluon fields, respectively. $J(x)$ is an external source for the composite
light quark fields. The masses of the light quarks are absorbed in the
external source $J(x)$ and would vanish in the chiral limit. The heavy quark
mass is denoted as $m_{Q}$. The indices $i,j,\cdots$ represent the color
indices in the adjoint representation.

Following the standard heavy quark effective theory (HQET) Manohar and Wise
(2000), we formulate the effective Lagrangian directly in terms of the
velocity-dependent fields $N_{v}(x)$ and $\mathcal{N}_{v}(x)$:

$\displaystyle N_{v}(x)$ $\displaystyle=$
$\displaystyle\frac{1+\not{v}}{2}e^{im_{Q}v\cdot x}Q(x),$
$\displaystyle\mathcal{N}_{v}(x)$ $\displaystyle=$
$\displaystyle\frac{1-\not{v}}{2}e^{im_{Q}v\cdot x}Q(x),$ (9)

where $N_{v}(x)$ is the large component of the heavy quark field and
$\mathcal{N}_{v}(x)$ is the component of the heavy quark field suppressed by
power of $1/m_{Q}$. The factor $e^{im_{Q}v\cdot x}$ rotates away the mass of
the heavy quark field. Integrating out $\mathcal{N}_{v}(x)$ at tree level
amounts to the replacement

$\displaystyle\mathcal{N}_{v}(x)=\frac{i\not{D}_{\perp}}{iv\cdot
D+2m_{Q}}N_{v},$ (10)

where $\not{D}_{\perp}=\not{D}-v\cdot D\not{v}$, is the perpendicular
component of $\not{D}$. Then the heavy quark part of the Lagrangian
$\mathcal{L}_{Q}=\bar{Q}(i\not{D}-m_{Q})Q$ becomes

$\begin{array}[]{ccl}\mathcal{L}_{Q}&=&\bar{N}_{v}\left(iv\cdot
D-\frac{\not{D}_{\perp}\not{D}_{\perp}}{iv\cdot D+2m_{Q}}\right)N_{v}\\\
&=&\bar{N}_{v}\left(iv\cdot\partial-\frac{(\not{\partial}_{\perp})^{2}}{2m_{Q}}\right)N_{v}+\mathcal{L}_{int}^{N_{v}G},\end{array}$
(11)

where
$(\not{\partial}_{\perp})^{2}=(\not{\partial}-v\cdot\partial\not{v})^{2}=\partial^{2}-(v\cdot\partial)^{2}$
and $\mathcal{L}_{int}^{N_{v}G}$ accounts for the interaction between the
heavy quark field $N_{v}$ and the gluon field. Up to the $1/m_{Q}$ order, the
generating functional becomes

$\displaystyle Z[J]$ $\displaystyle=$
$\displaystyle\int\mathcal{D}q\mathcal{D}\bar{q}\mathcal{D}N_{v}\mathcal{D}\bar{N}_{v}\exp\left\\{i\int
d^{4}x\left[\bar{q}(i\not{\partial}+J)q+\bar{N}_{v}\left(iv\cdot\partial-\frac{(\not{\partial}_{\perp})^{2}}{2m_{Q}}\right)N_{v}\right]\right\\}$
(12)
$\displaystyle{}~{}~{}\times\int\mathcal{D}G_{\mu}\varDelta_{F}(G_{\mu})\exp\left\\{i\int
d^{4}x\left[{}-\frac{1}{4}G_{\mu\nu}^{i}G^{\mu\nu}_{i}-\dfrac{1}{2\xi}[F^{i}(G_{\mu})]^{2}-g\mathcal{I}^{\mu}_{i}G^{i}_{\mu}+\mathcal{L}_{int}^{N_{v}G}\right]\right\\},$

where $\mathcal{I}_{i}^{\mu}=\bar{q}\frac{\lambda_{i}}{2}\gamma^{\mu}q$ is the
light quark current.

The chiral effective action for the heavy-light mesons can be obtained by
first integrating in the chiral field $U(x)$ and the heavy-light meson fields
$\sim(q\bar{N}_{v})$ to and then integrating out gluon fields and quark fields
from the generating functional (12). Taking the large $N_{c}$ limit and
keeping the leading order in the dynamical perturbation, we obtain the
effective action as

$\displaystyle S[U,\Pi_{2},\bar{\Pi}_{2}]$ $\displaystyle=$ $\displaystyle-
iN_{c}\mathrm{Tr}^{\prime}\ln\Bigg{[}(i\not{\partial}-\bar{\Sigma})I_{1}+J^{\prime}_{\Omega}$
(13)
$\displaystyle{}\qquad\qquad\quad\,+\left(iv\cdot\partial-\frac{(\not{\partial}_{\perp})^{2}}{2m_{Q}}\right)I_{4}$
$\displaystyle{}\qquad\qquad\quad\,-\Pi_{2}-\bar{\Pi}_{2}\Bigg{]},$

where $\Pi_{2}$ and $\bar{\Pi}_{2}$ are the heavy-light meson field and its
conjugate, respectively. $\bar{\Sigma}$ is the self-energy of the light quark
propagator. $J^{\prime}_{\Omega}$ is the chiral-rotated external source which
the chiral field is attached on. $I_{1}=\mathrm{diag}(1,1,0)$ and
$I_{4}=\mathrm{diag}(0,0,1)$ are the light- and heavy- projecting matrices in
the flavor space respectively. $\mathrm{Tr}^{\prime}$ represents a functional
trace over the flavor space, spinor space, and coordinate space. The details
of the derivation of Eq. (13) are given in Appendix A.

Since our purpose is to obtain the chiral effective Lagrangian for the chiral
partners given by Eq. (1), we only keep the $H$ and $G$ fields in $\Pi_{2}$:

$\displaystyle\Pi_{2}(x)$ $\displaystyle=$ $\displaystyle H^{q}(x)+G^{q}(x),$
(14)

where $q=1,2$ is light flavor indices. The chiral effective Lagrangian is
generated by expanding the action $S[U,\Pi_{2},\bar{\Pi}_{2}]$ with respect to
the fields $U$, $\Pi_{2}$ and $\bar{\Pi}_{2}$. The kinetic terms of the heavy-
light meson fields arise from

$\displaystyle S_{2}$ $\displaystyle\equiv$
$\displaystyle\frac{\delta^{2}S}{\delta\bar{\Pi}_{2}\delta\Pi_{2}}\Pi_{2}\bar{\Pi}_{2}$
$\displaystyle=$ $\displaystyle iN_{c}\int d^{4}x_{1}d^{4}x_{2}$
$\displaystyle{}\times\mathrm{Tr}\Bigg{[}\left(i\not{\partial}-\bar{\Sigma}\right)^{-1}\delta(x_{2}-x_{1})\bar{H}(x_{1})$
$\displaystyle\qquad\qquad\times\left(iv\cdot\partial-\frac{(\not{\partial}_{\perp})^{2}}{2m_{Q}}\right)^{-1}\delta(x_{1}-x_{2})H(x_{2})\Bigg{]},$

where $\bar{\Sigma}(x-y)=\Sigma(\nabla^{2})\delta(x-y)$ with $\Sigma$ being
the light-quark self-energy function in the coordinate space. We take the
argument of $\Sigma$ to be the covariant derivative in order to retain the
correct chiral transformation properties in the theory Yang  _et al._ (2002).
We have considered only the field $H$ in Eq. (LABEL:S2), while a similar
equation holds for $G$. Taking the derivative expansion up to the first order,
we obtain

$\displaystyle S_{2}$ $\displaystyle=$ $\displaystyle S^{(0)}_{2}+S^{(1)}_{2}$
$\displaystyle S^{(0)}_{2}$ $\displaystyle=$ $\displaystyle iN_{c}\int
d^{4}x\int\frac{d^{4}p}{(2\pi)^{4}}\left[-\frac{1}{p^{2}-\Sigma^{2}}+\frac{\Sigma}{(p^{2}-\Sigma^{2})v\cdot
p}\right]\mathrm{Tr}\left[\bar{H}(x)H(x)\right]$ $\displaystyle{}+iN_{c}\int
d^{4}x\int\frac{d^{4}p}{(2\pi)^{4}}\left[-\frac{1}{(p^{2}-\Sigma^{2})v\cdot
p}-\frac{2\Sigma}{(p^{2}-\Sigma^{2})^{2}}-\frac{2\Sigma^{\prime}(p^{2}+\Sigma^{2})}{(p^{2}-\Sigma^{2})^{2}}\right]\mathrm{Tr}\left[H(x)iv\cdot\partial\bar{H}(x)\right]$
$\displaystyle{}+iN_{c}\int
d^{4}x\int\frac{d^{4}p}{(2\pi)^{4}}\left[-\frac{2\Sigma^{\prime}(p^{2}+\Sigma^{2})}{(p^{2}-\Sigma^{2})^{2}}\right]\mathrm{Tr}\left[H(x)v\cdot
V_{\Omega}\bar{H}(x)\right],$ $\displaystyle S^{(1)}_{2}$ $\displaystyle=$
$\displaystyle iN_{c}\int
d^{4}x\int\frac{d^{4}p}{(2\pi)^{4}}\left[-\frac{1}{p^{2}-\Sigma^{2}}+\frac{\Sigma}{(p^{2}-\Sigma^{2})v\cdot
p}\right]\left[\frac{-p^{2}+(v\cdot
p)^{2}}{2m_{Q}}\right]\mathrm{Tr}\left[\bar{H}(x)H(x)\right]$ (16)
$\displaystyle{}+iN_{c}\int
d^{4}x\int\frac{d^{4}p}{(2\pi)^{4}}\biggr{[}{}-\frac{1}{(p^{2}-\Sigma^{2})v\cdot
p}-\frac{2\Sigma+2\Sigma^{\prime}(p^{2}+\Sigma^{2})-(4\Sigma\Sigma^{\prime}+2)v\cdot
p}{(p^{2}-\Sigma^{2})^{2}}\biggr{]}\biggr{[}\frac{-p^{2}+(v\cdot
p)^{2}}{2m_{Q}v\cdot p}\biggr{]}$
$\displaystyle\qquad\qquad\qquad\qquad\quad{}\times\mathrm{Tr}\left[H(x)iv\cdot\partial\bar{H}(x)\right]$
$\displaystyle{}+iN_{c}\int
d^{4}x\int\frac{d^{4}p}{(2\pi)^{4}}\left[-\frac{2\Sigma^{\prime}(p^{2}+\Sigma^{2})-4\Sigma\Sigma^{\prime}v\cdot
p}{(p^{2}-\Sigma^{2})^{2}}\right]\biggr{[}\frac{-p^{2}+(v\cdot
p)^{2}}{2m_{Q}v\cdot p}\biggr{]}\mathrm{Tr}\left[H(x)v\cdot
V_{\Omega}\bar{H}(x)\right],$

$S^{(0)}_{2}$ is the leading order in the $1/m_{Q}$ expansion (which survives
in the $m_{Q}\rightarrow\infty$ limit) and $S^{(1)}_{2}$ is the $1/m_{Q}$
order correction to $S^{(0)}_{2}$. The simplest terms of the heavy-light
mesons interacting with the Goldstone boson are generated by taking additional
derivative with respect to the chiral field $U(x)$. Since the $U$ field is
attached to the rotated external source, we actually take derivatives on the
action $S$ with respect to $J_{\Omega}$. The resultant expression is

$\displaystyle S_{3}$ $\displaystyle\equiv$
$\displaystyle\frac{\delta^{3}S}{\delta J_{\Omega}\delta\bar{H}\delta
H}J_{\Omega}H\bar{H},$ $\displaystyle=$ $\displaystyle{}-iN_{c}\int
d^{4}x_{1}d^{4}x_{2}d^{4}x_{3}$
$\displaystyle\qquad{}\times\mathrm{Tr}\Bigg{[}(i\not{\partial}-\Sigma)^{-1}\delta(x_{2}-x_{3})$
$\displaystyle\qquad\qquad{}\times
J_{\Omega}(x_{3})(i\not{\partial}-\Sigma)^{-1}\delta(x_{3}-x_{1})\bar{H}(x_{1})$
$\displaystyle\qquad\qquad{}\times\left(iv\cdot\partial-\frac{(\not{\partial}_{\perp})^{2}}{2m_{Q}}\right)^{-1}\delta(x_{1}-x_{2})H(x_{2})\Bigg{]}.$

The external source $J_{\Omega}$ contains the scalar part $S_{\Omega}$ and the
pseudo-scalar part $P_{\Omega}$, both of which are of $\mathcal{O}(p^{2})$ in
the chiral power counting. We only consider up to $\mathcal{O}(p)$ order
chiral Lagrangian, so that their contributions are ignored. Then we obtain

$\displaystyle S_{3}$ $\displaystyle=$ $\displaystyle
S^{(0)}_{3}+S^{(1)}_{3},$ $\displaystyle S^{(0)}_{3}$ $\displaystyle=$
$\displaystyle{}-iN_{c}\int
d^{4}x\int\frac{d^{4}p}{(2\pi)^{4}}\left[\frac{\Sigma^{2}+\frac{1}{3}p^{2}}{(p^{2}-\Sigma^{2})v\cdot
p}-\frac{2\Sigma}{(p^{2}-\Sigma^{2})^{2}}\right]\mathrm{Tr}\left[H(x)\gamma_{\mu}\gamma_{5}A^{\mu}_{\Omega}(x)\bar{H}(x)\right]$
$\displaystyle{}-iN_{c}\int
d^{4}x\int\frac{d^{4}p}{(2\pi)^{4}}\left[\frac{1}{(p^{2}-\Sigma^{2})v\cdot
p}+\frac{2\Sigma}{(p^{2}-\Sigma^{2})^{2}}\right]\mathrm{Tr}\left[H(x)v_{\mu}V^{\mu}_{\Omega}(x)\bar{H}(x)\right],$
$\displaystyle S^{(1)}_{3}$ $\displaystyle=$ $\displaystyle{}-iN_{c}\int
d^{4}x\int\frac{d^{4}p}{(2\pi)^{4}}\left[\frac{\Sigma^{2}+\frac{1}{3}p^{2}}{(p^{2}-\Sigma^{2})^{2}v\cdot
p}-\frac{2\Sigma-\frac{2}{3}v\cdot
p}{(p^{2}-\Sigma^{2})^{2}}\right]\biggr{[}\frac{-p^{2}+(v\cdot
p)^{2}}{2m_{Q}v\cdot
p}\biggr{]}\mathrm{Tr}\left[H(x)\gamma_{\mu}\gamma_{5}A^{\mu}_{\Omega}(x)\bar{H}(x)\right]$
(18) $\displaystyle{}-iN_{c}\int
d^{4}x\int\frac{d^{4}p}{(2\pi)^{4}}\left[\frac{1}{(p^{2}-\Sigma^{2})v\cdot
p}+\frac{2\Sigma-2v\cdot
p}{(p^{2}-\Sigma^{2})^{2}}\right]\biggr{[}\frac{-p^{2}+(v\cdot
p)^{2}}{2m_{Q}v\cdot
p}\biggr{]}\mathrm{Tr}\left[H(x)v_{\mu}V^{\mu}_{\Omega}(x)\bar{H}(x)\right].$

Again, $S^{(0)}_{3}$ represents the leading order term in the $1/m_{Q}$
expansion and $S^{(1)}_{3}$ is the $\mathcal{O}(1/m_{Q})$ term.

Summing up Eqs. (16) and (18), we obtain the expressions for the constants
$m_{H}$ and $g_{H}$ as

$\displaystyle m_{H}$ $\displaystyle=$
$\displaystyle\frac{iN_{c}}{Z_{H}}\int\frac{d^{4}p}{(2\pi)^{4}}\left[{}\frac{1}{p^{2}-\Sigma^{2}}-\frac{\Sigma}{(p^{2}-\Sigma^{2})v\cdot
p}\right]\biggr{[}1+\frac{-p^{2}+(v\cdot p)^{2}}{2m_{Q}v\cdot p}\biggr{]},$
$\displaystyle g_{H}$ $\displaystyle=$
$\displaystyle{}-\frac{iN_{c}}{Z_{H}}\int\frac{d^{4}p}{(2\pi)^{4}}\left[\frac{\Sigma^{2}+\frac{1}{3}p^{2}}{(p^{2}-\Sigma^{2})^{2}v\cdot
p}-\frac{2\Sigma}{(p^{2}-\Sigma^{2})^{2}}\right]\biggr{[}1+\frac{-p^{2}+(v\cdot
p)^{2}}{2m_{Q}v\cdot p}\biggr{]}+\frac{\frac{2}{3}[-p^{2}+(v\cdot
p)^{2}]}{2m_{Q}(p^{2}-\Sigma^{2})^{2}}.$ (19)

One can find that the coefficients of
$\mathrm{Tr}\left[H(x)iv\cdot\partial\bar{H}(x)\right]$ term and
$\mathrm{Tr}\left[H(x)v\cdot V_{\Omega}\bar{H}(x)\right]$ term are the same.
This means that the vector part of the chiral symmetry is reserved in our
approach, in agreement with the pattern of the chiral symmetry breaking in
QCD. According to our previous definition of the covariant derivative
$\nabla$, we can write them uniformly as
$\mathrm{Tr}\left[H(x)iv\cdot\nabla\bar{H}(x)\right]$. The coefficient of this
term is called the wave function renormalization factor which is given by

$\displaystyle Z_{H}$ $\displaystyle=$ $\displaystyle
iN_{c}\int\frac{d^{4}p}{(2\pi)^{4}}\biggr{[}{}-\frac{1}{(p^{2}-\Sigma^{2})v\cdot
p}-\frac{2\Sigma+2\Sigma^{\prime}(p^{2}+\Sigma^{2})}{(p^{2}-\Sigma^{2})^{2}}\biggr{]}\biggr{[}1+\frac{-p^{2}+(v\cdot
p)^{2}}{2m_{Q}v\cdot
p}\biggr{]}+\frac{(4\Sigma\Sigma^{\prime}+2)[-p^{2}+(v\cdot
p)^{2}]}{2m_{Q}(p^{2}-\Sigma^{2})^{2}}.$ (20)

Following the same procedure, we obtain the LECs for the field $G$ and the
coupling constant between the chiral partner fields $H$ and $G$, $g_{HG}$, as

$\displaystyle m_{G}$ $\displaystyle=$
$\displaystyle\frac{iN_{c}}{Z_{G}}\int\frac{d^{4}p}{(2\pi)^{4}}\left[\frac{1}{p^{2}-\Sigma^{2}}+\frac{\Sigma}{(p^{2}-\Sigma^{2})v\cdot
p}\right]\biggr{[}1+\frac{-p^{2}+(v\cdot p)^{2}}{2m_{Q}v\cdot p}\biggr{]},$
$\displaystyle g_{G}$ $\displaystyle=$
$\displaystyle{}-\frac{iN_{c}}{Z_{G}}\int\frac{d^{4}p}{(2\pi)^{4}}\left[\frac{\Sigma^{2}+\frac{1}{3}p^{2}}{(p^{2}-\Sigma^{2})^{2}v\cdot
p}+\frac{2\Sigma}{(p^{2}-\Sigma^{2})^{2}}\right]\biggr{[}1+\frac{-p^{2}+(v\cdot
p)^{2}}{2m_{Q}v\cdot p}\biggr{]}+\frac{\frac{2}{3}[-p^{2}+(v\cdot
p)^{2}]}{2m_{Q}(p^{2}-\Sigma^{2})^{2}}$ $\displaystyle Z_{G}$ $\displaystyle=$
$\displaystyle
iN_{c}\int\frac{d^{4}p}{(2\pi)^{4}}\biggr{[}{}-\frac{1}{(p^{2}-\Sigma^{2})v\cdot
p}+\frac{2\Sigma+2\Sigma^{\prime}(p^{2}+\Sigma^{2})}{(p^{2}-\Sigma^{2})^{2}}\biggr{]}\biggr{[}1+\frac{-p^{2}+(v\cdot
p)^{2}}{2m_{Q}v\cdot
p}\biggr{]}+\frac{(4\Sigma\Sigma^{\prime}+2)[-p^{2}+(v\cdot
p)^{2}]}{2m_{Q}(p^{2}-\Sigma^{2})^{2}}.$ $\displaystyle g_{HG}$
$\displaystyle=$
$\displaystyle{}-i\frac{N_{c}}{\sqrt{Z_{H}Z_{G}}}\int\frac{d^{4}p}{(2\pi)^{4}}\biggl{[}\frac{\Sigma^{2}+p^{2}}{(p^{2}-\Sigma^{2})^{2}v\cdot
p}\biggr{]}\biggr{[}1+\frac{-p^{2}+(v\cdot p)^{2}}{2m_{Q}v\cdot
p}\biggr{]}-\frac{2[-p^{2}+(v\cdot p)^{2}]}{2m_{Q}(p^{2}-\Sigma^{2})^{2}}.$
(21)

So far, we have obtained all the LECs of heavy-light mesons chiral effective
Lagrangian with $1/m_{Q}$ correction in Eqs. (19) - (21). The LECs are
expressed as integrals of the light-quark self-energy $\Sigma(-p^{2})$ and its
derivative, which can be calculated from QCD.

In deriving the effective action for the chiral effective Lagrangian, we have
taken a dynamical expansion which retains only the minimum of QCD dynamics
responsible for generating the dynamical chiral symmetry breaking. That’s why
the Gluon effects are taken into account only through the quark dressing
functions. This expansion has been found workable well for the calculation of
the pion decay constant Pagels and Stokar (1979). This simplification
unfortunately has a consequence that the mass splitting between the spin-0 and
spin-1 particles cannot be generated — the magnetic momentum operator of the
quark-gluon interaction cannot be generated. However, corrections could be
systematically calculated by taking into account higher order gluon Green
functions, which is beyond our present purpose.

##  IV Numerical results

Given the expressions of LECs shown in the previous section, one can calculate
the masses of the heavy-light mesons and the coupling constants, as long as
the light quark self-energy $\Sigma(-p^{2})$ is properly obtained. Following
our previous work Chen  _et al._ (2020), we use the self-energy obtained from
the Dyson-Schwinger equations as well as the lattice QCD to perform the
numerical calculations.

For the Dyson-Schwinger equation method, we use the differential form of the
gap equation Yang  _et al._ (2002):

$\displaystyle\left(\frac{\alpha_{s}(x)}{x}\right)^{\prime}\Sigma(x)^{\prime\prime}-\left(\frac{\alpha_{s}(x)}{x}\right)^{\prime\prime}\Sigma(x)^{\prime}$
$\displaystyle\qquad\qquad\qquad\quad{}-\frac{3C_{2}(R)}{4\pi}\frac{x\Sigma(x)}{x+\Sigma^{2}(x)}\left(\frac{\alpha_{s}(x)}{x}\right)^{\prime
2}=0,$

with boundary conditions

$\displaystyle\Sigma^{\prime}(0)$ $\displaystyle=$
$\displaystyle{}-\frac{3C_{2}(R)\alpha_{s}(0)}{8\pi\Sigma(0)},$
$\displaystyle\Sigma(\Lambda^{\prime})$ $\displaystyle=$
$\displaystyle\frac{3C_{2}(R)\alpha_{s}(\Lambda^{\prime})}{4\pi\Lambda^{\prime}}\int^{\Lambda^{\prime}}_{0}dx\frac{x\Sigma(x)}{x+\Sigma^{2}(x)},$
(22)

where $\alpha_{s}$ is the running coupling constant of QCD. $\Lambda^{\prime}$
is an ultraviolet cutoff regularizing the integral, which should be taken
$\Lambda^{\prime}\rightarrow\infty$ eventually. Since the low energy behavior
of the strong coupling constant $\alpha_{s}$ is not clear up to now, we take a
model description for $\alpha_{s}$ given in Ref. Dudal  _et al._ (2012), which
is called the refined Gribov-Zwanziger (G-Z) formalism: 222In our previous
work Chen  _et al._ (2020), another model was also adopted. However the model
gives similar results as those by using the self-energy function fitted from
Lattice QCD. So in this work, we shall not consider that model.

$\alpha_{s}(p^{2})=\alpha_{0}p^{2}\frac{p^{2}+M^{2}}{p^{4}+(M^{2}+m^{2})p^{2}+\lambda^{4}},$
(23)

where $M^{2}=4.303~{}{\rm GeV}^{2},(M^{2}+m^{2})=0.526~{}{\rm GeV}^{2}$ and
$\lambda^{4}=0.4929~{}{\rm GeV}^{4}$. $\alpha_{0}$ is a model parameter to be
determined. Although this model does not respect the UV behavior of QCD, since
the LECs are mostly controlled by the low energy behavior of QCD, this should
not be a problem in our study.

The LECs are calculated according to Eqs. (19)-(21). It is clear that the
integrals of the LECs have a physical ultraviolet cutoff $\Lambda_{c}$ which
should be of the order of the chiral symmetry breaking scale and serves as
another parameter in our calculations. When studying the low energy constants
in the heavy quark limit, we have found that there is a tension between the
pion decay constant $f_{\pi}$ and the mass splitting of the chiral partners
$\Delta m$, i.e., no parameter sets could fit both quantities perfectly Chen
_et al._ (2020). In this work, we find that the tension is largely released as
long as the $1/m_{Q}$ corrections are taken into account. Since $f_{\pi}$ is a
well-established quantity, we shall take $f_{\pi}\simeq 93$ MeV as an input to
determine $\Lambda_{c}$, and leave the parameter $a_{0}$ as a free one. The
parameter $a_{0}$ is scanned from $a_{0}=0.52$ through $a_{0}=0.62$, and we
find that $a_{0}=0.60$ gives the best fitted results. To give an intuitive
impression, we draw the running coupling constant $\alpha_{s}$ calculated with
the G-Z formalism at $a_{0}=0.60$ in Fig. 1. The light quark self-energy
$\Sigma(-p^{2})$ solved by the gap equation (IV) is shown in Fig. 2. In Table
1, we list the LECs with and without $1/m_{Q}$ corrections.

Figure 1: Running coupling constant of the G-Z model with $a_{0}=0.60$.  Figure 2: The lattice fittings of $M(-p^{2})$ and $Z(-p^{2})$ given in Ref. Oliveira  _et al._ (2019) and $\Sigma(-p^{2})$ from the gap equation with Model G-Z.  Table 1: The heavy-light meson masses and the coupling constants calculated from the G-Z Model with $a_{0}=0.6$.  Correction |  $\Sigma(0)$ (GeV) |  $f_{\pi}$ (GeV) |  ${}-\langle\bar{\psi}\psi\rangle$ (GeV3) |  $g_{G}$ |  $g_{H}$ |  $g_{HG}$ |  $m_{G}$ (GeV) |  $m_{H}$ (GeV) |  $\Delta m$ (GeV)
---|---|---|---|---|---|---|---|---|---
LO | 0.484 | 0.093 | 0.261 | 1.570 | -0.576 | 0.022 | 0.963 | 0.662 | 0.301
$1/m_{c}$ | 0.484 | 0.093 | 0.261 | 1.862 | -0.554 | -0.139 | 1.142 | 0.739 | 0.403
$1/m_{b}$ | 0.484 | 0.093 | 0.261 | 1.649 | -0.570 | -0.022 | 1.012 | 0.683 | 0.329

From Table 1, we find that the $1/m_{Q}$ corrections for the LECs except
$g_{HG}$ are ranging from $\sim 1\%-5\%$ for bottom mesons and from $\sim
4\%-19\%$ for charmed mesons. This agrees with the fact that, compare to the
charm quark sector, the heavy quark expansion is more valid in the bottom
quark sector. It is interesting to notice that the corrections to the LECs
associated to $H$ fields are typically smaller than those associated to $G$
fields by a few percents. It is found that the mass splitting of the chiral
partners $\Delta m$ which is explicitly related to the dynamical chiral
symmetry breaking of QCD in our method suffers a significant correction for
the charmed mesons, and the resulted value is consistent with the experimental
data (recalling that the averaged mass splitting for the $D$ meson sector is
$\Delta m=0.430$ GeV.) The mass splitting for the $B$ meson sector suffers the
$1/m_{b}$ correction by about $9\%$. According to our results, we may conclude
that, for the charmed mesons, the $1/m_{Q}$ corrections usually contribute
significantly, while for the bottom mesons, taking the heavy quark limit may
cause errors up to $\sim 10$ percents.

The masses $m_{H}$ and $m_{G}$ displayed in Table 1 are the “residue masses”
with the heavy quark mass rotated away. The physical masses for the $H$ and
$G$ doublets can be easily obtained by restoring the heavy quark mass $m_{c}$
or $m_{b}$. Using $m_{c}\approx 1.27$ GeV and $m_{b}\approx 4.66$ GeV Zyla
_et al._ (2020), we obtain the physical masses for the $D$ meson family to be
$\tilde{m}_{H}\approx 2.01$ GeV and $\tilde{m}_{G}\approx 2.41$ GeV; and for
the $B$ meson family to be $\tilde{m}_{H}\approx 5.34$ GeV and
$\tilde{m}_{G}\approx 5.67$ GeV. These results are consistent with the
experimental data for the charmed mesons (4) and the bottom mesons (6)
respectively.

The coupling constant $g_{H}$ directly governs the decay process $D^{\ast}\to
D\pi$. Given the value shown in Table 1, we find the decay width

$\displaystyle\Gamma(D^{\ast+}\to
D\pi)=\frac{3}{2}\frac{g^{2}_{H}|P_{\pi}|^{3}}{12\pi
f_{\pi}^{2}}=87.2~{}\text{KeV},$ (24)

which is very close to the experimental result $83.4\pm 1.8$ KeV Zyla  _et
al._ (2020).

We now turn to the calculation using the light quark self-energy obtained from
lattice QCD. In Ref. Oliveira  _et al._ (2019), the authors fitted the lattice
results for the quark wave function renormalization $Z(-p^{2})$ and the
running quark mass $M(-p^{2})$ (corresponding to $\Sigma(-p^{2})$ in our
terminology). For illustration, we plot $M(-p^{2})$ and $Z(-p^{2})$ in Fig. 2.
In the previous section, we present the formula for LECs with the wave
function renormalization $Z(-p^{2})$ ignored. It is straightforward to keep
$Z(-p^{2})$ in the formula which are given in Appendix B. For comparison, we
performed the calculations using the fitted functions from lattice QCD with
$Z(-p^{2})$ (see Table 2) and without $Z(-p^{2})$ (see Table 3). One can see
from Table 3 that the mass splitting is significantly deviated from the
expected values in the $Z(-p^{2})=1$ case. The reason may be understood that
in the lattice QCD, both $M(-p^{2})$ and $Z(-p^{2})$ are needed to describe
the dressing effects of the quark propagator, while in the G-Z model, the
running quark mass $M(-p^{2})$ (i.e., $\Sigma(-p^{2})$ in the G-Z model) alone
satisfies the Dyson-Schwinger equation for the quark propagator.

Table 2: LECs calculated from lattice fittings given in Ref Oliveira  _et al._ (2019).  Correction |  $M(0)$(GeV) |  $f_{\pi}$ (GeV) |  ${}-\langle\bar{\psi}\psi\rangle$ (GeV3) |  $g_{G}$ |  $g_{H}$ |  $g_{HG}$ |  $m_{G}$ (GeV) |  $m_{H}$ (GeV) |  $\Delta m$ (GeV)
---|---|---|---|---|---|---|---|---|---
LO | 0.418 | 0.093 | 0.307 | 0.820 | -0.095 | 0.805 | 1.012 | 0.660 | 0.352
$1/m_{c}$ | 0.418 | 0.093 | 0.307 | 1.116 | -0.032 | 0.696 | 1.246 | 0.806 | 0.440
$1/m_{b}$ | 0.418 | 0.093 | 0.307 | 0.901 | -0.078 | 0.776 | 1.076 | 0.700 | 0.376
Table 3: LECs calculated from lattice fittings given in Ref Oliveira  _et al._ (2019) with $Z(-p^{2})=1$.  Correction |  $M(0)$(GeV) |  $f_{\pi}$ (GeV) |  ${}-\langle\bar{\psi}\psi\rangle$ (GeV3) |  $g_{G}$ |  $g_{H}$ |  $g_{HG}$ |  $m_{G}$ (GeV) |  $m_{H}$ (GeV) |  $\Delta m$ (GeV)
---|---|---|---|---|---|---|---|---|---
LO | 0.418 | 0.093 | 0.307 | 1.016 | -0.137 | 0.819 | 1.183 | 0.611 | 0.572
$1/m_{c}$ | 0.418 | 0.093 | 0.307 | 1.556 | -0.087 | 0.749 | 1.657 | 0.710 | 0.947
$1/m_{b}$ | 0.418 | 0.093 | 0.307 | 1.163 | -0.123 | 0.800 | 1.312 | 0.638 | 0.674

From Table 2, we see that the $1/m_{Q}$ corrections to the LECs are even
larger compared to the G-Z model results. The mass splitting for $D$ and $B$
mesons is consistent with the experimental results respectively. The physical
masses are $\tilde{m}_{H}\approx 2.08$ GeV and $\tilde{m}_{G}\approx 2.52$ GeV
for $D$ mesons, and $\tilde{m}_{H}\approx 5.36$ GeV and $\tilde{m}_{G}\approx
5.74$ GeV for $B$ mesons. The masses of $D$ mesons are slightly larger than
the experimental data (Eq. (4)), while the masses of $B$ mesons are consistent
with the experimental data (Eq. (6)). We also notice from Table 2 that the
coupling constant $g_{H}$ is too small to fit the experimental data, which
indicates the shortcomings of the results based on lattice fittings. In
general, coupling constants in the effective Lagrangian are more sensitive to
the details of $M(-p^{2})$ and $Z(-p^{2})$ than the masses are. Since the
functions $Z(-p^{2})$ and $M(-p^{2})$ are fitted from the lattice data
corresponding to $M_{\pi}=295$ MeV Oliveira  _et al._ (2019), we are not
expecting them could reproduce every LEC properly.

##  V Summary

In this paper, we extend our previous work on deriving the chiral effective
Lagrangian for heavy-light mesons from QCD to include $1/m_{Q}$ corrections.
Using the quark self-energy and the wavefunction renormalization calculated
from Dyson-Schwinger equations as well as those fitted from lattice QCD, we
calculated the $1/m_{Q}$ corrections to the LECs in the Lagrangian. It is
found that for charmed mesons, the $1/m_{Q}$ corrections are significant and
typically give about $\sim 20\%$ contributions, while for bottom mesons, the
$1/m_{Q}$ corrections are usually ranging from a few percents up to $\sim
10\%$.

The $1/m_{Q}$ order corrections improve the leading order results. At the
leading order in the $1/m_{Q}$ expansion, the pion decay constant $f_{\pi}$
and the mass splitting of the chiral partners $\Delta m$ cannot be fitted to
the corresponding experimental data simultaneously. However, with the
$1/m_{Q}$ contributions included, our numerical results from both the G-Z
model and the lattice fittings are comparable to the corresponding
experimental data. In addition, in the G-Z model, most of the LECs at
$1/m_{Q}$ order are consistent with the existent experimental results.
Moreover, it is also find that the $D^{\ast}\to D\pi$ decay width is improved
to $\Gamma(D^{\ast+}\to D\pi)=87.2~{}\text{KeV}$ which is consistent to the
experimental result.

###### Acknowledgements.

The work of Y. L. M. was supported in part by National Science Foundation of
China (NSFC) under Grant No. 11875147 and No.11475071. Q. W. was supported by
the National Science Foundation of China (NSFC) under Grant No. 11475092.

##  Appendix A DERIVATION OF THE ACTION

We derive the effective action at the $1/m_{Q}$ order from QCD in this
appendix.

###  A.1 Integrate in the Nambu-Goldstone boson fields

Starting from Eq. (12) in Sec. III, we now introduce the pseudoscalar meson
field $U$ into the theory. Inserting the following constant into the
generating functional

$\displaystyle I$ $\displaystyle=$
$\displaystyle\int\mathcal{D}U\delta(U^{\dagger}U-1)\delta(\det U-1)$
$\displaystyle\quad{}\times\mathcal{F}[\mathcal{O}]\delta\left(\Omega\mathcal{O}^{\dagger}\Omega-\Omega^{\dagger}\mathcal{O}\Omega^{\dagger}\right),$

where
$\mathcal{F}[\mathcal{O}]\equiv\\{\prod\limits_{x}\det\mathcal{O}\int\mathcal{D}\sigma\delta(\mathcal{O}^{\dagger}\mathcal{O}-\sigma^{\dagger}\sigma)\delta(\sigma-\sigma^{\dagger})\\}^{-1}$,
$\mathcal{O}(x)=e^{-i[\theta(x)/N_{f}]}\mathrm{tr}_{l}[P_{R}B^{T}(x,x)]$, and
$B$ is the abbreviation of the bilocal composite light quark fields. Then,
integrating out the gluon fields, we obtain

$\displaystyle Z[J]$ $\displaystyle=$
$\displaystyle\int\mathcal{D}q\mathcal{D}\bar{q}\mathcal{D}N_{v}\mathcal{D}\bar{N}_{v}\mathcal{D}\bar{U}\delta(\mathcal{O}^{\dagger}-\mathcal{O})$
$\displaystyle{}\times\exp\biggl{\\{}i\int
d^{4}x\biggl{[}\bar{q}(i\not{\partial}+J_{\Omega})q+\bar{N}_{v}\left(iv\cdot\partial-\frac{\partial^{2}-(v\cdot\partial)^{2}}{2m_{Q}}\right)N_{v}\biggr{]}+i\Gamma_{I}[B]$
$\displaystyle\qquad\quad{}+\sum\limits_{n=2}^{\infty}\int
d^{4}x_{1}d^{4}x^{\prime}_{1}\cdots d^{4}x_{n}d^{4}x^{\prime}_{n}$
$\displaystyle\qquad\qquad\qquad{}\times\dfrac{(-i)^{n}(N_{c}g^{2})^{n-1}}{n!}\bar{G}^{\sigma_{1}\cdots\sigma_{n}}_{\rho_{1}\cdots\rho_{n}}(x_{1},x^{\prime}_{1},\cdots,x_{n},x^{\prime}_{n})\bar{\psi}^{\sigma_{1}}_{\alpha_{1}}(x_{1})\psi^{\rho_{1}}_{\alpha_{1}}(x^{\prime}_{1})\cdots\bar{\psi}^{\sigma_{n}}_{\alpha_{n}}(x_{n})\psi^{\rho_{n}}_{\alpha_{n}}(x^{\prime}_{n})\biggl{\\}},$

where the Fierz reordering has been made and the extended gluon Green
functions
$\bar{G}^{\sigma_{1}\cdots\sigma_{n}}_{\rho_{1}\cdots\rho_{n}}(x_{1},x^{\prime}_{1},\cdots,x_{n},x^{\prime}_{n})$
have been introduced,
$\mathcal{D}\bar{U}=\mathcal{D}U\delta(U^{\dagger}U-1)\delta(\det U-1)$ and
$e^{i\Gamma_{I}[B]}\equiv\mathcal{F}[\mathcal{O}]$. The details of Fierz
reordering and introducing the extended gluon Green functions could be found
in Refs Wang  _et al._ (2000a); Fu  _et al._ (2017). Due to the $1/m_{Q}$
corrections, $\bar{G}$’s include more rich types of gluon Green functions than
regular ones. The field $\psi$ stands for both the light and heavy quarks,
i.e., $\psi=(q,N_{v})$. We have made chiral rotation of the original light
quark fields and ignored the chiral anomaly. $J_{\Omega}(x)$ is the external
source under chiral rotation Wang  _et al._ (2000a)

$\displaystyle J_{\Omega}(x)$ $\displaystyle=$
$\displaystyle[\Omega(x)P_{R}+\Omega^{\dagger}(x)P_{L}]$ (26)
$\displaystyle\times[J(x)+i\not{\partial}][\Omega(x)P_{R}+\Omega^{\dagger}(x)P_{L}].$

###  A.2 Integrate in the heavy-light meson fields

The heavy-light meson fields may be introduced by inserting the following
constant

$\int\mathcal{D}\Phi\delta(N_{c}\Phi^{ab}_{\eta\xi}(x,x^{\prime})-\bar{\psi}^{a}_{\eta}(x)\psi^{b}_{\xi}(x^{\prime})),$
(27)

into Eq. (A.1), where $\Phi^{ab}_{\eta\xi}$ is a bilocal auxiliary field. The
$\delta$ function can be further expressed in the Fourier representation

$\displaystyle\delta{(N_{c}\Phi^{ab}_{\eta\xi}(x,x^{\prime})-\bar{\psi}^{a}_{\eta}(x)\psi^{b}_{\xi}(x^{\prime}))}\sim\int\mathcal{D}\Pi\exp^{i\int
d^{4}xd^{4}x^{\prime}\Pi^{ab}_{\eta\xi}(x,x^{\prime})[N_{c}\Phi^{ab}_{\eta\xi}(x,x^{\prime})-\bar{\psi}^{a}_{\eta}(x)\psi^{b}_{\xi}(x^{\prime})]}.$
(28)

Then the generating functional becomes

$\displaystyle Z[J]$ $\displaystyle=$
$\displaystyle\int\mathcal{D}q\mathcal{D}\bar{q}\mathcal{D}N_{v}\mathcal{D}\bar{N}_{v}\mathcal{D}\bar{U}\mathcal{D}\Phi\mathcal{D}\Pi\delta(\mathcal{O}^{\dagger}-\mathcal{O})$
(29) $\displaystyle{}\times\exp\biggl{\\{}i\int
d^{4}x\biggl{[}\bar{q}(i\not{\partial}+J_{\Omega})q+\bar{N}_{v}\left(iv\cdot\partial-\frac{\partial^{2}-(v\cdot\partial)^{2}}{2m_{Q}}\right)N_{v}\biggr{]}+i\Gamma_{I}[\Phi]+iN_{c}\bar{G}(\Phi)$
$\displaystyle{}\qquad\qquad+i\int
d^{4}xd^{4}x^{\prime}\biggl{[}N_{c}\Pi^{ab}_{\eta\xi}(x,x^{\prime})\Phi^{ab}_{\eta\xi}(x,x^{\prime})-\bar{\psi}(x)\Pi(x,x^{\prime})\psi(x^{\prime})\biggr{]}\biggr{\\}},$

where

$\displaystyle\bar{G}(\Phi)=\sum\limits_{n=2}^{\infty}\int
d^{4}x_{1}d^{4}x^{\prime}_{1}\cdots
d^{4}x_{n}d^{4}x^{\prime}_{n}\dfrac{(-i)^{n}(N_{c}g^{2})^{n-1}}{n!}\bar{G}^{\sigma_{1}\cdots\sigma_{n}}_{\rho_{1}\cdots\rho_{n}}(x_{1},x^{\prime}_{1},\cdots,x_{n},x^{\prime}_{n})\Phi^{\sigma_{1}\rho_{1}}(x_{1},x^{\prime}_{1})\cdots\Phi^{\sigma_{n}\rho_{n}}(x_{n},x_{n}^{\prime}).$

By integrating out the quark fields $\psi$ and $\bar{\psi}$, we obtain

$\displaystyle Z[J]$ $\displaystyle=$
$\displaystyle\int\mathcal{D}\bar{U}\mathcal{D}\Phi\mathcal{D}\Pi\delta(\mathcal{O}^{\dagger}-\mathcal{O})$
$\displaystyle{}\times\exp\biggl{\\{}iN_{c}\biggl{[}-i\mathrm{Tr}^{\prime}\ln[i\not{\partial}I_{1}+J^{\prime}_{\Omega}+\left(iv\cdot\partial-\frac{\partial^{2}-(v\cdot\partial)^{2}}{2m_{Q}}\right)I_{4}-\Pi]\biggr{]}+iN_{c}\mathrm{Tr}^{\prime}(\Pi\Phi^{T})+iN_{c}\bar{G}(\Phi)+i\Gamma_{I}[\Phi]\biggr{\\}},$

where we have defined the functional trace $\mathrm{Tr}^{\prime}$ taking over
the flavor space, spinor space and coordinate space,
$I_{1}=\mathrm{diag}(1,1,0)$ and $I_{4}=\mathrm{diag}(0,0,1)$ are the matrices
in the flavor space. $J^{\prime}_{\Omega}$ is the extension of $J_{\Omega}$ to
the whole flavor space.

In flavor space the matrix $\Pi$ can be decomposed as

$\displaystyle\Pi=\Pi_{1}+\Pi_{2}+\Pi_{3}+\Pi_{4},$ (31)

where

$\displaystyle\Pi_{1}=\begin{pmatrix}\Pi^{qq}_{2\times 2}&0\\\
0&0\end{pmatrix},\quad\Pi_{2}=\begin{pmatrix}0&0\\\ \Pi^{Qq}_{1\times
2}&0\end{pmatrix},$
$\displaystyle\quad\Pi_{3}=\begin{pmatrix}0&\Pi^{qQ}_{2\times 1}\\\
0&0\end{pmatrix},\quad\Pi_{4}=\begin{pmatrix}0&0\\\ 0&\Pi^{QQ}\end{pmatrix}.$
(32)

Given this decomposition, we recognize that $\Pi_{2}$ and $\Pi_{3}$ are the
bosonic interpolating fields for the heavy-light mesons. To get a local
effective Lagrangian, we take the following localization conditions which are
essentially the point coupling between quarks Tandy (1997)

$\displaystyle\Pi_{2}(x,y)$ $\displaystyle=$
$\displaystyle\Pi_{2}(x)\delta(x-y),$ $\displaystyle\Pi_{3}(x,y)$
$\displaystyle=$ $\displaystyle\Pi_{3}(x)\delta(x-y).$ (33)

It is easy to check that
$\bar{\Pi}_{2}(x,x)\equiv\gamma^{0}\Pi_{2}^{\dagger}(x,x)\gamma^{0}=\Pi_{3}(x,x)$,
then the generating functional can be written as

$\displaystyle Z[J]$ $\displaystyle=$
$\displaystyle\int\mathcal{D}\bar{U}\mathcal{D}\Phi\mathcal{D}\Pi\exp\Biggl{\\{}i\Gamma_{I}[\Phi]+iN_{c}\biggr{\\{}\mathrm{Tr}^{\prime}(\Pi\Phi^{T})+\bar{G}(\Phi)+\int
d^{4}x\mathrm{Tr}\left\\{\Xi\left[\left(-i\sin\frac{\theta}{N_{f}}+\cos\frac{\theta}{N_{f}}\right)\Phi^{T}\right]\right\\}$
(34)
$\displaystyle{}\qquad\qquad\qquad\qquad\qquad-i\mathrm{Tr}^{\prime}\ln\left[i\not{\partial}I_{1}+J^{\prime}_{\Omega}+\left(iv\cdot\partial-\frac{\partial^{2}-(v\cdot\partial)^{2}}{2m_{Q}}\right)I_{4}-\Pi_{1}-\Pi_{2}-\bar{\Pi}_{2}-\Pi_{4}\right]\biggr{\\}}\Biggr{\\}}$

where the $\delta(\mathcal{O}^{\dagger}-\mathcal{O})$ term has been
reexpressed as

$\displaystyle\delta(\mathcal{O}^{\dagger}-\mathcal{O})=\int\mathcal{D}\Xi
e^{iN_{c}\int d^{4}x\mathrm{Tr}\\{\Xi(x)[\Theta\Phi^{T}(x,x)]\\}}$ (35)

with $\Xi^{ab}(x)$ being a new auxiliary field and
$\Theta\equiv(-i\sin\theta(x)/N_{f}+\cos\theta(x)/N_{f}).$

Now, by integrating out the fields $\Phi$, $\Xi$, $\Pi_{1}$ and $\Pi_{4}$, we
can obtain the action, denoted as $S[U,\Pi_{2},\bar{\Pi}_{2}]$, for the chiral
effective theory with heavy-light mesons

$\displaystyle
Z[J]=\int\mathcal{D}\bar{U}\mathcal{D}\Pi_{2}\mathcal{D}\bar{\Pi}_{2}\exp\\{iS[U,\Pi_{2},\bar{\Pi}_{2}]\\},$
(36)

where

$\displaystyle e^{iS}$ $\displaystyle\equiv$
$\displaystyle\int\mathcal{D}\Phi\mathcal{D}\Pi_{1}\mathcal{D}\Pi_{4}\mathcal{D}\Xi\exp\Biggl{\\{}i\Gamma_{I}[\Phi]+iN_{c}\biggr{\\{}\mathrm{Tr}^{\prime}(\Pi\Phi^{T})+\bar{G}(\Phi)+\int
d^{4}x\mathrm{Tr}\\{\Xi[(-i\sin\frac{\theta}{N_{f}}+\cos\frac{\theta}{N_{f}})\Phi^{T}]\\}$
$\displaystyle{}\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\quad-i\mathrm{Tr}^{\prime}\ln\left[i\not{\partial}I_{1}+J^{\prime}_{\Omega}+\left(iv\cdot\partial-\frac{\partial^{2}-(v\cdot\partial)^{2}}{2m_{Q}}\right)I_{4}-\Pi_{1}-\Pi_{2}-\bar{\Pi}_{2}-\Pi_{4}\right]\biggr{\\}}\Biggr{\\}}.$

###  A.3 The action in the large $N_{c}$ Limit

Now, we take the large $N_{c}$ limit in the generating functional to obtain

$\displaystyle S[U,\Pi_{2},\bar{\Pi}_{2}]$ $\displaystyle=$ $\displaystyle
N_{c}\Bigg{\\{}\mathrm{Tr}^{\prime}(\Pi_{c}\Phi^{T}_{c})+\bar{G}(\Phi_{c})+\int
d^{4}x\mathrm{Tr}\left\\{\Xi_{c}\left[\left(-i\sin\frac{\theta(\Phi_{c})}{N_{f}}+\cos\frac{\theta(\Phi_{c})}{N_{f}}\right)\Phi^{T}_{c}\right]\right\\}$
(38)
$\displaystyle{}\qquad-i\mathrm{Tr}^{\prime}\ln\left[i\not{\partial}I_{1}+J^{\prime}_{\Omega}+\left(iv\cdot\partial-\frac{\partial^{2}-(v\cdot\partial)^{2}}{2m_{Q}}\right)I_{4}-\Pi_{1c}-\Pi_{2}-\bar{\Pi}_{2}-\Pi_{4c}\right]\Bigg{\\}},$

where $\Phi_{c},\Xi_{c},\Pi_{1c},\Pi_{4c}$ are classical fields satisfying the
saddle point equations

$\frac{\delta S}{\delta\Phi_{c}}=\frac{\delta S}{\delta\Xi_{c}}=\frac{\delta
S}{\delta\Pi_{1c}}=\frac{\delta S}{\delta\Pi_{4c}}=0.$ (39)

and $\Gamma_{I}$ term has been ignored because it is of $O(1/N_{c})$ Wang  _et
al._ (2000a).

These saddle point equations provide important information. For instance,
equations $\frac{\delta S}{\delta\Pi_{1c}}=0$ and $\frac{\delta
S}{\delta\Phi_{1c}}=0$ generate the coupled equations

$\Phi^{Tab}_{1c\eta\xi}(x,y)=-i\left[(i\not{\partial}+J^{\prime}_{\Omega}-\Pi_{1c})^{-1}\right]^{ab}_{\eta\xi}(x,y),$
(40)

with

$\displaystyle\Pi^{ab}_{1c\eta\xi}(x,y)$ $\displaystyle=$
$\displaystyle{}-\sum\limits_{n=1}^{\infty}\int d^{4}x_{1}\cdots
d^{4}x^{\prime}_{n}d^{4}x^{\prime}_{1}\cdots
d^{4}x^{\prime}_{n}\dfrac{(-i)^{n+1}(N_{c}g^{2})^{n}}{n!}$ (41)
$\displaystyle{}\qquad\qquad\times\bar{G}_{\eta\eta_{1}\cdots\eta_{n},\xi\xi_{1}\cdots\xi_{n}}^{aa_{1}\cdots
a_{n},bb_{1}\cdots
b_{n}}(x,y,x_{1},x^{\prime}_{1},\cdots,x_{n},x^{\prime}_{n})\Phi^{a_{1}b_{1}}_{1c\eta_{1}\xi_{1}}(x_{1},x^{\prime}_{1})\cdots\Phi^{a_{n}b_{n}}_{1c\eta_{n}\xi_{n}}(x_{n},x^{\prime}_{n}),$

where the term involving $\Xi$ has been omitted because it vanishes once the
external sources $J$ are eventually turned off Wang  _et al._ (2000a). When
$J_{\Omega}^{\prime}$ is turned off, the coupled equations (40) and (41) are
nothing but the DSEs for the quark propagators with $\Pi_{1c}$ being the self-
energy for light quarks. So that we rewrite $\Pi_{1c}(x,y)$ as
$\bar{\Sigma}(x-y)I_{1}$. Similarly, $\Pi_{4c}$ is just the self-energy for
the heavy quark. Since the contribution from the heavy quark self-energy is
less significant than the light ones, we simply ignore it.

Because it is difficult to solve classical fields from the saddle point
equations (39) without approximations, we can not get $U,\Pi_{2}$ and
$\bar{\Pi}_{2}$ dependence of the classical fields and further obtain the
chiral effective Lagrangian. Thus, we follow Ref. Yang  _et al._ (2002) and
keep only the “$\mathrm{Tr}\ln$” term in the action in the spirit of the
dynamical perturbation which works well in the calculation of pion decay
constant Pagels and Stokar (1979). Eventually, the effective action is
simplified to be

$\displaystyle S[U,\Pi_{2},\bar{\Pi}_{2}]$ $\displaystyle=$
$\displaystyle{}-iN_{c}\mathrm{Tr}^{\prime}\ln\left[(i\not{\partial}-\bar{\Sigma})I_{1}+J^{\prime}_{\Omega}+\left(iv\cdot\partial-\frac{\partial^{2}-(v\cdot\partial)^{2}}{2m_{Q}}\right)I_{4}-\Pi_{2}-\bar{\Pi}_{2}\right].$
(42)

##  Appendix B FORMULA FOR LECS INCLUDING $Z(-p^{2})$

In order to include the quark wave function renormalization $Z(-p^{2})$ into
our system, we consider the general quark propagator as

$\displaystyle S(p)$ $\displaystyle=$ $\displaystyle\frac{i}{A(-p^{2})\not{p}\
-B(-p^{2})}$ $\displaystyle=$ $\displaystyle i\,Z(-p^{2})\,\frac{\not{p}\
+M(-p^{2})}{p^{2}-M^{2}(-p^{2})},$

where $Z(-p^{2})=1/A(-p^{2})$ stands for the quark wave function
renormalization and $M(-p^{2})=B(-p^{2})/A(-p^{2})$ is the renormalization
group invariant running quark mass. After a series of calculations, we can get
the LECs with $Z(-p^{2})$ as follows:

$\displaystyle m_{H}$ $\displaystyle=$
$\displaystyle\frac{iN_{c}}{Z_{H}}\int\frac{d^{4}p}{(2\pi)^{4}}\biggr{\\{}\biggr{[}{}\frac{1}{p^{2}-M^{2}}-\frac{M}{(p^{2}-M^{2})v\cdot
p}\biggr{]}\biggr{[}1+\frac{-p^{2}+(v\cdot p)^{2}}{2m_{Q}v\cdot
p}\biggr{]}\biggr{\\}}Z,$ $\displaystyle g_{H}$ $\displaystyle=$
$\displaystyle-\frac{iN_{c}}{Z_{H}}\int\frac{d^{4}p}{(2\pi)^{4}}\biggr{\\{}\biggl{[}\frac{M^{2}+\frac{1}{3}p^{2}}{(p^{2}-M^{2})^{2}v\cdot
p}-\frac{2M}{(p^{2}-M^{2})^{2}}\biggr{]}\biggr{[}1+\frac{-p^{2}+(v\cdot
p)^{2}}{2m_{Q}v\cdot p}\biggr{]}+\frac{\frac{2}{3}[-p^{2}+(v\cdot
p)^{2}]}{2m_{Q}(p^{2}-M^{2})^{2}}\biggr{\\}}Z^{2},$ $\displaystyle m_{G}$
$\displaystyle=$
$\displaystyle\frac{iN_{c}}{Z_{G}}\int\frac{d^{4}p}{(2\pi)^{4}}\biggr{\\{}\biggr{[}\frac{1}{p^{2}-M^{2}}+\frac{M}{(p^{2}-M^{2})v\cdot
p}\biggr{]}\biggr{[}1+\frac{-p^{2}+(v\cdot p)^{2}}{2m_{Q}v\cdot
p}\biggr{]}\biggr{\\}}Z,$ $\displaystyle g_{G}$ $\displaystyle=$
$\displaystyle-\frac{iN_{c}}{Z_{G}}\int\frac{d^{4}p}{(2\pi)^{4}}\biggr{\\{}\biggl{[}\frac{M^{2}+\frac{1}{3}p^{2}}{(p^{2}-M^{2})^{2}v\cdot
p}+\frac{2M}{(p^{2}-M^{2})^{2}}\biggr{]}\biggr{[}1+\frac{-p^{2}+(v\cdot
p)^{2}}{2m_{Q}v\cdot p}\biggr{]}+\frac{\frac{2}{3}[-p^{2}+(v\cdot
p)^{2}]}{2m_{Q}(p^{2}-M^{2})^{2}}\biggr{\\}}Z^{2},$ $\displaystyle Z_{H}$
$\displaystyle=$ $\displaystyle
iN_{c}\int\frac{d^{4}p}{(2\pi)^{4}}\biggr{\\{}\biggr{[}\left(-\frac{1}{(p^{2}-M^{2})v\cdot
p}-\frac{2M+2M^{\prime}(p^{2}+M^{2})}{(p^{2}-M^{2})^{2}}\right)Z-\frac{2Z^{\prime}M}{p^{2}-M^{2}}\biggr{]}\biggr{[}1+\frac{-p^{2}+(v\cdot
p)^{2}}{2m_{Q}v\cdot p}\biggr{]}$
$\displaystyle\qquad\qquad\qquad+\frac{Z(4MM^{\prime}+2)[-p^{2}+(v\cdot
p)^{2}]}{2m_{Q}(p^{2}-M^{2})^{2}}+\frac{2Z^{\prime}[-p^{2}+(v\cdot
p)^{2}]}{2m_{Q}(p^{2}-M^{2})}\biggr{\\}}$ $\displaystyle Z_{G}$
$\displaystyle=$ $\displaystyle
iN_{c}\int\frac{d^{4}p}{(2\pi)^{4}}\biggr{\\{}\biggr{[}\left(-\frac{1}{(p^{2}-M^{2})v\cdot
p}+\frac{2M+2M^{\prime}(p^{2}+M^{2})}{(p^{2}-M^{2})^{2}}\right)Z+\frac{2Z^{\prime}M}{p^{2}-M^{2}}\biggr{]}\biggr{[}1+\frac{-p^{2}+(v\cdot
p)^{2}}{2m_{Q}v\cdot p}\biggr{]}$
$\displaystyle\qquad\qquad\qquad+\frac{Z(4MM^{\prime}+2)[-p^{2}+(v\cdot
p)^{2}]}{2m_{Q}(p^{2}-M^{2})^{2}}-\frac{2Z^{\prime}[-p^{2}+(v\cdot
p)^{2}]}{2m_{Q}(p^{2}-M^{2})}\biggr{\\}}$ $\displaystyle g_{HG}$
$\displaystyle=$
$\displaystyle-i\frac{N_{c}}{\sqrt{Z_{H}Z_{G}}}\int\frac{d^{4}p}{(2\pi)^{4}}\biggr{\\{}\biggl{[}\frac{M^{2}+p^{2}}{(p^{2}-M^{2})^{2}v\cdot
p}\biggr{]}\biggr{[}1+\frac{-p^{2}+(v\cdot p)^{2}}{2m_{Q}v\cdot
p}\biggr{]}-\frac{2[-p^{2}+(v\cdot
p)^{2}]}{2m_{Q}(p^{2}-M^{2})^{2}}\biggr{\\}}Z^{2}.$ (43)

## References

  * Wang  _et al._ (2000a) Q. Wang, Y.-P. Kuang, M. Xiao, and X.-L. Wang, Phys. Rev. D 61, 054011 (2000a).
  * Yang  _et al._ (2002) H. Yang, Q. Wang, Y.-P. Kuang, and Q. Lu, Phys. Rev. D 66, 014019 (2002).
  * Jiang  _et al._ (2010) S.-Z. Jiang, Y. Zhang, C. Li, and Q. Wang, Phys. Rev. D 81, 014001 (2010).
  * Jiang and Wang (2010) S.-Z. Jiang and Q. Wang, Phys. Rev. D 81, 094037 (2010).
  * Jiang  _et al._ (2013) S.-Z. Jiang, Q. Wang, and Y. Zhang, Phys. Rev. D 87, 094014 (2013).
  * Jiang  _et al._ (2015) S.-Z. Jiang, Z.-L. Wei, Q.-S. Chen, and Q. Wang, Phys. Rev. D 92, 025014 (2015).
  * Wang and Wang (2000) X.-L. Wang and Q. Wang, Commun. Theor. Phys. 34, 519 (2000).
  * Wang  _et al._ (2000b) X.-L. Wang, Z.-M. Wang, and Q. Wang, Commun. Theor. Phys. 34, 683 (2000b).
  * Ren  _et al._ (2017) K. Ren, H.-F. Fu, and Q. Wang, Phys. Rev. D 95, 074012 (2017).
  * Chen  _et al._ (2020) Q.-S. Chen, H.-F. Fu, Y.-L. Ma, and Q. Wang, Phys. Rev. D 102, 034034 (2020), arXiv:2001.06418 [hep-ph] .
  * Wise (1992) M. B. Wise, Phys. Rev. D 45, R2188 (1992).
  * Yan  _et al._ (1992) T.-M. Yan, H.-Y. Cheng, C.-Y. Cheung, G.-L. Lin, Y. Lin, and H.-L. Yu, Phys. Rev. D 46, 1148 (1992), [Erratum: Phys.Rev.D 55, 5851 (1997)].
  * Nowak  _et al._ (1993) M. A. Nowak, M. Rho, and I. Zahed, Phys. Rev. D 48, 4370 (1993).
  * Casalbuoni  _et al._ (1993) R. Casalbuoni, A. Deandrea, N. Di Bartolomeo, R. Gatto, F. Feruglio, and G. Nardulli, Phys. Lett. B 299, 139 (1993), arXiv:hep-ph/9211248 .
  * Cheng  _et al._ (1994) H.-Y. Cheng, C.-Y. Cheung, G.-L. Lin, Y. Lin, T.-M. Yan, and H.-L. Yu, Phys. Rev. D 49, 2490 (1994), arXiv:hep-ph/9308283 .
  * Balk  _et al._ (1994) S. Balk, J. Korner, and D. Pirjol, Nucl. Phys. B 428, 499 (1994), arXiv:hep-ph/9307230 .
  * Di Bartolomeo  _et al._ (1995) N. Di Bartolomeo, R. Gatto, F. Feruglio, and G. Nardulli, Phys. Lett. B 347, 405 (1995), arXiv:hep-ph/9411210 .
  * Cheung and Hwang (2016) C.-Y. Cheung and C.-W. Hwang, Eur. Phys. J. C 76, 19 (2016), arXiv:1508.07686 [hep-ph] .
  * Alhakami (2020) M. H. Alhakami, Phys. Rev. D 101, 016001 (2020), arXiv:1910.04409 [hep-ph] .
  * Manohar and Wise (2000) A. V. Manohar and M. B. Wise, Cambridge Monogr. Part. Phys. Nucl. Phys., Cosmol. 10, 1 (2000).
  * Nowak  _et al._ (2004) M. A. Nowak, M. Praszalowicz, M. Sadzikowski, and J. Wasiluk, Phys. Rev. D 70, 031503 (2004).
  * Bardeen and Hill (1994) W. A. Bardeen and C. T. Hill, Phys. Rev. D 49, 409 (1994).
  * Zyla  _et al._ (2020) P. Zyla _et al._ (Particle Data Group), Prog. Theor. Exp. Phys. 2020, 083C01 (2020).
  * Pagels and Stokar (1979) H. Pagels and S. Stokar, Phys. Rev. D 20, 2947 (1979).
  * Dudal  _et al._ (2012) D. Dudal, O. Oliveira, and J. Rodriguez-Quintero, Phys. Rev. D 86, 105005 (2012).
  * Oliveira  _et al._ (2019) O. Oliveira, W. de Paula, T. Frederico, and J. de Melo, Eur. Phys. J. C 79, 116 (2019).
  * Fu  _et al._ (2017) H.-F. Fu, Q. Wang, and L. Jiang, Phys. Rev. D 96, 094023 (2017), arXiv:1706.03181 [hep-th] .
  * Tandy (1997) P. C. Tandy, Prog. Part. Nucl. Phys. 39, 117 (1997).