#  Strong Vertices of Doubly Heavy Spin-3/2 Baryons with Light Pseudoscalar
Mesons

A. R. Olamaei1,4, K. Azizi2,,3,4, S. Rostami2
1 Department of Physics, Jahrom University, Jahrom, P. O. Box 74137-66171,
Jahrom, Iran
2 Department of Physics, University of Tehran, North Karegar Ave. Tehran
14395-547, Iran
3 Department of Physics, Doǧuş University, Acibadem-Kadiköy, 34722 Istanbul,
Turkey
4 School of Particles and Accelerators, Institute for Research in Fundamental
Sciences (IPM),
P. O. Box 19395-5531, Tehran, Iran

The strong coupling constants are basic quantities that carry information of
the strong interactions among the baryon and meson multiplets as well as
information on the natures and internal structures of the involved hadrons.
These parameters enter to the transition matrix elements of various decays as
main inputs and they play key roles in analyses of the experimental data
including various hadrons. We determine the strong coupling constants among
the doubly heavy spin-$3/2$ baryons, $\Xi^{*}_{QQ^{\prime}}$ and
$\Omega^{*}_{QQ^{\prime}}$, and light pseudoscalar mesons, $\pi$, $K$ and
$\eta$, using the light-cone QCD. The values obtained for the strong coupling
constants under study may be used in construction of the strong potentials
among the doubly heavy spin-3/2 baryons and light pseudoscalar mesons.

##  1 Introduction

Doubly heavy baryons composed of two heavy quarks and a single light quark are
interesting objects, investigation of which can help us better understand the
non-perturbative nature of QCD. Their investigation provides a good framework
for understanding and predicting the spectrum of heavy baryons. Theoretical
studies on different aspects of these baryons may shed light on the
experimental searches of these states. These baryons have been in the focus of
various theoretical studies [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14,
15, 16, 17, 18, 19, 20, 21, 27, 28, 29, 30, 31, 32, 22, 23, 24, 25, 26]. Some
interesting progress on the OPE calculations of their lifetimes can be found
in [33, 34]. Despite their predictions decades ago via quark model, we have
very limited experimental knowledge on these baryons. Nevertheless, the year
2017 was special in this regard, because it was marked by the discovery of the
doubly charmed $\Xi^{++}_{cc}$ baryon by the LHCb Collaboration [35]. This
state was then confirmed in the decay channel
$\Xi^{++}_{cc}\rightarrow\Xi^{+}_{c}\pi^{+}$ [36]. Although, the charmed-
bottom state $\Xi^{0}_{bc}$ has been searched in the $\Xi^{0}_{bc}\rightarrow
D^{0}pK^{-}$ decay by the LHCb collaboration, no evidence for a signal is
found [37]. However, the detection of $\Xi^{++}_{cc}$ has raised hopes for the
discovery of other members of doubly heavy baryons. It is wonderful that the
LHC opens new horizons in the discovery of heavy baryons and provides the
possibility to investigate their electromagnetic, weak and strong decays.

Beyond the quark models [38, 39, 40], many effective models like lattice QCD
[41], Quark Spin Symmetry [42, 43], QCD Sum Rules [22, 23, 24, 25, 26] and
Light-Cone Sum Rules [19, 20, 44, 45, 46, 47, 55, 56, 57, 58] have been
proposed to describe masses, residues, lifetimes, strong coupling constants,
and other properties of doubly heavy baryons. Given that the expectations to
discover more doubly heavy baryons are growing, we are witnessing the rise of
further theoretical investigations in this respect.

Concerning the strong decays of doubly heavy baryons, the strong coupling
constants are building blocks. They are entered to the amplitudes of the
strong decays and the widths of the related decays are calculated in terms of
these constants. They can also be used to construct the strong potential
energy among the hadronic multiplets. These couplings appear in the low-energy
(long-distance) region of QCD, where the running coupling gets larger and
perturbation theory breaks down. Therefore, to calculate such coupling
constants, a non-perturbative approach should be used. One of the
comprehensive and reliable methods for evaluating the non-perturbative effects
is the light-cone QCD sum rule (LCSR) which have given many successful
descriptions of the hadronic properties so far. This method is a developed
version of the standard technique of SVZ sum rule [59], using the conventional
distribution amplitudes (DAs) of the on-sell states. The main difference
between SVZ sum rule and LCSR is that the latter employs the operator product
expansion (OPE) near the light-cone $x^{2}\approx 0$ instead of short distance
$x\approx 0$, and the corresponding matrix elements are parametrized in terms
of hadronic DAs which are classified according to their twists [60, 61, 62].

In this study, we calculate the strong coupling constants among the doubly
heavy spin-3/2 baryons, $\Xi^{*}_{QQ^{\prime}}$ and
$\Omega^{*}_{QQ^{\prime}}$, and light pseudoscalar mesons $\pi$, $K$ and
$\eta$. Here $Q$ and $Q^{\prime}$ both can be $b$ or $c$ quark. The paper is
organized as follows: In the next section, we derive the sum rules for the
strong coupling constants using the LCSR. Numerical analysis and results are
presented in section 3. The final section 4, is devoted to the summary and
conclusion.

##  2 Light-Cone Sum Rules for $B^{*}_{1}B^{*}_{2}{\cal P}$ Strong Couplings

In this section we aim to derive the sum rules for the strong coupling
constants among the doubly heavy spin-3/2 baryons ($\Xi^{*}_{QQ^{\prime}}$ and
$\Omega^{*}_{QQ^{\prime}}$) and light pseudoscalar mesons ($\pi$, $K$ and
$\eta$) using the method of LCSR. The quark content of the doubly heavy
spin-3/2 baryons are presented in Table 1.

Baryon | $q$ | $Q$ | $Q^{\prime}$
---|---|---|---
$\Xi^{*}_{QQ^{\prime}}$ |  $u$ or $d$ |  $b$ or $c$ |  $b$ or $c$
$\Omega^{*}_{QQ^{\prime}}$ | $s$ |  $b$ or $c$ |  $b$ or $c$
Table 1: The quark content of the doubly heavy spin-3/2 baryons.

The first step to calculate the strong couplings is to write the proper
correlation function (CF) in terms of doubly heavy baryons’ interpolating
current, $\eta_{\mu}$. That is

$\displaystyle\Pi_{\mu\nu}(p,q)=i\int d^{4}xe^{ipx}\left<{\cal P}(q)|{\cal
T}\left\\{\eta_{\mu}(x)\bar{\eta}_{\nu}(0)\right\\}|0\right>,$ (1)

where ${\cal P}(q)$ represents the pseudoscalar meson carrying the four-
momentum $q$, and $p$ represents the outgoing doubly heavy baryon four-
momentum. The above CF can be calculated in two ways:

  * $\bullet$

One by inserting the complete set of hadronic states with the same quantum
numbers of the corresponding doubly heavy baryons, which is called the
physical or phenomenological side of the CF. It is calculated in the timelike
region and contains observables like strong coupling constants.

  * $\bullet$

In the second way, it is calculated in the deep Euclidean spacelike region by
the help of OPE and in terms of DAs of the on-shell mesons and other QCD
degrees of freedom. It is called the theoretical or QCD side of the CF.

These two sides are matched via a dispersion integral which leads to the sum
rules for the corresponding coupling constants. To suppress the contributions
of the higher states and continuum, Borel transformation and continuum
subtraction procedures are applied. In the following, we explain each of these
steps in detail.

### Constructing Interpolating Current

The procedure to construct the doubly heavy spin-3/2 interpolating current is
to make a diquark structure having spin one and then attach the remaining
spin-1/2 quark to that to build a spin-3/2 structure which has to contain the
quantum numbers of the related baryon. Following the line of [23] and [47] one
can find the corresponding interpolating current as follows. The diquark
structure has the form

$\displaystyle\eta_{\text{diquark}}=q_{1}^{T}C\Gamma q_{2},$ (2)

where $T$ and $C$ represent the transposition and charge conjugation operators
respectively and
$\Gamma=I,\gamma_{5},\gamma_{\mu},\gamma_{\mu}\gamma_{5},\sigma_{\mu\nu}$. The
dependence of the spinors on the spacetime $x$ is omitted for now. Attaching
the third quark, the general form of the corresponding structure is
$[q_{1}^{T}C\Gamma q_{2}]\Gamma^{\prime}q_{3}$. Considering the color indices
($a$,$b$ and $c$), the possible structures have the following forms:
$\varepsilon_{abc}(Q^{aT}C\Gamma Q^{{}^{\prime}b})\Gamma^{\prime}q^{c}$,
$\varepsilon_{abc}(q^{aT}C\Gamma Q^{b})\Gamma^{\prime}Q^{{}^{\prime}c}$ and
$\varepsilon_{abc}(q^{aT}C\Gamma Q^{{}^{\prime}b})\Gamma^{\prime}Q^{c}$, where
$Q^{(^{\prime})}$ and $q$ are the heavy and light quark spinors respectively
and the antisymmetric Levi-Civita tensor, $\varepsilon_{abc}$, makes the
interpolating current color singlet.

The $\Gamma$ matrices can be determined by investigating the diquark part of
the interpolating current. In the first form, as the diquark structure has
spin 1, it must be symmetric under the exchange of the heavy quarks
$Q\leftrightarrow Q^{\prime}$. Transposing yields to

$[\epsilon_{abc}Q^{aT}C\Gamma Q^{\prime b}]^{T}=-\epsilon_{abc}Q^{\prime
bT}\Gamma^{T}C^{-1}Q^{a}=\epsilon_{abc}Q^{\prime
bT}C(C\Gamma^{T}C^{-1})Q^{a},$ (3)

where we have used the identities $C^{T}=C^{-1}$ and $C^{2}=-1$, and consider
the anticommutation of the spinor components since they are Grassmann numbers.
Since $C=i\gamma_{2}\gamma_{0}$, the quantity $C\Gamma^{T}C^{-1}$ would be

$C\Gamma^{T}C^{-1}=\begin{cases}\Gamma&\text{for
$\Gamma=1,\gamma_{5},\gamma_{\mu}\gamma_{5}$,}\\\ -\Gamma&\text{for
$\Gamma=\gamma_{\mu},\sigma_{\mu\nu}$.}\end{cases}$ (4)

Therefore, after switching the dummy indices in the LHS of Eq. (3) one gets

$[\epsilon_{abc}Q^{aT}C\Gamma Q^{\prime b}]^{T}=\pm\epsilon_{abc}Q^{\prime
aT}C\Gamma Q^{b},$ (5)

where the $+$ and $-$ signs are for $\Gamma=\gamma_{\mu},\sigma_{\mu\nu}$ and
$\Gamma=1,\gamma_{5},\gamma_{5}\gamma_{\mu}$ respectively. On the other hand,
since the RHS of the above equation has to be symmetric with respect to the
exchange of heavy quarks, we have

$[\epsilon_{abc}Q^{aT}C\Gamma Q^{\prime b}]^{T}=\pm\epsilon_{abc}Q^{aT}C\Gamma
Q^{\prime b},$ (6)

where the $+$ and $-$ signs are for $\Gamma=\gamma_{\mu},\sigma_{\mu\nu}$ and
$\Gamma=1,\gamma_{5},\gamma_{5}\gamma_{\mu}$ respectively. Moreover, we know
that the structure $\epsilon_{abc}Q^{aT}C\Gamma Q^{\prime b}$ is a number
($1\times 1$ matrix) and therefore the transposition leaves it unchanged. This
leads us to conclude that the only possible choice for the $\Gamma$ matrices
would be $\Gamma=\gamma_{\mu}$ or $\sigma_{\mu\nu}$.

The above mentioned symmetry of exchanging heavy quarks should be respected by
the two remaining structures, $\varepsilon_{abc}(q^{aT}C\Gamma
Q^{b})\Gamma^{\prime}Q^{{}^{\prime}c}$ and $\varepsilon_{abc}(q^{aT}C\Gamma
Q^{{}^{\prime}b})\Gamma^{\prime}Q^{c}$, properly which leads to the following
form:

$\displaystyle\varepsilon_{abc}\big{[}(q^{aT}C\Gamma
Q^{b})\Gamma^{\prime}Q^{{}^{\prime}c}+(q^{aT}C\Gamma
Q^{{}^{\prime}b})\Gamma^{\prime}Q^{c}\big{]},$ (7)

where $\Gamma=\gamma_{\mu}$ or $\sigma_{\mu\nu}$. Consequently, the
interpolating current can be written in two possible forms as

$\displaystyle\varepsilon_{abc}\Big{\\{}\big{(}Q^{aT}C\gamma_{\mu}Q^{{}^{\prime}b}\big{)}\Gamma^{\prime}_{1}q^{c}+\big{(}q^{aT}C\gamma_{\mu}Q^{b}\big{)}\Gamma^{\prime}_{1}Q^{{}^{\prime}c}+\big{(}q^{aT}C\gamma_{\mu}Q^{{}^{\prime}b}\big{)}\Gamma^{\prime}_{1}Q^{c}\Big{\\}},$
(8)

and

$\displaystyle\varepsilon_{abc}\Big{\\{}\big{(}Q^{aT}C\sigma_{\mu\nu}Q^{{}^{\prime}b}\big{)}\Gamma^{\prime}_{2}q^{c}+\big{(}q^{aT}C\sigma_{\mu\nu}Q^{b}\big{)}\Gamma^{\prime}_{2}Q^{{}^{\prime}c}+\big{(}q^{aT}C\sigma_{\mu\nu}Q^{{}^{\prime}b}\big{)}\Gamma^{\prime}_{2}Q^{c}\Big{\\}}.$
(9)

To determine $\Gamma^{\prime}_{1}$ and $\Gamma^{\prime}_{2}$, one should
consider Lorentz and parity symmetries. As Eqs. (8) and (9) must have the
Lorentz vector structure, $\Gamma^{\prime}_{1}=1$ or $\gamma_{5}$ and
$\Gamma^{\prime}_{2}=\gamma_{\nu}$ or $\gamma_{\nu}\gamma_{5}$. But parity
considerations exclude the $\gamma_{5}$ matrix and therefore
$\Gamma^{\prime}_{1}=1$ and $\Gamma^{\prime}_{2}=\gamma_{\nu}$. Moreover, Eq.
(9) won’t survive if one consider all three quarks the same, and therefore the
only possible choice comes from Eq. (8) as

$\displaystyle\eta_{\mu}(x)$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{3}}\varepsilon_{abc}\Big{\\{}\big{[}Q^{aT}(x)C\gamma_{\mu}Q^{{}^{\prime}b}(x)\big{]}q^{c}(x)+\big{[}q^{aT}(x)C\gamma_{\mu}Q^{b}(x)\big{]}Q^{{}^{\prime}c}(x)$
(10)
$\displaystyle+\big{[}q^{aT}(x)C\gamma_{\mu}Q^{{}^{\prime}b}(x)\big{]}Q^{c}(x)\Big{\\}}.$

Not that these currents couple not only to the positive-parity but to
negative-parity doubly-heavy baryons [48, 49, 50, 51, 52, 53]. Considering the
negative-parity partners bring uncertainties in the numerical results
discussed also in [54]. Here, we do not consider the contribution of negative-
parity baryons.

### Physical Side

At hadronic (low energy) level, first we insert the complete set of hadronic
states with the same quantum numbers of the corresponding initial and final
doubly heavy spin-3/2 baryons and then perform the Fourier transformation by
integrating over four-$x$. By isolating the ground state one gets

$\displaystyle\Pi^{\text{Phys.}}_{\mu\nu}(p,q)=\frac{\langle
0|\eta_{\mu}|B^{*}_{2}(p,r)\rangle\langle B^{*}_{2}(p,r){\cal
P}(q)|B^{*}_{1}(p+q,s)\rangle\langle
B^{*}_{1}(p+q,s)|\bar{\eta}_{\mu}|0\rangle}{(p^{2}-m_{2}^{2})[(p+q)^{2}-m_{1}^{2}]}+\cdots~{},$
(11)

where $B^{*}_{1}(p+q,s)$ and $B^{*}_{2}(p,r)$ are the incoming and outgoing
doubly heavy spin-3/2 baryons with the masses $m_{1}$ and $m_{2}$
respectively. $s$ and $r$ are their spins and dots represent the higher states
and continuum. The matrix element $\langle 0|\eta_{\mu}|B^{*}_{i}(p,s)\rangle$
is defined as

$\displaystyle\langle 0|\eta_{\mu}|B^{*}_{i}(p,s)\rangle$ $\displaystyle=$
$\displaystyle\lambda_{B^{*}_{i}}u_{\mu}(p,s),$ (12)

where $u_{\mu}(p,s)$ is the Rarita–Schwinger spinor and $\lambda_{B^{*}_{i}}$
is the residue for the baryon $B^{*}_{i}$. The matrix element $\langle
B^{*}_{2}(p,r){\cal P}(q)|B^{*}_{1}(p+q,s)\rangle$ can be determined using the
Lorentz and parity considerations as

$\displaystyle\langle B^{*}_{2}(p,r){\cal
P}(q)|B^{*}_{1}(p+q,s)\rangle=g_{B^{*}_{1}B^{*}_{2}{\cal
P}}\bar{u}_{\alpha}(p,r)\gamma_{5}u^{\alpha}(p+q,s),$ (13)

where $g_{B^{*}_{1}B^{*}_{2}{\cal P}}$ is the strong coupling constant of the
doubly heavy spin-3/2 baryons $B^{*}_{1}$ and $B^{*}_{2}$ with the light
pseudoscalar meson ${\cal P}$. After substituting the above matrix elements,
(12) and (13), in (11), since the initial and final baryons are unpolarized,
one needs to sum over their spins using the following completeness relation:

$\displaystyle\sum_{s}u_{\mu}(p,s)\bar{u}_{\nu}(p,s)=-({\hbox
to0.0pt{/\hss}p+m})\Bigg{(}g_{\mu\nu}-{1\over
3}\gamma_{\mu}\gamma_{\nu}-{2p_{\mu}p_{\nu}\over
3m^{2}}+{p_{\mu}\gamma_{\nu}-p_{\nu}\gamma_{\mu}\over 3m}\Bigg{)}~{},$ (14)

which may lead to the corresponding physical side of the CF. But here we face
two major problems: First, not all emerging Lorentz structures are
independent; and second, since the interpolating current $\eta_{\mu}$ also
couples to the spin-1/2 doubly heavy baryon states, there are some unwanted
contributions from them which must be removed properly. Imposing the condition
$\gamma^{\mu}\eta_{\mu}=0$, these contributions can be written as

$\displaystyle\langle
0|\eta_{\mu}|B(p,s=1/2)\rangle=A\Big{(}\gamma_{\mu}-{4\over
m_{\frac{1}{2}}}p_{\mu}\Big{)}u(p,s=1/2)~{}.$ (15)

To fix the above-mentioned problems, we re-order the Dirac matrices in a way
that help us eliminate the spin-1/2 states’ contributions easily. The ordering
we choose is $\gamma_{\mu}\hbox to0.0pt{/\hss}p\hbox
to0.0pt{/\hss}q\gamma_{\nu}\gamma_{5}$ which leads us to the final form of the
physical side

$\displaystyle\Pi^{\text{Phys.}}_{\mu\nu}(p,q)$ $\displaystyle=$
$\displaystyle{\lambda_{B_{1}^{*}}\lambda_{B_{2}^{*}}\over[(p+q)^{2}-m_{1}^{2})](p^{2}-m_{2}^{2})}\Big{\\{}\frac{2~{}g_{B^{*}_{1}B^{*}_{2}{\cal
P}}}{3m_{1}^{2}}q_{\mu}q_{\nu}\hbox to0.0pt{/\hss}p\hbox
to0.0pt{/\hss}q\gamma_{5}+\mbox{\rm structures }$ (16) beginning with
$\gamma_{\mu}$ and ending with $\gamma_{\nu}\gamma_{5}$, or terms that are
$\displaystyle\mbox{\rm proportional to $p_{\mu}$ or $(p+q)_{\nu}$}+\mbox{\rm
other structures}\Big{\\}}$ $\displaystyle+\int
ds_{1}ds_{2}e^{-(s_{1}+s_{2})/2M^{2}}\rho^{\text{Phys.}}(s_{1},s_{2}),$

where $\rho^{\text{Phys.}}(s_{1},s_{2})$ represents the spectral density for
the higher states and continuum. One should note that there are many
structures like $g_{\mu\nu}\hbox to0.0pt{/\hss}p\hbox
to0.0pt{/\hss}q\gamma_{5}$ that emerge in the QCD side of the CF which are not
shown in the Eq. (16) expansion and can be considered to calculate the strong
coupling constant $g_{B^{*}_{1}B^{*}_{2}{\cal P}}$. We select the structure
with large number of momenta, $q_{\mu}q_{\nu}\hbox to0.0pt{/\hss}p\hbox
to0.0pt{/\hss}q\gamma_{5}$, which leads to more stable and reliable results
and it is free of the unwanted doubly heavy spin-1/2 contributions.

Now, to suppress the contributions of the higher states and continuum, we
perform the double Borel transformation with respect to the squared momenta
$p_{1}^{2}=(p+q)^{2}$ and $p_{2}^{2}=p^{2}$, which leads to

$\displaystyle{\cal B}_{p_{1}}(M_{1}^{2}){\cal
B}_{p_{2}}(M_{2}^{2})\Pi^{\text{Phys.}}_{\mu\nu}(p,q)$ $\displaystyle\equiv$
$\displaystyle\Pi^{\text{Phys.}}_{\mu\nu}(M^{2})$ $\displaystyle=$
$\displaystyle\frac{2~{}g_{B^{*}_{1}B^{*}_{2}{\cal
P}}}{3m_{1}^{2}}\lambda_{B_{1}}\lambda_{B_{2}}e^{-m_{1}^{2}/M_{1}^{2}}e^{-m_{2}^{2}/M_{2}^{2}}q_{\mu}q_{\nu}\hbox
to0.0pt{/\hss}p\hbox to0.0pt{/\hss}q\gamma_{5}~{}+\cdots~{},$

where dots represent suppressed higher states and continuum contributions,
$M_{1}^{2}$ and $M_{2}^{2}$ are the Borel parameters and
$M^{2}=M^{2}_{1}M^{2}_{2}/(M^{2}_{1}+M^{2}_{2})$. The Borel parameters are
chosen to be equal since the mass of the initial and final baryons are the
same, therefore $M^{2}_{1}=M^{2}_{2}=2M^{2}$. After matching the physical and
QCD sides, we will perform the continuum subtraction supplied by quark-hadron
duality assumption.

### QCD Side

After evaluating the physical side of the CF, the next stage is to calculate
the QCD side in the deep Euclidean spacelike region, where
$-(p+q)^{2}\rightarrow\infty$ and $-p^{2}\rightarrow\infty$. The CF in this
way is calculated in terms of QCD degrees of freedom as well as non-local
matrix elements of pseudoscalar mesons expressed in terms of the DAs of
different twists [63, 64, 65].

To proceed, according to Eq. (16), we choose the relevant structure
$q_{\mu}q_{\nu}\hbox to0.0pt{/\hss}p\hbox to0.0pt{/\hss}q\gamma_{5}$ from the
QCD side and express the CF as

$\displaystyle\Pi^{\text{QCD}}_{\mu\nu}(p,q)=\Pi(p,q)q_{\mu}q_{\nu}\hbox
to0.0pt{/\hss}p\hbox to0.0pt{/\hss}q\gamma_{5},$ (18)

where $\Pi(p,q)$ is an invariant function of $(p+q)^{2}$ and $p^{2}$. To
calculate $\Pi(p,q)$, we insert the interpolating current (10) into the CF
(1), and using the Wick theorem we find the QCD side of the CF as:

$\displaystyle\Big{(}\Pi^{\text{QCD}}_{\mu\nu}\Big{)}_{\rho\sigma}(p,q)$
$\displaystyle=$
$\displaystyle\frac{i}{3}\epsilon_{abc}\epsilon_{a^{\prime}b^{\prime}c^{\prime}}\int
d^{4}xe^{iq.x}\langle{\cal
P}(q)|\bar{q}^{c^{\prime}}_{\alpha}(0)q^{c}_{\beta}(x)|0\rangle$ (19)
$\displaystyle\times\Bigg{\\{}\delta_{\rho\alpha}\delta_{\beta\sigma}\text{Tr}\Big{[}\tilde{S}^{aa^{\prime}}_{Q}(x)\gamma_{\mu}S^{bb^{\prime}}_{Q^{\prime}}(x)\gamma_{\nu}\Big{]}+\delta_{\alpha\rho}\Big{(}\gamma_{\nu}\tilde{S}^{aa^{\prime}}_{Q}(x)\gamma_{\mu}S^{bb^{\prime}}_{Q^{\prime}}(x)\Big{)}_{\beta\sigma}$
$\displaystyle+\delta_{\alpha\rho}\Big{(}\gamma_{\nu}\tilde{S}^{bb^{\prime}}_{Q^{\prime}}(x)\gamma_{\mu}S^{aa^{\prime}}_{Q}(x)\Big{)}_{\beta\sigma}+\delta_{\beta\sigma}\Big{(}S^{bb^{\prime}}_{Q^{\prime}}(x)\gamma_{\nu}\tilde{S}^{aa^{\prime}}_{Q}(x)\gamma_{\mu}\Big{)}_{\rho\alpha}$
$\displaystyle-\delta_{\beta\sigma}\Big{(}S^{aa^{\prime}}_{Q}(x)\gamma_{\nu}\tilde{S}^{bb^{\prime}}_{Q^{\prime}}(x)\gamma_{\mu}\Big{)}_{\rho\alpha}+\Big{(}\gamma_{\nu}\tilde{S}^{aa^{\prime}}_{Q}(x)\gamma_{\mu}\Big{)}_{\beta\alpha}\Big{(}S^{bb^{\prime}}_{Q^{\prime}}(x)\Big{)}_{\rho\sigma}$
$\displaystyle-\Big{(}C\gamma_{\mu}S^{aa^{\prime}}_{Q}(x)\Big{)}_{\alpha\sigma}\Big{(}S^{bb^{\prime}}_{Q^{\prime}}(x)\gamma_{\nu}C\Big{)}_{\rho\beta}-\Big{(}C\gamma_{\mu}S^{bb^{\prime}}_{Q^{\prime}}(x)\Big{)}_{\alpha\sigma}\Big{(}S^{aa^{\prime}}_{Q}(x)\gamma_{\nu}C\Big{)}_{\rho\beta}$
$\displaystyle+\Big{(}\gamma_{\nu}\tilde{S}^{bb^{\prime}}_{Q^{\prime}}(x)\gamma_{\mu}\Big{)}_{\beta\alpha}\Big{(}S^{aa^{\prime}}_{Q}(x)\Big{)}_{\rho\sigma}\Bigg{\\}},$

where $S^{aa^{\prime}}_{Q}(x)$ is the propagator of the heavy quark $Q$ with
color indices $a$ and $a^{\prime}$, and $\tilde{S}=CS^{T}C$. $\mu$ and $\nu$
are Minkowski and $\rho$ and $\sigma$ are Dirac indices respectively. The
propagators in the above equation contain both the perturbative and the non-
perturbative contributions. As we previously mentioned, the non-local matrix
elements $\langle{\cal
P}(q)|\bar{q}^{c^{\prime}}_{\alpha}(0)q^{c}_{\beta}(x)|0\rangle$, are written
in terms of DAs of the corresponding light pseudoscalar meson ${\cal P}(q)$.

The next stage is to insert the explicit expression of the heavy quark
propagator in Eq. (19) as [68]:

$\displaystyle S_{Q}^{aa^{\prime}}(x)$ $\displaystyle=$
$\displaystyle{m_{Q}^{2}\over
4\pi^{2}}{K_{1}(m_{Q}\sqrt{-x^{2}})\over\sqrt{-x^{2}}}\delta^{aa^{\prime}}-i{m_{Q}^{2}\hbox
to0.0pt{/\hss}{x}\over
4\pi^{2}x^{2}}K_{2}(m_{Q}\sqrt{-x^{2}})\delta^{aa^{\prime}}$ (20)
$\displaystyle-
ig_{s}\int{d^{4}k\over(2\pi)^{4}}e^{-ikx}\int_{0}^{1}du\Bigg{[}{\hbox
to0.0pt{/\hss}k+m_{Q}\over
2(m_{Q}^{2}-k^{2})^{2}}\sigma^{\lambda\tau}G_{\lambda\tau}^{aa^{\prime}}(ux)$
$\displaystyle+{u\over
m_{Q}^{2}-k^{2}}x^{\lambda}\gamma^{\tau}G_{\lambda\tau}^{aa^{\prime}}(ux)\Bigg{]}+\cdots,$

where $m_{Q}$ is the heavy quark mass, $K_{1}$ and $K_{2}$ are the modified
Bessel functions of the second kind and $G_{\lambda\tau}^{aa^{\prime}}$ is the
gluon field strength tensor. It is defined as

$\displaystyle G^{aa^{\prime}}_{\lambda\tau}\equiv
G^{A}_{\lambda\tau}t^{aa^{\prime}}_{A},$ (21)

with $\lambda$ and $\tau$ being the Minkowski indices.
$t^{aa^{\prime}}_{A}=\lambda^{aa^{\prime}}_{A}/2$ where $\lambda_{A}$ are the
Gell-Mann matrices with $A=1,\cdots,8$ and $a,a^{\prime}$ the color indices.
The free propagator contribution is determined by the first two terms and the
rest which $\sim G_{\lambda\tau}^{aa^{\prime}}$ are due to the interaction
with the gluon field.

Inserting Eq. (20) into (19) leads to several contributions. The first one is
due to replacing both the heavy quark propagators with their free part:

$\displaystyle S_{Q}^{(\text{pert.})}(x)$ $\displaystyle=$
$\displaystyle{m_{Q}^{2}\over
4\pi^{2}}{K_{1}(m_{Q}\sqrt{-x^{2}})\over\sqrt{-x^{2}}}-i{m_{Q}^{2}\hbox
to0.0pt{/\hss}{x}\over 4\pi^{2}x^{2}}K_{2}(m_{Q}\sqrt{-x^{2}}).$ (22)

The DAS in this case are two-particle DAs. Replacing one heavy quark
propagator (say $S^{aa^{\prime}}_{Q}$) with its gluonic part

$\displaystyle S^{aa^{\prime}(\text{non-p.})}_{Q}(x)$ $\displaystyle=$
$\displaystyle-
ig_{s}\int{d^{4}k\over(2\pi)^{4}}e^{-ikx}\int_{0}^{1}duG^{aa^{\prime}}_{\lambda\tau}(ux)\Delta^{\lambda\tau}_{Q}(x),$
(23)

where

$\displaystyle\Delta^{\lambda\tau}_{Q}(x)$ $\displaystyle=$
$\displaystyle\dfrac{1}{2(m_{Q}^{2}-k^{2})^{2}}\Big{[}(\hbox
to0.0pt{/\hss}k+m_{Q})\sigma^{\lambda\tau}+2u(m_{Q}^{2}-k^{2})x^{\lambda}\gamma^{\tau}\Big{]},$
(24)

and the other with its free part, leads to the one gluon exchange between the
heavy quark $Q$ and the light pseudoscalar meson ${\cal P}$. The non-local
matrix elements of this contribution can be calculated in terms of three-
particle DAs of meson $\cal P$. Replacing both heavy quark propagators with
their gluonic parts involves ${\cal P}$ meson four-particle DAs which are not
yet determined and we ignore them in the present work. Instead, we consider
the two-gluon condensate contributions in this study.

The non-local matrix elements can be expanded using proper Fierz identities
like

$\displaystyle\bar{q}_{\alpha}^{c^{\prime}}q_{\beta}^{c}\to-\frac{1}{12}(\Gamma_{J})_{\beta\alpha}\delta^{cc^{\prime}}\bar{q}\Gamma^{J}q,$
(25)

where $\Gamma^{J}=\mathbf{1,\ }\gamma_{5},\ \gamma_{\mu},\
i\gamma_{5}\gamma_{\mu},\ \sigma_{\mu\nu}/\sqrt{2}$. It puts them in a form to
be determined in terms of the corresponding DAs of different twists which can
be found in Refs. [63, 64, 65].

Putting all things together, one can calculate different contributions to the
corresponding strong coupling. The leading order contribution corresponding to
no gluon exchange which can be obtained by replacing both heavy quark
propagators with their free part is as follows:

$\displaystyle\Big{(}\Pi^{\text{QCD(0)}}_{\mu\nu}\Big{)}_{\rho\sigma}(p,q)$
$\displaystyle=$ $\displaystyle\frac{i}{6}\int d^{4}xe^{iq.x}\langle{\cal
P}(q)|\bar{q}(0)\Gamma^{J}q(x)|0\rangle\Bigg{\\{}\text{Tr}\Big{[}\tilde{S}^{(\text{pert.})}_{Q}(x)\gamma_{\mu}S^{(\text{pert.})}_{Q^{\prime}}(x)\gamma_{\nu}\Big{]}\Big{(}\Gamma_{J}\Big{)}_{\rho\sigma}$
(26) $\displaystyle+$
$\displaystyle\text{Tr}\Big{[}\Gamma_{J}\gamma_{\nu}\tilde{S}^{(\text{pert.})}_{Q}(x)\gamma_{\mu}\Big{]}\Big{(}S^{(\text{pert.})}_{Q^{\prime}}(x)\Big{)}_{\rho\sigma}+\text{Tr}\Big{[}\Gamma_{J}\gamma_{\nu}\tilde{S}^{(\text{pert.})}_{Q^{\prime}}(x)\gamma_{\mu}\Big{]}\Big{(}S^{(\text{pert.})}_{Q}(x)\Big{)}_{\rho\sigma}$
$\displaystyle+$
$\displaystyle\Big{(}\Gamma_{J}\gamma_{\nu}\tilde{S}^{(\text{pert.})}_{Q}(x)\gamma_{\mu}S^{(\text{pert.})}_{Q^{\prime}}(x)\Big{)}_{\rho\sigma}+\Big{(}\Gamma_{J}\gamma_{\nu}\tilde{S}^{(\text{pert.})}_{Q^{\prime}}(x)\gamma_{\mu}S^{(\text{pert.})}_{Q}(x)\Big{)}_{\rho\sigma}$
$\displaystyle+$
$\displaystyle\Big{(}S^{(\text{pert.})}_{Q^{\prime}}(x)\gamma_{\nu}\tilde{S}^{(\text{pert.})}_{Q}(x)\gamma_{\mu}\Gamma_{J}\Big{)}_{\rho\sigma}-\Big{(}S^{(\text{pert.})}_{Q^{\prime}}(x)\gamma_{\nu}\tilde{\Gamma}_{J}\gamma_{\mu}S^{(\text{pert.})}_{Q}(x)\Big{)}_{\rho\sigma}$
$\displaystyle-$
$\displaystyle\Big{(}S^{(\text{pert.})}_{Q}(x)\gamma_{\nu}\tilde{S}^{(\text{pert.})}_{Q^{\prime}}(x)\gamma_{\mu}\Gamma_{J}\Big{)}_{\rho\sigma}-\Big{(}S^{(\text{pert.})}_{Q}(x)\gamma_{\nu}\tilde{\Gamma}_{J}\gamma_{\mu}S^{(\text{pert.})}_{Q^{\prime}}(x)\Big{)}_{\rho\sigma}\Bigg{\\}},$

where the superscript $(0)$ indicates no gluon exchange. The contribution of
the exchange of one gluon between the heavy quark $Q$ and the light
pseudoscalar meson ${\cal P}$ is obtained as

$\displaystyle\Big{(}\Pi^{\text{QCD(1)}}_{\mu\nu}\Big{)}_{\rho\sigma}(p,q)$
$\displaystyle=$
$\displaystyle\frac{-ig_{s}}{96}\int\frac{d^{4}k}{(2\pi)^{2}}e^{-ik.x}\int_{0}^{1}du\langle{\cal
P}(q)|\bar{q}(0)\Gamma^{J}G_{\lambda\tau}q(x)|0\rangle$ (27)
$\displaystyle\times$
$\displaystyle\Bigg{\\{}\text{Tr}\Big{[}\tilde{\Delta}^{\lambda\tau}_{Q}(x)\gamma_{\mu}S^{(\text{pert.})}_{Q^{\prime}}(x)\gamma_{\nu}\Big{]}\Big{(}\Gamma_{J}\Big{)}_{\rho\sigma}+\text{Tr}\Big{[}\gamma_{\nu}\tilde{\Delta}^{\lambda\tau}_{Q}(x)\gamma_{\mu}\Gamma_{J}\Big{]}\Big{(}S^{(\text{pert.})}_{Q^{\prime}}(x)\Big{)}_{\rho\sigma}$
$\displaystyle+$
$\displaystyle\text{Tr}\Big{[}\Gamma_{J}\gamma_{\nu}\tilde{S}^{(\text{pert.})}_{Q^{\prime}}(x)\gamma_{\mu}\Big{]}\Big{(}\Delta^{\lambda\tau}_{Q}(x)\Big{)}_{\rho\sigma}+\Big{(}\Gamma_{J}\gamma_{\nu}\tilde{\Delta}^{\lambda\tau}_{Q}(x)\gamma_{\mu}S^{(\text{pert.})}_{Q^{\prime}}(x)\Big{)}_{\rho\sigma}$
$\displaystyle+$
$\displaystyle\Big{(}\Gamma_{J}\gamma_{\nu}\tilde{S}^{(\text{pert.})}_{Q^{\prime}}(x)\gamma_{\mu}\Delta^{\lambda\tau}_{Q}(x)\Big{)}_{\rho\sigma}+\Big{(}S^{(\text{pert.})}_{Q^{\prime}}(x)\gamma_{\nu}\tilde{\Delta}^{\lambda\tau}_{Q}(x)\gamma_{\mu}\Gamma_{J}\Big{)}_{\rho\sigma}$
$\displaystyle-$
$\displaystyle\Big{(}S^{(\text{pert.})}_{Q^{\prime}}(x)\gamma_{\nu}\tilde{\Gamma}_{J}\gamma_{\mu}\Delta^{\lambda\tau}_{Q}(x)\Big{)}_{\rho\sigma}-\Big{(}\Delta^{\lambda\tau}_{Q}(x)\gamma_{\nu}\tilde{S}^{(\text{pert.})}_{Q^{\prime}}(x)\gamma_{\mu}\Gamma_{J}\Big{)}_{\rho\sigma}$
$\displaystyle-$
$\displaystyle\Big{(}\Delta^{\lambda\tau}_{Q}(x)\gamma_{\nu}\tilde{\Gamma}_{J}\gamma_{\mu}S^{(\text{pert.})}_{Q^{\prime}}(x)\Big{)}_{\rho\sigma}\Bigg{\\}},$

where the superscript $(1)$ indicates one gluon exchange. By exchanging $Q$
and $Q^{\prime}$ in the above equation one can simply find the contribution of
one gluon exchange between the heavy quark $Q^{\prime}$ and the light
pseudoscalar meson ${\cal P}$.

The general configurations appear in the calculation of the QCD side of the CF
(26) and (27) have the form

$\displaystyle T_{[~{}~{},\mu,\mu\nu,...]}(p,q)$ $\displaystyle=$
$\displaystyle i\int d^{4}x\int_{0}^{1}dv\int{\cal D}\alpha
e^{ip.x}\big{(}x^{2}\big{)}^{n}[e^{i(\alpha_{q}+v\alpha_{g})q.x}\mathcal{G}(\alpha_{i}),e^{iq.x}f(u)]$
(28) $\displaystyle\times$
$\displaystyle[1,x_{\mu},x_{\mu}x_{\nu},...]K_{n_{1}}(m_{1}\sqrt{-x^{2}})K_{n_{2}}(m_{2}\sqrt{-x^{2}}).$

The expressions in the brackets on the RHS correspond to different
configurations which might arise in the calculation and

$\int\mathcal{D}\alpha=\int_{0}^{1}d\alpha_{q}\int_{0}^{1}d\alpha_{\bar{q}}\int_{0}^{1}d\alpha_{g}\delta(1-\alpha_{q}-\alpha_{\bar{q}}-\alpha_{g}).$
(29)

On the LHS, the blank subscript indicates no $x_{\mu}$ in the corresponding
configuration. There are several representations for the modified Bessel
function of the second kind and we use the cosine representation as

$K_{n}(m_{Q}\sqrt{-x^{2}})=\frac{\Gamma(n+1/2)~{}2^{n}}{\sqrt{\pi}m_{Q}^{n}}\int_{0}^{\infty}dt~{}\cos(m_{Q}t)\frac{(\sqrt{-x^{2}})^{n}}{(t^{2}-x^{2})^{n+1/2}}.$
(30)

It is shown that choosing this representation increases the radius of
convergence of the Borel transformed CF [44]. To perform the Fourier integral
over $x$, we write the $x$ configurations in the exponential representation as

$\displaystyle(x^{2})^{n}$ $\displaystyle=$
$\displaystyle(-1)^{n}\frac{\partial^{n}}{\partial\beta^{n}}\big{(}e^{-\beta
x^{2}}\big{)}\arrowvert_{\beta=0},$ $\displaystyle x_{\mu}e^{iP.x}$
$\displaystyle=$ $\displaystyle(-i)\frac{\partial}{\partial P^{\mu}}e^{iP.x}.$
(31)

### Borel Transformation and Continuum Subtraction

Performing the Fourier transformation, gives us the version of CF that has to
be Borel transformed to suppress the divergences arise due to the dispersion
integral. Having two independent momenta $(p+q)$ and $p$, we use the double
Borel transformation with respect to the square of these momenta as

${\cal B}_{p_{1}}(M_{1}^{2}){\cal
B}_{p_{2}}(M_{2}^{2})e^{b(p+uq)^{2}}=M^{2}\delta(b+\frac{1}{M^{2}})\delta(u_{0}-u)e^{\frac{-q^{2}}{M_{1}^{2}+M_{2}^{2}}},$
(32)

where $u_{0}=M_{1}^{2}/(M_{1}^{2}+M_{2}^{2})$. To be specific, we select the
following configuration:

$\displaystyle T_{\mu\nu}(p,q)$ $\displaystyle=$ $\displaystyle i\int
d^{4}x\int_{0}^{1}dv\int{\cal D}\alpha
e^{i[p+(\alpha_{q}+v\alpha_{g})q].x}\mathcal{G}(\alpha_{i})\big{(}x^{2}\big{)}^{n}$
(33) $\displaystyle\times$ $\displaystyle
x_{\mu}x_{\nu}K_{n_{1}}(m_{Q_{1}}\sqrt{-x^{2}})K_{n_{2}}(m_{Q_{2}}\sqrt{-x^{2}}),$

where after Fourier and Borel transformation leads to

$\displaystyle T_{\mu\nu}(M^{2})$ $\displaystyle=$
$\displaystyle\frac{i\pi^{2}2^{4-n_{1}-n_{2}}e^{\frac{-q^{2}}{M_{1}^{2}+M_{2}^{2}}}}{M^{2}m_{Q_{1}}^{2n_{1}}m_{Q_{2}}^{2n_{2}}}\int\mathcal{D}\alpha\int_{0}^{1}dv\int_{0}^{1}dz\frac{\partial^{n}}{\partial\beta^{n}}e^{-\frac{m_{Q_{1}}^{2}\bar{z}+m_{Q_{2}}^{2}z}{z\bar{z}(M^{2}-4\beta)}}z^{n_{1}-1}\bar{z}^{n_{2}-1}$
(34) $\displaystyle\times$
$\displaystyle(M^{2}-4\beta)^{n_{1}+n_{2}-1}\delta[u_{0}-(\alpha_{q}+v\alpha_{g})]\Big{[}p_{\mu}p_{\nu}+(v\alpha_{g}+\alpha_{q})(p_{\mu}q_{\nu}+q_{\mu}p_{\nu})$
$\displaystyle+(v\alpha_{g}+\alpha_{q})^{2}q_{\mu}q_{\nu}+\frac{M^{2}}{2}g_{\mu\nu}\Big{]}.$

In Ref. [44] one can find the details of calculations for hadrons containing
different numbers of heavy quarks (zero to five).

The dispersion integral that matches the physical and QCD sides of the CF
contains contributions from both the ground state and also the excited and
continuum ones which is calculated in Eq. (16) as a double dispersion integral
over $\rho^{\text{Phys.}}(s_{1},s_{2})$. To suppress the latter, one has to
perform a proper version of subtraction which enhances the contribution of the
ground state as well. There are several versions of subtraction. Using quark-
hadron duality one can approximate $\rho^{\text{Phys.}}(s_{1},s_{2})$, the
physical spectral density, with its theoretical counterpart,
$\rho^{\text{QCD}}(s_{1},s_{2})$. For more details see also Refs. [66, 67, 69,
70, 68, 71].

To proceed, for the generic factor $(M^{2})^{N}e^{-m^{2}/M^{2}}$, the
replacement [72]

$\left(M^{2}\right)^{N}e^{-m^{2}/M^{2}}\to\frac{1}{\Gamma(N)}\int_{m^{2}}^{s_{0}}dse^{-s/M^{2}}\left(s-m^{2}\right)^{N-1},$
(35)

is used for $N>0$ and it is left unchanged for $N<0$, where the energy
threshold for higher states and continuum is considered as $\sqrt{s_{0}}$.
Then we restrict the boundaries of $z$ integral in a way to more suppress the
unwanted excited states and continuum. Solving the equation

$\displaystyle
e^{-\frac{m_{Q_{1}}^{2}\bar{z}+m_{Q_{2}}^{2}z}{M^{2}z\bar{z}}}=e^{-s_{0}/M^{2}},$
(36)

for $z$ gives us the proper boundaries as follows:

$\displaystyle
z_{\text{max}(\text{min})}=\frac{1}{2s_{0}}\Big{[}(s_{0}+m_{Q_{1}}^{2}-m_{Q_{2}}^{2})+(-)\sqrt{(s0+m_{Q_{1}}^{2}-m_{Q_{2}}^{2})^{2}-4m_{Q_{1}}^{2}s_{0}}\Big{]}.$
(37)

Replacing them in the $z$ integral as

$\displaystyle\int_{0}^{1}dz\rightarrow\int_{z_{\text{min}}}^{z_{\text{max}}}dz,$
(38)

enhances the ground state and suppresses the higher states and continuum as
much as possible.

Finally, selecting the coefficient of the corresponding structure
$q_{\mu}q_{\nu}\hbox to0.0pt{/\hss}p\hbox to0.0pt{/\hss}q\gamma_{5}$, gives us
the invariant function to be used to find the corresponding strong coupling
constants. The Borel transformed subtracted version of the invariant function
for the vertex $B^{*}_{1}B^{*}_{2}{\cal P}$ can be written as

$\displaystyle\Pi_{B^{*}_{1}B^{*}_{2}{\cal
P}}(M^{2},s_{0})=T_{B^{*}_{1}B^{*}_{2}{\cal
P}}(M^{2},s_{0})+T_{B^{*}_{1}B^{*}_{2}{\cal P}}^{GG}(M^{2},s_{0}),$ (39)

where the above functions for the vertex $\Xi^{*}_{bb}\Xi^{*}_{bb}\pi^{0}$, as
an example, are given as

$\displaystyle T_{\Xi^{*}_{bb}\Xi^{*}_{bb}\pi^{0}}(M^{2},s_{0})$
$\displaystyle=$
$\displaystyle\frac{e^{-\frac{q^{2}}{M_{1}^{2}+M_{2}^{2}}}}{9\sqrt{2}\pi^{2}}\Bigg{\\{}6u_{0}f_{\pi}m_{b}m_{\pi}^{2}i_{1}({\cal
V}^{\parallel}(\alpha_{i}),1)\Big{[}\zeta_{-1,0}(m_{b})-\zeta_{0,0}(m_{b})\Big{]}$
$\displaystyle+$ $\displaystyle 6u_{0}f_{\pi}m_{b}m_{\pi}^{2}i_{1}({\cal
V}^{\perp}(\alpha_{i}),1)\Big{[}\zeta_{-1,0}(m_{b})-\zeta_{0,0}(m_{b})\Big{]}$
$\displaystyle+$ $\displaystyle\mu_{\pi}\Bigg{(}-3\Big{[}i_{2}(\alpha_{g}{\cal
T}(\alpha_{i}),v)-2i_{2}(\alpha_{g}{\cal
T}(\alpha_{i}),v^{2})+i_{2}(\alpha_{q}{\cal T}(\alpha_{i}),1)$
$\displaystyle-$ $\displaystyle 2i_{2}(\alpha_{q}{\cal
T}(\alpha_{i}),v)\Big{]}\Big{[}\tilde{\zeta}^{(1)}_{0,0}(m_{b};s_{0},M^{2})-\tilde{\zeta}^{(1)}_{1,0}(m_{b};s_{0},M^{2})\Big{]}$
$\displaystyle+$
$\displaystyle(-1+\rho_{\pi}^{2})u_{0}^{2}\phi_{\sigma}(u_{0})\Big{[}-\tilde{\zeta}^{(1)}_{1,0}(m_{b};s_{0},M^{2})+\tilde{\zeta}^{(1)}_{2,0}(m_{b};s_{0},M^{2})\Big{]}\Bigg{)}\Bigg{\\}},$

and

$\displaystyle T_{\Xi^{*}_{bb}\Xi^{*}_{bb}\pi^{0}}^{GG}(M^{2},s_{0})$
$\displaystyle=$ $\displaystyle g_{s}^{2}\langle
GG\rangle\frac{e^{-\frac{q^{2}}{M_{1}^{2}+M_{2}^{2}}}}{324\sqrt{2}M^{6}\pi^{2}}\Bigg{\\{}6u_{0}f_{\pi}m_{b}^{3}m_{\pi}^{2}i_{1}({\cal
V}^{\parallel}(\alpha_{i}),1)\Big{[}\zeta_{-2,0}(m_{b})+\zeta_{-1,-1}(m_{b})\Big{]}$
(41) $\displaystyle+$ $\displaystyle
6u_{0}f_{\pi}m_{b}^{3}m_{\pi}^{2}i_{1}({\cal
V}^{\perp}(\alpha_{i}),1)\Big{[}\zeta_{-2,0}(m_{b})+\zeta_{-1,-1}(m_{b})\Big{]}$
$\displaystyle+$ $\displaystyle\mu_{\pi}\Bigg{(}-3\Big{[}i_{2}(\alpha_{g}{\cal
T}(\alpha_{i}),v)-2i_{2}(\alpha_{g}{\cal
T}(\alpha_{i}),v^{2})+i_{2}(\alpha_{q}{\cal T}(\alpha_{i}),1)$
$\displaystyle-$ $\displaystyle 2i_{2}(\alpha_{q}{\cal
T}(\alpha_{i}),v)\Big{]}\Big{\\{}m_{b}^{2}\Big{[}\tilde{\zeta}^{(1)}_{-1,0}(m_{b};s_{0},M^{2})+\tilde{\zeta}^{(1)}_{0,-1}(m_{b};s_{0},M^{2})\Big{]}$
$\displaystyle+$
$\displaystyle\tilde{\zeta}^{(2)}_{0,0}(m_{b};s_{0},M^{2})\Big{\\}}-2(-1+\rho_{\pi}^{2})u_{0}^{2}\phi_{\sigma}(u_{0})\Big{[}m_{b}^{2}\tilde{\zeta}^{(1)}_{-1,0}(m_{b};s_{0},M^{2})$
$\displaystyle+$ $\displaystyle
m_{b}^{2}\tilde{\zeta}^{(1)}_{0,-1}(m_{b};s_{0},M^{2})+2\tilde{\zeta}^{(2)}_{0,0}(m_{b};s_{0},M^{2})-m_{b}^{2}\tilde{\zeta}^{(1)}_{0,0}(m_{b};s_{0},M^{2})$
$\displaystyle-$ $\displaystyle
m_{b}^{2}\tilde{\zeta}^{(1)}_{1,-1}(m_{b};s_{0},M^{2})-2\tilde{\zeta}^{(2)}_{1,0}(m_{b};s_{0},M^{2})\Big{]}\Bigg{)}\Bigg{\\}}.$

Here the function $\zeta_{m,n}(m_{1},m_{2})$ is defined as

$\displaystyle\zeta_{m,n}(m_{1},m_{2})=\int_{z_{\text{min}}}^{z_{\text{max}}}dzz^{m}\bar{z}^{n}e^{-\frac{m_{1}^{2}}{M^{2}z}-\frac{m_{2}^{2}}{M^{2}\bar{z}}}~{},$
(42)

where $m$ and $n$ are integers and
$\zeta_{m,n}(m_{1},m_{1})=\zeta_{m,n}(m_{1})$. Also
$\tilde{\zeta}^{(N)}_{m,n}(m_{1},m_{1};s_{0},M^{2})$ is defined as follows:

$\displaystyle\tilde{\zeta}^{(N)}_{m,n}(m_{1},m_{2};s_{0},M^{2})=\frac{1}{\Gamma(N)}\int_{z_{\text{min}}}^{z_{\text{max}}}dz\int_{\frac{m_{1}^{2}}{z}+\frac{m_{2}^{2}}{\bar{z}}}^{s_{0}}dse^{-s/M^{2}}z^{m}\bar{z}^{n}\Big{(}s-\frac{m_{1}^{2}}{z}-\frac{m_{2}^{2}}{\bar{z}}\Big{)}^{N-1}~{},$
(43)

where
$\tilde{\zeta}^{(N)}_{m,n}(m_{1},m_{1};s_{0},M^{2})=\tilde{\zeta}^{(N)}_{m,n}(m_{1};s_{0},M^{2})$.
The functions $i_{1}$ and $i_{2}$ are also defined as

$\displaystyle i_{1}(\phi(\alpha_{i}),f(v))=\int{\cal
D}\alpha_{i}\int_{0}^{1}dv\phi(\alpha_{q},\alpha_{\bar{q}},\alpha_{g})f(v)\theta(k-u_{0})~{},$
(44) $\displaystyle i_{2}(\phi(\alpha_{i}),f(v))=\int{\cal
D}\alpha_{i}\int_{0}^{1}dv\phi(\alpha_{q},\alpha_{\bar{q}},\alpha_{g})f(v)\delta(k-u_{0})~{},$
(45)

where $k=\alpha_{q}-v\alpha_{g}$.

Finally, matching both the physical and QCD sides of the CF leads to the sum
rules for the strong coupling constant as

$\displaystyle g_{B^{*}_{1}B^{*}_{2}{\cal
P}}=\frac{3m_{1}^{2}}{2}\dfrac{e^{m_{1}^{2}/M_{1}^{2}}e^{m_{2}^{2}/M_{2}^{2}}}{\lambda_{B^{*}_{1}}\lambda_{B^{*}_{2}}}\Pi_{B^{*}_{1}B^{*}_{2}{\cal
P}}(M^{2},s_{0}).$ (46)

##  3 Numerical Results

In this section, we present the numerical results of the sum rules for the
strong coupling constants of light pseudoscalar meson with doubly heavy
spin–$3/2$ baryons. For numerical evaluations, two sets of parameters are
being employed as follows: One which corresponds to light pseudoscalar mesons
that are their masses and decay constants as well as their non-perturbative
parameters appearing in their DAs of different twists, which are listed in
Tables (2) and (3) respectively. The other is the masses and residues of the
doubly heavy spin-3/2 baryons and are listed in Table (4).

Parameters | Values
---|---
$m_{s}$ | $93^{+11}_{-5}~{}\mbox{MeV}$
$m_{c}$ | $1.27^{+0.02}_{-0.02}~{}\mbox{GeV}$
$m_{b}$ | $4.18^{+0.03}_{-0.02}~{}\mbox{GeV}$
$m_{\pi^{0}}$ | $134.9768\pm 0.0005~{}\mbox{MeV}$
$m_{\pi^{\pm}}$ | $139.57039\pm 0.00018~{}\mbox{MeV}$
$m_{\eta}$ | $957.862\pm 0.017~{}\mbox{MeV}$
$m_{K^{0}}$ | $497.611\pm 0.013~{}\mbox{MeV}$
$m_{K^{\pm}}$ | $493.677\pm 0.016~{}\mbox{MeV}$
$f_{\pi}$ | $131~{}\mbox{MeV}$
$f_{\eta}$ | $130~{}\mbox{MeV}$
$f_{K}$ | $160~{}\mbox{MeV}$
Table 2: The meson masses and leptonic decay constants along with the quark masses [73]. meson | $a_{2}$ | $\eta_{3}$ | $w_{3}$ | $\eta_{4}$ | $w_{4}$
---|---|---|---|---|---
$\pi$ | 0.44 | $0.015$ | -3 | 10 | 0.2
K | 0.16 | $0.015$ | -3 | 0.6 | 0.2
$\eta$ | 0.2 | 0.013 | -3 | 0.5 | 0.2
Table 3: Input parameters for DAs of twist 2, 3 and 4 at the renormalization scale $\mu=1\text{GeV}$ [63, 64]. Baryon |  Mass $[\text{GeV}]$ | Residue $[\text{GeV}^{3}]$
---|---|---
$\Xi^{*}_{cc}$ | $3.69\pm 0.16$ | $0.12\pm 0.01$
$\Xi^{*}_{bc}$ | $7.25\pm 0.20$ | $0.15\pm 0.01$
$\Xi^{*}_{bb}$ | $10.14\pm 1.0$ | $0.22\pm 0.03$
$\Omega^{*}_{cc}$ | $3.78\pm 0.16$ | $0.14\pm 0.02$
$\Omega^{*}_{bc}$ | $7.3\pm 0.2$ | $0.18\pm 0.02$
$\Omega^{*}_{bb}$ | $10.5\pm 0.2$ | $0.25\pm 0.03$

Table 4: The masses and residues of spin-$3/2$ doubly heavy baryons [23].

The first step in the numerical analysis is to fix the working intervals of
the auxiliary parameters. These are continuum threshold $s_{0}$ and Borel
parameter $M^{2}$. To this end, we impose the conditions of: OPE convergence,
pole dominance and relatively weak dependence of the results on auxiliary
parameters. The threshold $s_{0}$, depends on the energy of the first excited
state at each channel. We have not any experimental information on the excited
states of the doubly heavy baryons, however, our analyses show that in the
interval $m_{B^{*}}+0.3~{}\text{GeV}\leq\sqrt{s_{0}}\leq
m_{B^{*}}+0.5~{}\text{GeV}$ the standard requirements are satisfied. To fix
the Borel parameter $M^{2}$ we follow the steps: for its upper limit we impose
the condition of pole dominance and for determination of its lower limit we
demand the OPE convergence condition. These requirements lead to the working
regions $14$ GeV2 $\leq M^{2}$ $\leq 18$ GeV 2, $8$ GeV2 $\leq M^{2}$ $\leq
11$ GeV 2 and $4$ GeV2 $\leq M^{2}$ $\leq 6$ GeV 2 for $bb$, $bc$ and $cc$
channels, respectively.

In Fig. (1), as an example, we show the dependence of the coupling constant
$g$ on the continuum threshold at some fixed values of the Borel parameter for
the vertices $\Xi^{*}_{bb}\Xi^{*}_{bb}\pi^{\pm}$,
$\Omega^{*}_{bb}\Omega^{*}_{bb}\eta$, $\Omega^{*}_{bb}\Xi^{*}_{bb}K^{\pm}$,
and $\Omega^{*}_{bb}\Xi^{*}_{bb}K^{0}$. As is seen, the strong coupling
constants depend very weakly on $s_{0}$.

Figure 1: Dependence of the strong coupling constant $g$ on the continuum
threshold at different values of the Borel parameter $M^{2}=14,16$ and $18$
GeV2.

As another example, we show the dependence of $g$ on the Borel parameter
$M^{2}$ for the vertices $\Xi^{*}_{bb}\Xi^{*}_{bb}\pi^{\pm}$,
$\Omega^{*}_{bb}\Omega^{*}_{bb}\eta$, $\Omega^{*}_{bb}\Xi^{*}_{bb}K^{\pm}$,
and $\Omega^{*}_{bb}\Xi^{*}_{bb}K^{0}$ at fixed values of $s_{0}$ in Fig (2).
From these plots, we see that the strong coupling constants show some
dependence on the Borel parameter, which constitute the main parts of the
uncertainties. These uncertainties, however, remain inside the limits allowed
by the method.

Figure 2: Dependence of the strong coupling constant $g$ on the Borel
parameter at different values of the continuum threshhold $s_{0}$. The lines
correspond to $s_{0}=109$, $113$ and $117$ GeV2 for
$\Xi^{*}_{bb}\Xi^{*}_{bb}\pi^{\pm}$ channel and $s_{0}=116$, $121$ and $125$
GeV2 for $\Omega^{*}_{bb}\Omega^{*}_{bb}\eta$,
$\Omega^{*}_{bb}\Xi^{*}_{bb}K^{\pm}$, and
$\Omega^{*}_{bb}\Xi^{*}_{bb}K^{0}$ones.

Taking into account all the input values as well as the working intervals for
the auxiliary parameters, we extract the numerical values for the strong
coupling constants from the light-cone sum rules (46). The numerical results
of $g$ for all the considered vertices in $cc$, $bc$ and $bb$ channels are
presented in Table (5).

Vertex |  $M^{2}$(GeV2) |  $s_{0}$ (GeV2) | strong coupling constant
---|---|---|---
Decays to $\pi$
$\Xi^{*}_{bb}\Xi^{*}_{bb}\pi^{0}$ | $14\leq M^{2}\leq 18$ | $109\leq s_{0}\leq 117$ | $50.0^{\>7.9}_{\>7.8}$
$\Xi^{*}_{bb}\Xi^{*}_{bb}\pi^{\pm}$ | $14\leq M^{2}\leq 18$ | $109\leq s_{0}\leq 117$ | $70.8^{\>11.4}_{\>10.3}$
$\Xi^{*}_{bc}\Xi^{*}_{bc}\pi^{0}$ | $7\leq M^{2}\leq 10$ | $57\leq s_{0}\leq 63$ | $13.4^{\>1.4}_{\>1.4}$
$\Xi^{*}_{bc}\Xi^{*}_{bc}\pi^{\pm}$ | $7\leq M^{2}\leq 10$ | $57\leq s_{0}\leq 63$ | $18.9^{\>2.0}_{\>2.0}$
$\Xi^{*}_{cc}\Xi^{*}_{cc}\pi^{0}$ | $3\leq M^{2}\leq 6$ | $16\leq s_{0}\leq 19$ | $3.8^{\>0.3}_{\>0.4}$
$\Xi^{*}_{cc}\Xi^{*}_{cc}\pi^{\pm}$ | $3\leq M^{2}\leq 6$ | $16\leq s_{0}\leq 19$ | $5.5^{\>0.5}_{\>0.6}$
Decays to $\eta$
$\Omega^{*}_{bb}\Omega^{*}_{bb}\eta$ | $14\leq M^{2}\leq 18$ | $116\leq s_{0}\leq 125$ | $125.7^{\>22.3}_{\>23.7}$
$\Omega^{*}_{bc}\Omega^{*}_{bc}\eta$ | $7\leq M^{2}\leq 10$ | $57\leq s_{0}\leq 64$ | $18.4^{\>2.1}_{\>2.0}$
$\Omega^{*}_{cc}\Omega^{*}_{cc}\eta$ | $3\leq M^{2}\leq 6$ | $16\leq s_{0}\leq 20$ | $5.9^{\>0.5}_{\>0.6}$
Decays to $K$
$\Omega^{*}_{bb}\Xi^{*}_{bb}K^{0}$ | $14\leq M^{2}\leq 18$ | $116\leq s_{0}\leq 125$ | $142.9^{\>26.6}_{\>23.4}$
$\Omega^{*}_{bb}\Xi^{*}_{bb}K^{\pm}$ | $14\leq M^{2}\leq 18$ | $116\leq s_{0}\leq 125$ | $142.3^{\>26.4}_{\>23.2}$
$\Omega^{*}_{bc}\Xi^{*}_{bc}K^{0}$ | $7\leq M^{2}\leq 10$ | $57\leq s_{0}\leq 64$ | $27.5^{\>3.0}_{\>2.9}$
$\Omega^{*}_{bc}\Xi^{*}_{bc}K^{\pm}$ | $7\leq M^{2}\leq 10$ | $57\leq s_{0}\leq 64$ | $27.3^{\>3.0}_{\>2.9}$
$\Omega^{*}_{cc}\Xi^{*}_{cc}K^{0}$ | $3\leq M^{2}\leq 6$ | $16\leq s_{0}\leq 20$ | $8.7^{\>0.8}_{\>0.9}$
$\Omega^{*}_{cc}\Xi^{*}_{cc}K^{\pm}$ | $3\leq M^{2}\leq 6$ | $16\leq s_{0}\leq 20$ | $8.6^{\>0.7}_{\>0.9}$
Table 5: Working regions of the Borel parameter $M^{2}$ and continuum
threshold $s_{0}$ as well as the numerical values for different strong
coupling constants extracted from the analyses.

There are three sources of errors in the results. The first is from the input
parameters we have used in the calculations. The second, the intrinsic error,
is due to the employment of the quark-hadron duality assumption and finally,
the third one, is the error because of the variations with respect to the
auxiliary parameters in their working intervals. The results contain errors
overall up to $20\%$.

##  4 Summary and Conclusions

In the present study, we investigated the strong coupling constants of the
doubly heavy spin-$3/2$ baryons, $\Omega^{*}_{QQ^{\prime}}$ and
$\Xi^{*}_{QQ^{\prime}}$, with the light pseudoscalar mesons, $\pi$, $K$ and
$\eta$, by applying the LCSR formalism. We used the interpolating currents of
the doubly heavy baryons and DAs of the mesons to construct the sum rules for
the strong coupling constants. By fixing the auxiliary parameters, entered the
calculations through the Borel transformation as well as continuum
subtraction, we extracted the numerical results for different possible
vertices considering the quark contents of the participating particles and
other considerations.

The obtained results indicate that the heavier the doubly heavy baryons, the
larger their strong couplings to the same pseudoscalar mesons. The significant
difference between the c- and b-baryons can be explained using the heavy-quark
spin symmetry that tells us the partial strong decay widths of heavy quark
hadronic systems are independent of the heavy quark flavor and spin [74].
Considering the Fermi golden rule and the fact that the amplitude is
proportional to the strong coupling constant, defining the corresponding
vertex, immediately implies that the heavier doubly-heavy baryon has the
larger strong coupling constant. Also comparing with the results presented in
Ref. [47], one finds out that the values of the strong coupling constants
among the doubly heavy spin-3/2 baryons and light pseudoscalar mesons are
mostly larger than that of the spin-1/2 baryons with the same mesons taken
into account in Ref. [47].

Our results may help experimental groups in analyses of the data including the
hadrons taken into account in the present study. The results may also be used
in the construction of the strong potentials among the doubly heavy baryons
and pseudoscalar mesons.

## ACKNOWLEDGEMENTS

K. Azizi and S. Rostami are thankful to Iran Science Elites Federation
(Saramadan) for the financial support provided under the grant number
ISEF/M/99171.

## References

  * [1] W. Wang, F. S. Yu and Z. X. Zhao, “Weak decays of doubly heavy baryons: the $1/2\rightarrow 1/2$ case,”  Eur. Phys. J. C 77, no. 11, 781 (2017).
  * [2] W. Wang, Z. P. Xing and J. Xu, “Weak Decays of Doubly Heavy Baryons: SU(3) Analysis,”  Eur. Phys. J. C 77, no. 11, 800 (2017).
  * [3] T. Gutsche, M. A. Ivanov, J. G. Körner and V. E. Lyubovitskij, “Decay chain information on the newly discovered double charm baryon state $\Xi_{cc}^{++}$,”  Phys. Rev. D 96, no. 5, 054013 (2017).
  * [4] H. S. Li, L. Meng, Z. W. Liu and S. L. Zhu, “Radiative decays of the doubly charmed baryons in chiral perturbation theory,”  Phys. Lett. B 777, 169 (2018).
  * [5] L. Y. Xiao, K. L. Wang, Q. f. Lu, X. H. Zhong and S. L. Zhu, “Strong and radiative decays of the doubly charmed baryons,”  Phys. Rev. D 96, no. 9, 094005 (2017).
  * [6] N. Sharma and R. Dhir, “Estimates of W-exchange contributions to $\Xi_{cc}$ decays,”  Phys. Rev. D 96, no. 11, 113006 (2017).
  * [7] Y. L. Ma and M. Harada, “Chiral partner structure of doubly heavy baryons with heavy quark spin-flavor symmetry,”  J. Phys. G 45, no. 7, 075006 (2018).
  * [8] X. H. Hu, Y. L. Shen, W. Wang and Z. X. Zhao, “Weak decays of doubly heavy baryons: ”decay constants”,”  Chin. Phys. C 42, no. 12, 123102 (2018).
  * [9] Y. J. Shi, W. Wang, Y. Xing and J. Xu, “Weak Decays of Doubly Heavy Baryons: Multi-body Decay Channels,”  Eur. Phys. J. C 78, no. 1, 56 (2018).
  * [10] D. L. Yao, “Masses and sigma terms of doubly charmed baryons up to $\mathcal{O}(p^{4})$ in manifestly Lorentz-invariant baryon chiral perturbation theory,”  Phys. Rev. D 97, no. 3, 034012 (2018).
  * [11] Z. X. Zhao, “Weak decays of doubly heavy baryons: the $1/2\rightarrow 3/2$ case,”  Eur. Phys. J. C 78, no. 9, 756 (2018).
  * [12] M. Z. Liu, Y. Xiao and L. S. Geng, “Magnetic moments of the spin-1/2 doubly charmed baryons in covariant baryon chiral perturbation theory,”  Phys. Rev. D 98, no. 1, 014040 (2018).
  * [13] Z. P. Xing and Z. X. Zhao, “Weak decays of doubly heavy baryons: the FCNC processes,” Phys. Rev. D 98, no. 5, 056002 (2018).
  * [14] R. Dhir and N. Sharma, “Weak decays of doubly heavy charm $\Omega_{cc}^{+}$ baryon,” Eur. Phys. J. C 78, no. 9, 743 (2018).
  * [15] A. V. Berezhnoy, A. K. Likhoded and A. V. Luchinsky, “Doubly heavy baryons at the LHC,” Phys. Rev. D 98, no. 11, 113004 (2018).
  * [16] L. J. Jiang, B. He and R. H. Li, “Weak decays of doubly heavy baryons: $\mathcal{B}_{cc}\rightarrow\mathcal{B}_{c}V$,” Eur. Phys. J. C 78, no. 11, 961 (2018).
  * [17] Q. A. Zhang, “Weak Decays of Doubly Heavy Baryons: W-Exchange,” Eur. Phys. J. C 78, no. 12, 1024 (2018).
  * [18] T. Gutsche, M. A. Ivanov, J. G. Körner, V. E. Lyubovitskij and Z. Tyulemissov, “Ab initio three-loop calculation of the $W$-exchange contribution to nonleptonic decays of double charm baryons,” Phys. Rev. D 99, no. 5, 056013 (2019).
  * [19] Y. J. Shi, Y. Xing and Z. X. Zhao, “Light-cone sum rules analysis of $\Xi_{QQ^{\prime}q}\to\Lambda_{Q^{\prime}}$ weak decays,” Eur. Phys. J. C 79, no. 6, 501 (2019).
  * [20] X. H. Hu and Y. J. Shi, “Light-Cone Sum Rules Analysis of $\Xi_{QQ^{\prime}}\to\Sigma_{Q^{\prime}}$ Weak Decays,” Eur. Phys. J. C 80, no. 1, 56 (2020).
  * [21] S. J. Brodsky, F. K. Guo, C. Hanhart and U.-G. Meißner, “Isospin splittings of doubly heavy baryons,” Phys. Lett. B 698, 251 (2011).
  * [22] T. M. Aliev, K. Azizi and M. Savci, “Doubly Heavy Spin–1/2 Baryon Spectrum in QCD,”  Nucl. Phys. A 895, 59 (2012).
  * [23] T. M. Aliev, K. Azizi and M. Savci, “The masses and residues of doubly heavy spin-3/2 baryons,”  J. Phys. G 40, 065003 (2013).
  * [24] K. Azizi, T. M. Aliev and M. Savci, “Properties of doubly and triply heavy baryons,”  J. Phys. Conf. Ser. 556, no. 1, 012016 (2014).
  * [25] T. M. Aliev, K. Azizi and M. Savcı, “Mixing angle of doubly heavy baryons in QCD,” Phys. Lett. B 715, 149-151 (2012.
  * [26] U. Özdem, “Magnetic dipole moments of the spin-$\frac{3}{2}$ doubly heavy baryons,” Eur. Phys. J. A 56, no.2, 34 (2020).
  * [27] M. Mattson et al. [SELEX Collaboration], “First Observation of the Doubly Charmed Baryon $\Xi^{+}_{cc}$,”  Phys. Rev. Lett. 89, 112001 (2002).
  * [28] A. Ocherashvili et al. [SELEX Collaboration], “Confirmation of the double charm baryon $\Xi^{+}_{cc}$(3520) via its decay to $pD^{+}K^{-}$,”  Phys. Lett. B 628, 18 (2005).
  * [29] S. P. Ratti, “New results on c-baryons and a search for cc-baryons in FOCUS,”  Nucl. Phys. Proc. Suppl. 115, 33 (2003).
  * [30] B. Aubert et al. [BaBar Collaboration], “Search for doubly charmed baryons $\Xi_{cc}^{+}$ and $\Xi_{cc}^{++}$ in BABAR,”  Phys. Rev. D 74, 011103 (2006).
  * [31] R. Aaij et al. [LHCb Collaboration], “Search for the doubly charmed baryon $\Xi_{cc}^{+}$,”  JHEP 1312, 090 (2013).
  * [32] Y. Kato et al. [Belle Collaboration], “Search for doubly charmed baryons and study of charmed strange baryons at Belle,”  Phys. Rev. D 89, no. 5, 052003 (2014).
  * [33] A. I. Onishchenko, “Doubly heavy systems: Decays and OPE,”  arXiv:hep-ph/9912424 [hep-ph].
  * [34] C. H. Chang, T. Li, X. Q. Li and Y. M. Wang, “Lifetime of doubly charmed baryons,”  Commun. Theor. Phys. 49, 993-1000 (2008).
  * [35] R. Aaij et al. [LHCb Collaboration], “Observation of the doubly charmed baryon $\Xi_{cc}^{++}$,”  Phys. Rev. Lett. 119, no. 11, 112001 (2017).
  * [36] R. Aaij et al. [LHCb Collaboration], “First Observation of the Doubly Charmed Baryon Decay $\Xi_{cc}^{++}\rightarrow\Xi_{c}^{+}\pi^{+}$,”  Phys. Rev. Lett. 121, no. 16, 162002 (2018).
  * [37] R. Aaij et al. [LHCb], “Search for the doubly heavy $\mathit{\Xi}_{bc}^{0}$ baryon via decays to $D^{0}pK^{-}$,”  arXiv:2009.02481 [hep-ex].
  * [38] M. Gell-Mann, A schematic model of baryons and mesons Physics Letters 8 (3) (1964) 214 – 215.
  * [39] G. Zweig, An SU(3) model for strong interaction symmetry and its breaking. Version 1 https://inspirehep.net/record/11881/files/CM-P00042883.pdf.
  * [40] D. Lichtenberg, S. P. Rosen, G. Zweig, DEVELOPMENTS IN THE QUARK THEORY OF HADRONS https://inspirehep.net/record/4674/files/cern-th-412.pdf.
  * [41] Z. S. Brown, W. Detmold, S. Meinel and K. Orginos, “Charmed bottom baryon spectroscopy from lattice QCD,” Phys. Rev. D 90, no.9, 094507 (2014).
  * [42] E. Hernandez, J. Nieves and J. M. Verde-Velasco, “Heavy quark symmetry constraints on semileptonic form-factors and decay widths of doubly heavy baryons,” Phys. Lett. B 663 (2008), 234-241.
  * [43] J. M. Flynn, E. Hernandez and J. Nieves, “Triply Heavy Baryons and Heavy Quark Spin Symmetry,” Phys. Rev. D 85 (2012), 014012.
  * [44] K. Azizi, A. R. Olamaei and S. Rostami, “Beautiful mathematics for beauty-full and other multi-heavy hadronic systems,” Eur. Phys. J. A 54, no. 9, 162 (2018).
  * [45] A. R. Olamaei, K. Azizi and S. Rostami, “Strong coupling constants of the doubly heavy $\Xi_{QQ}$ Baryons with $\pi$ Meson,” Eur. Phys. J. C 80, no.7, 613 (2020).
  * [46] K. Azizi, A. R. Olamaei and S. Rostami, “Strong interaction of doubly heavy spin-3/2 baryons with light vector mesons,” Eur. Phys. J. C 80, no.12, 1196 (2020).
  * [47] S. Rostami, K. Azizi and A. R. Olamaei, “Strong Coupling Constants of the Doubly Heavy Spin-1/2 Baryons with Light Pseudoscalar Mesons,” Chin. Phys. C 45, no.2, 023120 (2021).
  * [48] K. Azizi and H. Sundu, “Mass and magnetic dipole moment of negative parity heavy baryons with spin–3/2,” Eur. Phys. J. Plus 132, no.1, 22 (2017).
  * [49] K. Azizi, N. Er and H. Sundu, “Positive and negative parity hyperons in nuclear medium,” Phys. Rev. D 92, no.5, 054026 (2015).
  * [50] K. Azizi, Y. Sarac and H. Sundu, “Strong couplings of negative and positive parity nucleons to the heavy baryons and mesons,” Phys. Rev. D 92, no.1, 014022 (2015).
  * [51] T. M. Aliev, K. Azizi, T. Barakat and M. Savcı, “Diagonal and transition magnetic moments of negative parity heavy baryons in QCD sum rules,” Phys. Rev. D 92, no.3, 036004 (2015).
  * [52] T. M. Aliev, K. Azizi and M. Savci, “Magnetic moments of negative parity heavy baryons in QCD,” arXiv:1504.00172 [hep-ph].
  * [53] K. Azizi and H. Sundu, “Radiative transition of negative to positive parity nucleon,” Phys. Rev. D 91, no.9, 093012 (2015).
  * [54] A. Khodjamirian, C. Klein, T. Mannel and Y. M. Wang, “Form Factors and Strong Couplings of Heavy Baryons from QCD Light-Cone Sum Rules,” JHEP 09, 106 (2011).
  * [55] T. M. Aliev and K. Şimşek, “Strong vertices of doubly heavy spin-3/2-spin-1/2 baryons with light mesons in light-cone QCD sum rules,” [arXiv:2011.07150 [hep-ph].
  * [56] T. M. Aliev and K. Şimşek, “Strong coupling constants of doubly heavy baryons with vector mesons in QCD,” Eur. Phys. J. C 80, no.10, 976 (2020).
  * [57] H. I. Alrebdi, T. M. Aliev and K. Şimşek, “Determination of the strong vertices of doubly heavy baryons with pseudoscalar mesons in QCD,” Phys. Rev. D 102, no.7, 074007 (2020).
  * [58] T. M. Aliev, T. Barakat and K. Şimşek, “Strong $B_{QQ^{\prime}}^{*}B_{QQ^{\prime}}V$ vertices and the radiative decays of $B_{QQ}^{*}\to B_{QQ}\gamma$ in the light-cone sum rules,” arXiv:2101.10264 [hep-ph].
  * [59] M. A. Shifman, A. I. Vainshtein and V. I. Zakharov, “QCD and Resonance Physics. Theoretical Foundations,” Nucl. Phys. B 147, 385-447 (1979).
  * [60] I. Balitsky, V. M. Braun and A. Kolesnichenko, “Radiative Decay $\Sigma^{+}\rightarrow p\gamma$ in Quantum Chromodynamics,” Nucl. Phys. B 312, 509-550 (1989).
  * [61] A. Khodjamirian, R. Ruckl, S. Weinzierl and O. I. Yakovlev, “Perturbative QCD correction to the $B\rightarrow\pi$ transition form-factor,”  Phys. Lett. B 410, 275 (1997)
  * [62] V. M. Braun, “Light cone sum rules,” arXiv:hep-ph/9801222 [hep-ph].
  * [63] P. Ball and R. Zwicky, “SU(3) breaking of leading-twist K and K* distribution amplitudes: A Reprise,” Phys. Lett. B 633, 289-297 (2006).
  * [64] P. Ball and R. Zwicky, “New results on $B\to\pi,K,\eta$ decay formfactors from light-cone sum rules,”  Phys. Rev. D 71, 014015 (2005).
  * [65] P. Ball, “Theoretical update of pseudoscalar meson distribution amplitudes of higher twist: The Nonsinglet case,” JHEP 9901, 010 (1999).
  * [66] H. D. Li, C. D. Lü, C. Wang, Y. M. Wang and Y. B. Wei, “QCD calculations of radiative heavy meson decays with subleading power corrections,” JHEP 04, 023 (2020).
  * [67] A. Khodjamirian, B. Melić, Y. M. Wang and Y. B. Wei, “The $D^{*}D\pi$ and $B^{*}B\pi$ couplings from light-cone sum rules,” JHEP 03, 016 (2021).
  * [68] T. M. Aliev, K. Azizi and M. Savci, “Strong coupling constants of heavy spin–3/2 baryons with light pseudoscalar mesons,” Nucl. Phys. A 870-871, 58-71 (2011)..
  * [69] V. M. Belyaev, V. M. Braun, A. Khodjamirian and R. Ruckl, “D* D $\pi$ and B* B $\pi$ couplings in QCD,” Phys. Rev. D 51, 6177 (1995).
  * [70] T. M. Aliev, K. Azizi and M. Savci, “Strong coupling constants of light pseudoscalar mesons with heavy baryons in QCD,” Phys. Lett. B 696, 220 (2011).
  * [71] P. Ball, V. M. Braun, “Next-to-leading-order corrections to meson masses in the heavy quark effective theory,” Phys. Rev. D 49, 2472 (1994).
  * [72] S. S. Agaev, K. Azizi and H. Sundu, “Application of the QCD light cone sum rule to tetraquarks: the strong vertices $X_{b}X_{b}\rho$ and $X_{c}X_{c}\rho$,” Phys. Rev. D 93, no.11, 114036 (2016).
  * [73] P.A. Zyla et al. [Particle Data Group], “Review of Particle Physics,” Prog. Theor. Exp. Phys. 2020, 083C01 (2020).
  * [74] N. Isgur and M. B. Wise, “Spectroscopy with heavy-quark symmetry,” Phys. Rev. Lett. 66, 1130 (1991).


Please address correspondence to Thomas J. Faulkenberry, Department of
Psychological Sciences, Box T-0820, Tarleton State University, Stephenville,
TX 76402. Email<EMAIL_ADDRESS>All raw data and R scripts are
available to download at http://github.com/tomfaulkenberry/bayesRegression/.
The author would like to thank Lisa Fazio for providing data for Model 2.

# A hierarchical Bayesian model for measuring individual-level and group-level
numerical representations

Thomas J. Faulkenberry  Department of Psychological Sciences, Tarleton State
University, USA

###### Abstract

A popular method for indexing numerical representations is to compute an
individual estimate of a response time effect, such as the SNARC effect or the
numerical distance effect. Classically, this is done by estimating individual
linear regression slopes and then either pooling the slopes to obtain a group-
level slope estimate, or using the individual slopes as predictors of other
phenomena. In this paper, I develop a hierarchical Bayesian model for
simultaneously estimating group-level and individual-level slope parameters. I
show examples of using this modeling framework to assess two common effects in
numerical cognition: the SNARC effect and the numerical distance effect.
Finally, I demonstrate that the Bayesian approach can result in better
measurement fidelity than the classical approach, especially with small
samples.

Keywords: hierarchical regression, Bayesian modeling, numerical
representations, SNARC effect, distance effect

The purpose of this paper is to introduce a hierarchical Bayesian model for
measuring two types of mental representations that are often assessed in
numerical cognition: the SNARC effect (Dehaene ., 1990) and the numerical
distance effect (Moyer Landauer, 1967). The SNARC effect (Spatial-Numerical
Associations of Response Codes) refers to the finding that response times
(RTs) for small numerical stimuli are faster with the left hand, whereas RTs
for large numbers are faster with the right hand. The numerical distance
effect (NDE) refers to the finding that RTs decrease as the numerical distance
between stimulus numbers increases. Both effects are thought to reflect
particular aspects of mental representation of number. As the SNARC effect
occurs even when number magnitude is not salient to a given task (e.g., asking
whether a given number is odd or even), the SNARC effect is often taken as an
index of automatic processing of number magnitude (Wood ., 2008). Similarly,
the NDE is thought to represent a noisy, analog representation of number
(Dehaene, 1992), and thus is purported to occur because as numerical distance
increases, the amount of representational overlap (and hence response
competition) decreases, resulting in faster RTs (Gevers ., 2006).

As indices of numerical representation, both the SNARC effect and the NDE have
been used to study individual differences in number representation (Holloway
Ansari, 2009; Viarouge ., 2014). For example, Viarouge . (2014) showed that
participants with stronger SNARC effects tended to exhibit slower mental
rotation processing speed. Similarly, Pinhas . (2014) found that the SNARC
effect was positively associated with spatial operational momentum. The NDE
has been used in a similar manner. Holloway Ansari (2009) found that children
with a larger NDE tended to score lower on standardized assessments of
mathematical fluency and calculation skill. Fazio . (2014) showed a similar
pattern, namely that the NDE was negatively correlated with mathematics
achivement, but not reading achievement. As such, it is clear that accurate
measurement of the SNARC effect and NDE is a desirable objective for a wide
range of research paradigms within mathematical cognition.

The classical method for measuring the SNARC effect and NDE for individuals is
based on a procedure known as regression coefficient analysis (RCA; Lorch
Myers, 1990). In this procedure, a dependent measure such as RT is regressed
against a predictor for each individual participant. The estimated slope and
intercept for each individual is then recorded, at which point any number of
analyses may be performed. One of the first examples of this procedure in
numerical cognition comes from Fias (1996), who adapted this method to measure
the SNARC effect. Participants’ left-hand and right-hand RTs for each number
stimulus (e.g., 1, 2, 8, 9) were aggregated via the median. Then, the
difference between right- and left-hand RTs ($dRT$) was computed for each
number. The SNARC effect was then presented as a negative correlation between
number and dRT. That is, for small numbers (e.g., 1,2), the left hand is
faster, so $dRT>0$. For large numbers, the right hand is faster, so $dRT<0$.
This negative correlation is captured by a linear regression, which is used to
calculate the slope on $dRT$ for each participant. Analyses of the NDE are
similar, but instead of using $dRT$, one typically regresses $RT$ against a
predictor such as numerical distance (e.g., Holloway Ansari, 2009) or ratio
(e.g., Fazio ., 2014). As with the SNARC effect, the NDE typically presents as
a negative slope. These slopes can then be submitted to a variety of analyses,
including group-level tests (e.g., does the SNARC/NDE effect vary by group?)
or individual-level correlations (i.e., is the SNARC/NDE effect correlated
with mathematical performance?).

In the following sections, I describe a hierarchical Bayesian approach to
measuring the SNARC effect and the NDE. Certainly, a complete treatment of
Bayesian methods is beyond the scope of this paper, but the interested reader
is advised consult the excellent books by Gelman . (2013), McElreath (2015),
or Kruschke (2015) for more details. However, I will briefy outline some
advantages to using such an approach to measuring numerical representations.
One such advantage is that in a hierarchical model, variability from
participants and items is modeled simultaneously, resulting in in both
participant-level and group-level parameter estimates. As such, the
hierarchical model allows the researcher to avoid aggregating single estimates
across individuals, which can be problematic (Heathcote ., 2000; Haider
Frensch, 2002; Pratte ., 2010). Another advantage is both philosophical and
pragmatic. Bayesian modeling provides a principled means to combine prior
knowledge with new data via Bayes’ theorem. The result is a quantification of
_posterior_ belief in terms of a probability distribution that reflects
uncertainty after seeing data. From this distribution, one can estimate a
variety of summary statistics, including the posterior mode (i.e., the
parameter estimate with highest density), and Bayesian credible intervals,
which give a range of values containing a particular probability
mass.111Bayesian credible intervals are not the same thing as frequentist
_confidence_ intervals, which are defined as a range of parameters that would
contain a population parameter in some specified percentage of a large number
of repeated samples. While people often misinterpret confidence intervals
using probability statements (Hoekstra ., 2014), Bayesian credible intervals
_can_ be directly interpreted as probability statements One particular type of
credible interval used in a Bayesian context is the highest posterior density
interval (HPDI), which is the _narrowest_ interval containing a given
probability mass.

Generally, the models will be constructed as follows. As with the classical
method described above, the goal is to estimate a slope on $dRT$ (for the
SNARC effect) or $RT$ (for the NDE). However, the proposed models are
hierarchical, which means that individual-level slopes are drawn from a group-
level distribution that specifies how the slopes are distributed in the
population. Thus, individual-level and group-level estimates are modeled
simuntaneously, which improves upon the noisy slope estimates usually obtained
from classical regression with small samples (Gelman Hill, 2013). Further,
since the model is Bayesian, one can specify group-level prior distributions
that will moderate extreme individual-level parameter estimates (i.e.,
shrinkage), which will result in increased measurement accuracy for our
parameter estimates.

##  1 Model 1: estimating the SNARC effect

The model is a hierarchical Bayesian linear regression model, where $dRT$ is
predicted by stimulus number. A graphical representation of the model can be
seen in Figure 1. Formally, we define

$dRT_{ij}=\alpha_{i}+\beta_{i}(j)+\varepsilon_{ij}$ (1)

where $i=$ subject number and $j=$ stimulus number. The residuals
$\varepsilon_{ij}$ are assumed to be normally distributed with mean 0 and
precision $1/\sigma^{2}$. Thus, we can express the likelihood for the data as

$dRT_{ij}\sim\text{Normal}(\mu_{ij},1/\sigma^{2})$ (2)

where $\mu_{ij}=\alpha_{i}+\beta_{i}(j)$. The prior for each individual-level
intercept $\alpha_{i}$ is set to be uniformly distributed between -200 and
200. The hierarchical structure is instantiated on slope. First, I set the
prior for each individual-level slope $\beta_{i}$ to be normally distributed
with two hyperparameters: mean $b$ and precision $1/\sigma_{b}^{2}$. In turn,
this requires priors on the hyperparameters – the prior for $b$ is uniformly
distributed between -20 and 20. Note that this range is based on inspection of
the slopes obtained in Viarouge . (2014). The prior for group-level slop
precision $1/\sigma_{b}^{2}$ (as well as group-level residual precision
$1/\sigma^{2}$) is set as a Gamma distribution with both shape and scale equal
to 0.01.

$dRT_{ij}$$\mu_{ij}$$\beta_{i}$$\alpha_{i}$$b$$\sigma_{b}^{2}$$\sigma^{2}$number
$j$subject $i$ $\displaystyle b$ $\displaystyle\sim\text{Uniform}(-20,20)$
$\displaystyle 1/\sigma_{b}^{2}$ $\displaystyle\sim\text{Gamma}(0.01,0.01)$
$\displaystyle\alpha_{i}$ $\displaystyle\sim\text{Uniform}(-200,200)$
$\displaystyle\beta_{i}$ $\displaystyle\sim\text{Gaussian}(b,\sigma_{b}^{2})$
$\displaystyle\mu_{ij}$ $\displaystyle=\alpha_{i}+\beta_{i}(j)$ $\displaystyle
1/\sigma^{2}$ $\displaystyle\sim\text{Gamma}(0.01,0.01)$ $\displaystyle
dRT_{ij}$ $\displaystyle\sim\text{Gaussian}(\mu_{ij},\sigma^{2})$ Figure 1:
Graphical model of the hierarchical Bayesian linear regression model for the
SNARC effect. Following the convention of Lee Wagenmakers (2014), nodes
represent variables of interest (observed=shaded, latent=unshaded), with
dependencies represented via the graph structure.

###  1.1 Fitting the model

####  1.1.1 Data

I fit the model to data collected from 35 participants in a number parity
task. The numbers 1, 2, 8, and 9 were presented in the center of a computer
screen, after which participants were asked to quickly indicate via a button
press whether the number was even or odd. The procedure mirrored that of
Experiment 1 of Pinhas . (2014). Participants completed 112 trials of the task
under two counterbalanced response rules (either even=left or even=right).
This resulted in a collection of 3,920 trials. We removed 259 error trials and
an additional 12 trials for which RT exceeded 3 seconds (a total of 6.9% of
trials). The remaining 3,649 trials were collapsed into $35\times 8\times 2$
cells by computing median $RT$ for each of the conditions defined by crossing
the factors of subject, stimulus number, and response hand. Then, $dRT$ was
computed for each combination of subject and number by subtracting left-hand
RT from right-hand RT.

####  1.1.2 Results

The regression model parameters were estimated using R (R Core Team, 2016) and
JAGS (Plummer, 2017). Posterior sampling consisted of 3 MCMC chains, each
containing 100,000 draws. The first 5000 draws of each chain were discarded as
“burn-in” samples, leaving 285,000 samples remaining. These remaining samples
were thinned by a factor of 10, leaving a final sample of 28,500 posterior
draws for each parameter. Visual inspection of trace plots indicated that all
chains converged appropriately. Additionally, the Gelman-Rubin statistic
$\hat{R}=1.001$ for each parameter, indicating that the Markov chains for each
parameter converged to the appropriate stationary distribution (Gelman Rubin,
1992; Gelman ., 2013).

Since the model is hierarchical, I was able to estimate posterior
distributions for each $\alpha_{i}$ and $\beta_{i}$ (i.e., each participant’s
intercept and slope, respectively). Further, I estimated the posterior
distribution of $b$, which is the group level mean slope. The flexibility of
this model allows one to ask many questions about the SNARC effect, both at
the group level and the individual level. To illustrate, I will investigate
whether there was an overall SNARC effect for the group (c.f., Fias, 1996).
This can be answered by looking at the posterior distribution of the group-
level $b$.

The posterior distribution of the group-level slope $b$ is depicted in Figure
2. As can be seen in the figure, the mass of the distribution is centered over
the posterior mode $b=-11.5$. Further, the 95% HPDI for the slope $b$ is
$[-15.6,-7.3]$. Finally, since the posterior distribution is a _probability
distribution_ , we can compute the probability that $b<0$ (that is, the
probability that there is a non-zero SNARC effect). This probability is
greater than 0.999. Other Bayesian tools can be applied to this model, such as
computing a Bayes factor comparing a null-SNARC model to our obtained model
via the Savage-Dickey density ratio (Wagenmakers ., 2010). Briefly, this
method amounts to computing the density of the value $b=0$ in the posterior
distribution divided by the density of $b=0$ in the prior (i.e.,
$\text{Uniform}(-20,20)$). Doing this resulted in a Bayes factor of
approximately 280,000 to 1 in favor of the alternative hypothesis, indicating
overwhelming support for an overall SNARC effect.

Figure 2: Posterior distribution of the group-level slope $b$ for Model 1,
indicating a substantial SNARC effect across the group.

##  2 Model 2: estimating the numerical distance effect

As with Model 1, this model is also a hierarchical Bayesian linear regression
model, where $RT$ is predicted by numerical distance. For this specific
instantiation of the model, I will index numerical distance via the _ratio_
between compared numbers. For example, Fazio . (2014) used four ratio bins as
predictors of RT. A graphical representation of the model can be seen in
Figure 3. Formally, we define

$RT_{ij}=\alpha_{i}+\beta_{i}(j)+\varepsilon_{ij}$ (3)

where $i=$ participant number and $j=1,\dots,4$ equals ratio bin number. The
residuals $\varepsilon_{ij}$ are again assumed to be normally distributed with
mean 0 and precision $1/\sigma^{2}$. The likelihood for the data is

$RT_{ij}\sim\text{Normal}(\mu_{ij},1/\sigma^{2})$ (4)

where $\mu_{ij}=\alpha_{i}+\beta_{i}(j)$. The prior for each individual-level
intercept $\alpha_{i}$ is set to be uniformly distributed between 0 and 2000.
The prior for each individual-level slope $\beta_{i}$ to be normally
distributed with two hyperparameters: mean $b$ (with uniform prior between
-100 and 100) and precision $1/\sigma_{b}^{2}$ (with a $\text{Gamma}(.01,.01)$
prior). Finally, the prior for the group-level residual precision is
$\text{Gamma}(.01,.01)$.

$RT_{ij}$$\mu_{ij}$$\beta_{i}$$\alpha_{i}$$b$$\sigma_{b}^{2}$$\sigma^{2}$number
$j$subject $i$ $\displaystyle b$ $\displaystyle\sim\text{Uniform}(-100,100)$
$\displaystyle 1/\sigma_{b}^{2}$ $\displaystyle\sim\text{Gamma}(0.01,0.01)$
$\displaystyle\alpha_{i}$ $\displaystyle\sim\text{Uniform}(0,2000)$
$\displaystyle\beta_{i}$ $\displaystyle\sim\text{Gaussian}(b,\sigma_{b}^{2})$
$\displaystyle\mu_{ij}$ $\displaystyle=\alpha_{i}+\beta_{i}(j)$ $\displaystyle
1/\sigma^{2}$ $\displaystyle\sim\text{Gamma}(0.01,0.01)$ $\displaystyle
dRT_{ij}$ $\displaystyle\sim\text{Gaussian}(\mu_{ij},\sigma^{2})$ Figure 3:
Graphical model of the hierarchical Bayesian linear regression model for the
numerical distance effect. Following the convention of Lee Wagenmakers (2014),
nodes represent variables of interest (observed=shaded, latent=unshaded), with
dependencies represented via the graph structure.

###  2.1 Fitting the model

####  2.1.1 Data

I fit the model to the RT data from Fazio . (2014). Fifty-five 5th graders
completed a symbolic number comparison task, in which they were asked to
choose which of two Arabic numerals was larger. Each child completed 40 trials
using stimulus numbers ranging from 5 to 21. Of the 40 trials, 10 came from
each of 4 ratio bins: 1.15 - 1.28, 1.28 - 1.43, 1.48 - 1.65, and 2.46 - 2.71.
The ratio for each stimulus pair was defined as the quotient obtained when
dividing the larger number by the smaller number. In all, participants
completed 2,200 trials. I removed 95 error trials and 6 additional trials for
which RT exceeded 5 seconds (a total of 4.6% of trials). The remaining 2,099
trials were collapsed into $55\times 4$ cells by computing median $RT$ for
each of the conditions defined by crossing the factors of subject and ratio
bin.

####  2.1.2 Results

The regression model was fit using the same procedure as in Model 1. As
before, visual inspection of trace plots indicated that all chains converged
appropriately, with $\hat{R}=0.001$ for all parameters.

The posterior distribution of the group-level slope $b$ is depicted in Figure
4. As can be seen in the figure, the mass of the distribution is centered over
the posterior mode $b=-65.9$, with 95% HPDI for the slope $b$ equal to
$[-83.4,-50.7]$. In addition to estimating the group-level slope, we can test
the existence of the numerical distance effect via a Bayes factor, which as
before I computed using the Savage-Dickey density ratio. This Bayes factor was
approximately 1.7 million to 1 in favor of the alternative hypothesis,
indicating overwhelming support for an overall numerical distance effect.

Figure 4: Posterior distribution of the group-level slope $b$ for Model 2,
indicating a substantial numerical distance effect across the group.

##  3 Comparing the model to the classical approach

In this section, I will demonstrate that the hierarchical Bayesian model
developed in this paper can provide better measurement fidelity for the group-
level regression slopes that are desired when assessing group-level SNARC
effects or NDEs, particularly for small samples. As mentioned earlier, the
reason for this increased accuracy is because of a property that is unique to
Bayesian inference – namely the property of _shrinkage_. When one specifies a
group-level prior distribution for the slope (e.g., the
$\text{Uniform}(-20,20)$ distribution that I used in Model 1), this prior is
then combined with the data likelihood via Bayes theorem by multiplication. In
our case, extreme individual slope estimates (those beyond -20 and 20) vanish
by virtue of being multiplied by the prior probability of obtaining those
estimates (i.e., probability = 0). The resulting posterior distribution is
then _shrunk_ away from these parameter estimates.

To demonstrate this, I conducted a simulation. I randomly generated $dRT$
values for a small samples of $n=15$ simulated participants as follows. First,
I generated 15 random slopes $b_{i}$, where $b_{i}\sim\text{Gaussian}(-10,1)$.
That is, I assumed that individual slopes are drawn from a normal distribution
centered at -10, with standard deviation 1. Similarly, I randomly generated 15
random intercepts $a_{i}$, where $a_{i}\sim\text{Uniform}(-200,200)$. Then, I
generated

$dRT_{ij}=a_{i}+b_{i}(j)+\varepsilon_{ij}$

where $i=1,\dots,15$, $j=1,2,8,9$, and
$\varepsilon_{ij}\sim\text{Gaussian}(0,100)$. Then, I computed a classical
linear regression for each participant $i=1,\dots,15$, recording each
estimated slope $\hat{b}_{i}$. Finally, I fit the hierarchical Bayesian model
(Model 1) to the overall set of data to compute the posterior distribution of
$b$.

The results can be seen in Figure 5. Notice that while both the frequentist
95% confidence interval and the Bayesian 95% highest posterior density
interval contain the population mean $b=-10$, the Bayesian estimate is much
less variable. This is because compared to the classical linear regression
method (e.g., Lorch & Myers, 1990), the extreme parameter estimates are shrunk
toward the posterior mean. The same story is repeated with hypothesis testing,
as well. Indeed, performing a $t$-test on the individual regression slopes (as
in Fias, 1996) results in a non-significant SNARC effect, $t(14)=-1.82$,
$p=0.09$, which results in a Type II error. However, a Bayesian hypothesis
test (the Savage-Dickey method) yields a Bayes factor that favors the SNARC
effect by a factor of 75 to 1.

Figure 5: Slope estimates from a simulated data set of $n=15$ participants.
The solid line represents the posterior density of $b$, the group-level slope
estimate for the SNARC effect. The histogram represents the slope estimates
obtained from separate classical linear regressions for each individual.

##  4 Discussion

The purpose of this paper was to introduce a hierarchical Bayesian linear
regression framework for measuring group-level and individual-level numerical
representations. The models were then applied to two common markers of
numerical representations: the SNARC effect, which indexes spatial-numerical
associations, and the numerical distance effect (NDE), which indexes
representations of numerical magnitude. The models developed here represent a
straightforward Bayesian extension of the classical linear regression approach
to measuring these effects that was introduced by Lorch Myers (1990) and first
applied to the SNARC effect by Fias (1996).

Both models were built from a hierarchical Bayesian framework, which is
advantageous for several reasons. First, the models assume that the regression
slopes (the primary object of interest in these analyses) are drawn from a
group-level distribution. This hierarchical definition allows both group-level
and individual-level slope estimates to be modeled simultaneously from the
data. As such, we can answer questions at either the group-level (i.e., is
there a SNARC effect?) or participant-level (i.e., are the slopes associated
with ohter measures, such as math achievement?). For example, one might use
this framework as a stepping stone to a more complex model, where the group-
level slopes $b$ are hypothesized to differ by some independent variable. A
concrete instance of this comes from studies in which the authors found group
differences in the SNARC effect (e.g., Fischer ., 2010; Cipora ., 2015). One
could perform similar studies by building an additional linear model on these
group-level parameters, with priors appropriate to such effects (e.g., a
Cauchy prior, as in Rouder ., 2009). Importantly, since the inference would be
done in a Bayesian framework, it is possible to measure evidence for null
effects too (Wagenmakers, 2007), so the model could be used to test for
_invariances_ as well as differences.

Another advantage to using a hierarchical Bayesian model for numerical
representations is that such models tend to have better measurement fidelity.
Indeed, one of the advantages of the Bayesian framework in general is the
notion of shrinkage, where extreme parameter estimates are “shrunk” toward the
mean by virtue of the prior (Gelman ., 2013). This property is particularly
salient for small samples, where classical frequentist methods tend to perform
poorly. I demonstrated exactly this phenomenon in the simulation above:
because of some extreme individual-level slope estimates, the frequentist
confidence interval was quite wide, and as a result, we could not detect the
SNARC effect. However, the Bayesian highest posterior density interval
provided a much more accurate estimate of the true population slope.
Critically, in this simulation, only the Bayesian method produced the correct
inference.

A final advantage to this modeling framework that I will mention is that all
assumptions of the model are made explicit. This may be a new approach to
some, especially those who are accustomed to classical inferential software
packages for which the statistical assumptions are kept “under the hood.”
However, I think the approach presented in this paper can be very useful to a
wide variety of problems in cognitive psychology, as the researcher can take
the models presented here and modify them for any desired context. Indeed, the
prior knowledge of a given field can be easily and coherently integrated into
the model without too much work.

In summary, the models developed in this paper will provide a flexible tool
that can be used to estimate group-level and individual-level numerical
representations quickly and coherently. Further, the methods developed are not
specific to numerical cognition, so application to a wide variety of problems
should not be too difficult. Regardless of the application, the hierarchical
Bayesian linear regression framework presented here provides a powerful,
coherent, and accurate measurement model for applied work in cognitive
psychology.

## References

  * Cipora . (2015) cipora2015Cipora, K., Hohol, M., Nuerk, HC., Willmes, K., Brożek, B., Kucharzyk, B. Necka, E.  2015\.  Professional mathematicians differ from controls in their spatial-numerical associations Professional mathematicians differ from controls in their spatial-numerical associations.  Psychological Research804710–726.  10.1007/s00426-015-0677-6
  * Dehaene (1992) dehaene1992Dehaene, S.  1992\.  Varieties of numerical abilities Varieties of numerical abilities.  Cognition441-21–42.  10.1016/0010-0277(92)90049-n
  * Dehaene . (1990) dehaene1990Dehaene, S., Dupoux, E. Mehler, J.  1990\.  Is numerical comparison digital? Analogical and symbolic effects in two-digit number comparison. Is numerical comparison digital? Analogical and symbolic effects in two-digit number comparison.  Journal of Experimental Pychology: Human Perception and Performance163626–641.  10.1037/0096-1523.16.3.626
  * Fazio . (2014) fazio2014Fazio, LK., Bailey, DH., Thompson, CA. Siegler, RS.  2014\.  Relations of different types of numerical magnitude representations to each other and to mathematics achievement Relations of different types of numerical magnitude representations to each other and to mathematics achievement.  Journal of Experimental Child Psychology12353–72.  10.1016/j.jecp.2014.01.013
  * Fias (1996) fias1996Fias, W.  1996\.  The Importance of Magnitude Information in Numerical Processing: Evidence from the SNARC Effect The importance of magnitude information in numerical processing: Evidence from the SNARC effect.  Mathematical Cognition2195–110.  10.1080/135467996387552
  * Fischer . (2010) fischer2010Fischer, MH., Mills, RA. Shaki, S.  2010\.  How to cook a SNARC: Number placement in text rapidly changes spatial–numerical associations How to cook a SNARC: Number placement in text rapidly changes spatial–numerical associations.  Brain and Cognition723333–336.  10.1016/j.bandc.2009.10.010
  * Gelman . (2013) gelman2013Gelman, A., Carlin, JB., Stern, HS. Rubin, DB.  2013\.  Bayesian data analysis Bayesian data analysis (3rd ).  Chapman & Hall/CRC, Boca Raton, FL.
  * Gelman Hill (2013) gelmanHill2013Gelman, A. Hill, J.  2013\.  Data Analysis Using Regression and Multilevel/Hierarchical Models Data analysis using regression and multilevel/hierarchical models.  Cambridge University Press.
  * Gelman Rubin (1992) gelman1992Gelman, A. Rubin, DB.  1992\.  Inference from Iterative Simulation Using Multiple Sequences Inference from iterative simulation using multiple sequences.  Statistical Science74457–472.  10.1214/ss/1177011136
  * Gevers . (2006) gevers2006Gevers, W., Verguts, T., Reynvoet, B., Caessens, B. Fias, W.  2006\.  Numbers and space: A computational model of the SNARC effect. Numbers and space: A computational model of the SNARC effect.  Journal of Experimental Psychology: Human Perception and Performance32132–44.  10.1037/0096-1523.32.1.32
  * Haider Frensch (2002) haider2002Haider, H. Frensch, PA.  2002\.  Why aggregated learning follows the power law of practice when individual learning does not: Comment on Rickard (1997, 1999), Delaney et al. (1998), and Palmeri (1999). Why aggregated learning follows the power law of practice when individual learning does not: Comment on Rickard (1997, 1999), Delaney et al. (1998), and Palmeri (1999).  Journal of Experimental Psychology: Learning, Memory, and Cognition282392–406.  10.1037/0278-7393.28.2.392
  * Heathcote . (2000) heathcote2000Heathcote, A., Brown, S. Mewhort, DJK.  2000\.  The power law repealed: The case for an exponential law of practice The power law repealed: The case for an exponential law of practice.  Psychonomic Bulletin & Review72185–207.  10.3758/bf03212979
  * Hoekstra . (2014) hoekstra2014Hoekstra, R., Morey, RD., Rouder, JN. Wagenmakers, EJ.  2014\.  Robust misinterpretation of confidence intervals Robust misinterpretation of confidence intervals.  Psychonomic Bulletin & Review2151157–1164.  10.3758/s13423-013-0572-3
  * Holloway Ansari (2009) holloway2009Holloway, ID. Ansari, D.  2009\.  Mapping numerical magnitudes onto symbols: The numerical distance effect and individual differences in children’s mathematics achievement Mapping numerical magnitudes onto symbols: The numerical distance effect and individual differences in children’s mathematics achievement.  Journal of Experimental Child Psychology103117–29.  10.1016/j.jecp.2008.04.001
  * Kruschke (2015) kruschke2015Kruschke, JK.  2015\.  Doing Bayesian Data Analysis Doing Bayesian data analysis (2nd ).  Academic Press, London.
  * Lee Wagenmakers (2014) lee2014Lee, MD. Wagenmakers, EJ.  2014\.  Bayesian cognitive modeling: A practical course Bayesian cognitive modeling: A practical course.  Cambridge university press.
  * Lorch Myers (1990) lorch1990Lorch, RF. Myers, JL.  1990\.  Regression analyses of repeated measures data in cognitive research. Regression analyses of repeated measures data in cognitive research.  Journal of Experimental Psychology: Learning, Memory, and Cognition161149–157.  10.1037/0278-7393.16.1.149
  * McElreath (2015) mcelreath2015McElreath, R.  2015\.  Statistical Rethinking Statistical rethinking.  Chapman & Hall/CRC, Boca Raton, FL.
  * Moyer Landauer (1967) moyer1967Moyer, RS. Landauer, TK.  1967\.  Time required for judgements of numerical inequality Time required for judgements of numerical inequality.  Nature2151519-1520.  10.1038/2151519a0
  * Pinhas . (2014) pinhas2014Pinhas, M., Shaki, S. Fischer, MH.  2014\.  Heed the signs: Operation signs have spatial associations Heed the signs: Operation signs have spatial associations.  The Quarterly Journal of Experimental Psychology6781527–1540.  10.1080/17470218.2014.892516
  * Plummer (2017) jagsPlummer, M.  2017\.  JAGS: Just Another Gibbs Sampler JAGS: Just Another Gibbs Sampler [].  https://mcmc-jags.sourceforge.net
  * Pratte . (2010) pratte2010Pratte, MS., Rouder, JN. Morey, RD.  2010\.  Separating mnemonic process from participant and item effects in the assessment of ROC asymmetries. Separating mnemonic process from participant and item effects in the assessment of ROC asymmetries.  Journal of Experimental Psychology: Learning, Memory, and Cognition361224–232.  10.1037/a0017682
  * R Core Team (2016) rR Core Team.  2016\.  R: A Language and Environment for Statistical Computing R: A language and environment for statistical computing [].  Vienna, Austria.  https://www.R-project.org/
  * Rouder . (2009) rouder2009Rouder, JN., Speckman, PL., Sun, D., Morey, RD. Iverson, G.  2009\.  Bayesian $t$ tests for accepting and rejecting the null hypothesis Bayesian $t$ tests for accepting and rejecting the null hypothesis.  Psychonomic Bulletin & Review162225–237.  10.3758/pbr.16.2.225
  * Viarouge . (2014) viarouge2014Viarouge, A., Hubbard, EM. McCandliss, BD.  2014\.  The Cognitive Mechanisms of the SNARC Effect: An Individual Differences Approach The cognitive mechanisms of the SNARC effect: An individual differences approach.  PLoS ONE94e95756.  10.1371/journal.pone.0095756
  * Wagenmakers (2007) wagenmakers2007Wagenmakers, EJ.  2007\.  A practical solution to the pervasive problems ofp values A practical solution to the pervasive problems ofp values.  Psychonomic Bulletin & Review145779–804.  https://doi.org/10.3758/bf03194105 10.3758/bf03194105
  * Wagenmakers . (2010) wagenmakers2010Wagenmakers, EJ., Lodewyckx, T., Kuriyal, H. Grasman, R.  2010\.  Bayesian hypothesis testing for psychologists: A tutorial on the Savage-Dickey method Bayesian hypothesis testing for psychologists: A tutorial on the Savage-Dickey method.  Cognitive Psychology603158–189.  10.1016/j.cogpsych.2009.12.001
  * Wood . (2008) wood2008Wood, G., Willmes, K., Nuerk, HC. Fischer, MH.  2008\.  On the cognitive link between space and number: A meta-analysis of the SNARC effect On the cognitive link between space and number: A meta-analysis of the SNARC effect.  Psychology Science Quarterly504489–525.