#  Study of $I=0$ bottomonium bound states and resonances in $S$, $P$, $D$ and
$F$ waves with lattice QCD static-static-light-light potentials

(1)Pedro Bicudo<EMAIL_ADDRESS>(1)Nuno Cardoso
<EMAIL_ADDRESS>(2)Lasse Mueller<EMAIL_ADDRESS>frankfurt.de (2),(3)Marc Wagner<EMAIL_ADDRESS>(1)CeFEMA, Dep.
Física, Instituto Superior Técnico, Universidade de Lisboa, Av. Rovisco Pais,
1049-001 Lisboa, Portugal  (2)Johann Wolfgang Goethe-Universität Frankfurt am
Main, Institut für Theoretische Physik, Max-von-Laue-Straße 1, D-60438
Frankfurt am Main, Germany  (3)Helmholtz Research Academy Hesse for FAIR,
Campus Riedberg,
Max-von-Laue-Straße 12, D-60438 Frankfurt am Main, Germany

###### Abstract

In this paper we study $I=0$ bottomonium in $S$, $P$, $D$ and $F$ waves
considering five coupled channels, one confined quarkonium and four open
$B^{(*)}\bar{B}^{(*)}$ and $B_{s}^{(*)}\bar{B}_{s}^{(*)}$ meson-meson
channels. To this end we use and extend a recently developed novel approach
utilizing lattice QCD string breaking potentials for the study of quarkonium
bound states and resonances. This approach is based on the Born Oppenheimer
approximation and the unitary emergent wave method and allows to compute the
poles of the T matrix. We compare our results to existing experimental results
for $I=0$ bottomonium and discuss masses, decay widths and the assignment of
angular momentum quantum numbers. Moreover, we determine the quarkonium and
meson-meson composition of these states to clarify, which of them are ordinary
quarkonium, and which of them should rather be interpreted as tetraquarks.

##  I Introduction

In the last decade a whole new class of hadrons has been discovered
experimentally, so-called tetraquarks, which are composed of two quarks and
two antiquarks [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18,
19, 20, 21]. In particular, a larger number of such tetraquarks observed in
the last couple of years at Belle, BESIII and LHCb have at least one heavy
quark, as anticipated in Refs. [22, 23]. From the onset of QCD the existence
of tetraquarks was expected [24], However, a quantitative first principles
prediction of their properties, e.g. quark composition, quantum numbers,
masses and decay widths, using e.g. lattice QCD, remains to be achieved. The
main reason, why this has not been successful yet, is that the majority of
observed heavy tetraquarks are resonances high in the quarkonium spectrum (the
only exception is the $T_{cc}$ tetraquark recently found by LHCb at CERN [25,
26]). Studying such resonances with lattice QCD is possible in principle using
the Lüscher phase shift method [27], but practically feasible only for a
single or a small number of decay channels. For several open channels, as it
is the case for some of the recently observed heavy tetraquarks, following
this path seems tremendously difficult.

Because of these difficulties we recently started to develop another approach
[28, 29], utilizing lattice QCD results for static potentials, to study
bottomonium as well as tetraquark resonances with the same non-exotic quantum
numbers high in the spectrum. We start with lattice QCD potentials computed
with static quarks and light quarks in the context of string breaking [30, 31,
32], which provide information on the interactions between a two-quark
quarkonium channel and several four-quark meson-meson channels. We use the
Born-Oppenheimer diabatic approximation as in Refs. [33, 34], i.e. include the
kinetic energy of the heavy quarks, and study the dynamically coupled
quarkonium and meson-meson channels with techniques from quantum mechanics.

In our previous work [28, 29] we applied this method to systems with an $S$
wave bottonomium channel coupled to a $B^{(*)}\bar{B}^{(*)}$ and a
$B_{s}^{(*)}\bar{B}_{s}^{(*)}$ channel. In this work we extend our studies to
$P$ wave, $D$ wave and $F$ wave quarkonium channels, again coupled to
$B^{(*)}\bar{B}^{(*)}$ and $B_{s}^{(*)}\bar{B}_{s}^{(*)}$ channels. Because of
the intrinsic parity of quarks, the relative orbital angular momentum of our
meson-meson channels differs by one unit from the quarkonium orbital angular
momentum. Compared to the $S$ wave case, this doubles the number of meson-
meson channels, and, thus, we have to take into account five coupled channels.
As mentioned in the previous paragraph, such a large number of coupled
channels seems currently inaccessible for fully dynamical lattice QCD studies
of resonances.

Using our approach we were already able to address claims in the literature,
from studies with different hadronic models, on the nature of some of the
bottonomium resonances observed at Belle. For example, it has been discussed,
which of the observed excited bottonomium resonances are $S$ wave or $D$ wave
states [35, 36, 37], since the total angular momentum is not yet
experimentally determined, and whether they have a large meson-meson content
[38, 39, 40, 41]. Of particular interest is the nature of the newly discovered
resonance $\Upsilon(10753)$ recently observed at Belle with mass $(10.753\pm
0.007)\,\text{GeV}$ [42]. Model and effective field theory calculations
suggest for instance this resonance to be either a tetraquark [43, 44], a
hybrid meson [45, 46, 47] or the more canonical and so far missing
$\Upsilon(3D)$ [35, 36, 37]. With our lattice QCD based approach we found a
pole of the T matrix in the $S$ wave channel corresponding to the mass
$10.773\,\text{GeV}$ [29], very similar to the Belle measurement. Moreover, we
found a large meson-meson component for that state. Thus, we proposed that the
recently observed $\Upsilon(10753)$ is a dynamical state, composed mostly of a
meson-meson pair.

By extending our work to $P$ wave, $D$ wave and $F$ wave channels, we can now
check, whether there is also a candidate for the $\Upsilon(10753)$ in the $D$
wave channel. Moreover, we will obtain a complete picture of the $I=0$
bottomonium spectrum. In particular we can identify states with a large meson-
meson component, which can be interpreted as tetraquarks. Finally, it is
interesting to explore the possible existence of a bottomonium counterpart of
$X(3872)$ discovered at Belle and CDF [48, 49], which is extremely close to
the $D\bar{D}^{*}$ and $D^{*}\bar{D}$ thresholds. In our case, this would
correspond to a resonance very close to either the $B^{(*)}\bar{B}^{(*)}$ or
the $B_{s}^{(*)}\bar{B}_{s}^{(*)}$ thresholds.

This paper is organized as follows. In Section II we detail our approach,
where we use lattice QCD potentials computed for string breaking in a coupled
channel Schrödinger equation. In particular we carry out a partial wave
decomposition and derive $5\times 5$ Schrödinger equations for the $P$ wave,
$D$ wave and $F$ wave channels. Using the emergent wave method, we also show,
how the T matrix can be computed, leading to both the phase shifts and, after
a pole search, to masses and decay widths. In Section III we show and discuss
corresponding numerical results. We also determine the quarkonium and meson-
meson composition for all states. Finally, in Section IV, we conclude on the
points raised in this introduction and discuss possibilities for future
research within our formalism.

##  II Quarkonium resonances from lattice QCD static potentials

###  II.1 Quantum numbers

We consider a heavy quark-antiquark pair and either no light quarks
($\bar{Q}Q$) or a light quark-antiquark pair with isospin $I=0$
($\bar{Q}Q(\bar{u}u+\bar{d}d)\equiv\bar{M}M$ or
$\bar{Q}Q\bar{s}s\equiv\bar{M}_{s}M_{s}$). In the limit, where the heavy
quarks are static, such a system can be characterized by the following quantum
numbers:

  * •

$J^{PC}$: total angular momentum, parity and charge conjugation.

  * •

$S^{PC}_{Q}$: spin of the heavy quark-antiquark pair and corresponding parity
and charge conjugation.

  * •

${\tilde{J}}^{PC}$: total angular momentum excluding the heavy quark spins and
corresponding parity and charge conjugation. (For quarkonium
${\tilde{J}}^{PC}$ coincides with the orbital angular momentum $L^{PC}$ of the
heavy quark-antiquark pair.)

Moreover, the majority of observables, in particular energy levels, do not
depend on $S^{PC}_{Q}$. Thus, the relevant quantum numbers in our context are
${\tilde{J}}^{PC}$ and not, as usual, $J^{PC}$.

For rather heavy $b$ quarks we expect that our approach, which is based on
static symmetries and quantum numbers, will yield reasonably accurate results,
possibly even for $c$ quarks.

###  II.2 The coupled channel Schrödinger equation

Now we discuss the coupled channel Schrödinger equation, which we have derived
in detail in our previous papers [28, 29]. We take the two lowest meson decay
channels into account, where each channel contains two negative parity heavy-
light mesons, either $\bar{M}M$ or $\bar{M}_{s}M_{s}$. The corresponding light
spin is $S_{q}^{PC}=1^{--}$ [28]. In addition to these decay channels there
is, of course, also the quarkonium channel $\bar{Q}Q$. This amounts to a
seven-component wave function
$\psi(\textbf{r})=(\psi_{\bar{Q}Q}(\textbf{r}),\vec{\psi}_{\bar{M}M}(\textbf{r}),\vec{\psi}_{\bar{M}_{s}M_{s}}(\textbf{r}))$,
where the first component represents the $\bar{Q}Q$ channel and the remaining
six components the spin-1 triplets of the $\bar{M}M$ and the
$\bar{M}_{s}M_{s}$ channel, respectively. r denotes the relative coordinates
of the heavy quark-antiquark pair.

The coupled channel Schrödinger equation reads

$\displaystyle\Bigg{(}-\frac{1}{2}\mu^{-1}\bigg{(}\partial_{r}^{2}+\frac{2}{r}\partial_{r}-\frac{\textbf{L}^{2}}{r^{2}}\bigg{)}+V(\textbf{r})+\left(\begin{array}[]{ccc}E_{\text{threshold}}&0&0\\\
0&2m_{M}&0\\\ 0&0&2m_{M_{s}}\\\
\end{array}\right)-E\Bigg{)}\psi(\textbf{r})=0,$ (4)

where
$\mu^{-1}=\textrm{diag}(1/\mu_{Q},1/\mu_{M},1/\mu_{M},1/\mu_{M},1/\mu_{M_{s}},1/\mu_{M_{s}},1/\mu_{M_{s}})$
is a $7\times 7$ diagonal matrix with the reduced masses of the heavy quarks
and the heavy mesons, $\mu_{Q}=m_{Q}/2$, $\mu_{M}=m_{M}/2$ and
$\mu_{M_{s}}=m_{M_{s}}/2$. In the static limit the pseudoscalar and the vector
heavy-light meson masses are identical. For finite heavy quark mass there is a
small difference in these meson masses, e.g. for heavy $b$ quarks
$m_{B^{*}}-m_{B}\approx 45\,\text{MeV}$ and
$m_{B_{s}^{*}}-m_{B_{s}}=49\,\text{MeV}$. We take the spin averaged masses for
$m_{M}$ and $m_{M_{s}}$ (see Section III.2).
$\textbf{L}=\textbf{r}\times\textbf{p}$ denotes the orbital angular momentum
operator and $E_{\text{threshold}}$ is the threshold energy corresponding to
two negative parity static-light mesons in the same lattice setup, where the
static potentials are computed (for more details see Section III.2 and Ref.
[29]). The potential matrix $V(\textbf{r})$ is given by

$\displaystyle
V(\textbf{r})=\left(\begin{array}[]{ccc}V_{\bar{Q}Q}(r)&V_{\textrm{mix}}(r)\left(1\otimes\textbf{e}_{r}\right)&(1/\sqrt{2})V_{\textrm{mix}}(r)\left(1\otimes\textbf{e}_{r}\right)\\\
V_{\textrm{mix}}(r)\left(\textbf{e}_{r}\otimes 1\right)&V_{\bar{M}M}(r)&0\\\
(1/\sqrt{2})V_{\textrm{mix}}(r)\left(\textbf{e}_{r}\otimes
1\right)&0&V_{\bar{M}M}(r)\\\ \end{array}\right)$ (8)

with

$\displaystyle
V_{\bar{M}M}(r)=V_{\bar{M}M,\parallel}(r)\Big{(}\textbf{e}_{r}\otimes\textbf{e}_{r}\Big{)}+V_{\bar{M}M,\perp}(r)\Big{(}1-\textbf{e}_{r}\otimes\textbf{e}_{r}\Big{)},$
(9)

where we have assumed that meson-meson interactions vanish. $V_{\bar{Q}Q}(r)$,
$V_{\bar{M}M,\parallel}(r)$, $V_{\bar{M}M,\perp}(r)$ and $V_{\textrm{mix}}(r)$
can be expressed in terms of QCD static potentials, which can be computed with
lattice QCD (see Refs. [30, 28] for details). $V_{\bar{Q}Q}(r)$ represents the
potential of a heavy quark-antiquark pair, $V_{\bar{M}M,\parallel}(r)$,
$V_{\bar{M}M,\perp}(r)$ the interaction between a pair of heavy-light mesons
and $V_{\textrm{mix}}(r)$ describes the mixing of the quarkonium channel and
the meson-meson channels. We use the same potentials both for the $\bar{M}M$
channel and the $\bar{M}_{s}M_{s}$ channel and expect this to be reasonable,
because the light quark dependence of static potentials is known to be rather
mild. We confirmed this expectation by carrying out a consistency check with a
recent $2+1$-flavor lattice study of string breaking [31] (for more details
see our previous work [29]). Note that, the mixing between the $\bar{Q}Q$
channel and the $\bar{M}_{s}M_{s}$ channel is suppressed by $1/\sqrt{2}$,
because there is only a single strange quark flavor in comparison to the two
degenerate light flavors corresponding to the $\bar{M}M$ channel.

We use lattice QCD data from Ref. [30] to determine continuous functions
$V_{\bar{Q}Q}(r)$, $V_{\bar{M}M,\parallel}(r)$, $V_{\bar{M}M,\perp}(r)$ and
$V_{\textrm{mix}}(r)$. The data points for $V_{\bar{Q}Q}(r)$ and
$V_{\textrm{mix}}(r)$ are consistently parameterized by

$\displaystyle V_{\bar{Q}Q}(r)=E_{0}-\frac{\alpha}{r}+\sigma
r+\sum_{j=1}^{2}c_{\bar{Q}Q,j}\,r\exp\bigg{(}-\frac{r^{2}}{2\lambda^{2}_{\bar{Q}Q,j}}\bigg{)}$
(10) $\displaystyle
V_{\textrm{mix}}(r)=\sum_{j=1}^{2}c_{\textrm{mix},j}\,r\exp\bigg{(}-\frac{r^{2}}{2\lambda^{2}_{\textrm{mix},j}}\bigg{)}$
(11)

with parameters $E_{0}$, $\alpha$, $\sigma$, $c_{\bar{Q}Q,j}$,
$\lambda_{\bar{Q}Q,j}$, $c_{\textrm{mix},j}$, $\lambda_{\textrm{mix},j}$
collected in Table 1. The data points for $V_{\bar{M}M,\parallel}(r)$ are
consistent with

$\displaystyle V_{\bar{M}M,\parallel}(r)=0.$ (12)

For $V_{\bar{M}M,\perp}(r)$ no lattice data is available yet. We assume

$\displaystyle V_{\bar{M}M,\perp}(r)=0$ (13)

(see also Ref. [28]). In Fig. 1 we show the data points for $V_{\bar{Q}Q}(r)$,
$V_{\bar{M}M,\parallel}(r)$ and $V_{\textrm{mix}}(r)$ together with the
parameterizations (10) to (12).

potential | parameter | value
---|---|---
$V_{\bar{Q}Q}(r)$ | $E_{0}$ | $-1.599(269)\,\textrm{GeV}\phantom{1.^{-1}}$
| $\alpha$ | $+0.320(94)\phantom{1.0\,\textrm{GeV}^{-1}}$
| $\sigma$ | $+0.253(035)\,\textrm{GeV}^{2\phantom{-}}\phantom{1.}$
| $c_{\bar{Q}Q,1}$ | $+0.826(882)\,\textrm{GeV}^{2\phantom{-}}\phantom{1.}$
| $\lambda_{\bar{Q}Q,1}$ | $+0.964(47)\,\textrm{GeV}^{-1}\phantom{1.0}$
| $c_{\bar{Q}Q,2}$ | $+0.174(1.004)\,\textrm{GeV}^{2\phantom{-}}$
| $\lambda_{\bar{Q}Q,2}$ | $+2.663(425)\,\textrm{GeV}^{-1}\phantom{1.}$
$V_{\textrm{mix}}(r)$ | $c_{\textrm{mix},1}$ | $-0.988(32)\,\textrm{GeV}^{2\phantom{-}}\phantom{1.0}$
| $\lambda_{\textrm{mix},1}$ | $+0.982(18)\,\textrm{GeV}^{-1}\phantom{1.0}$
| $c_{\textrm{mix},2}$ | $-0.142(7)\,\textrm{GeV}^{2\phantom{-}}\phantom{1.00}$
| $\lambda_{\textrm{mix},2}$ | $+2.666(46)\,\textrm{GeV}^{-1}\phantom{1.0}$
Table 1: Parameters of the potential parameterizations (10) and (11). Figure
1: Data points for $V_{\bar{Q}Q}(r)$, $V_{\bar{M}M,\parallel}(r)$ and
$V_{\textrm{mix}}(r)$. The curves represent the parameterizations (10) to
(12).

###  II.3 The coupled channel Schrödinger equation for total angular momentum
${\tilde{J}}$

Now we specialize the coupled channel Schrödinger equation (4) to study
quarkonium bound states and resonances and meson-meson scattering for
arbitrary given ${\tilde{J}}$. For details we refer to Section II D of Ref.
[28], where the mathematical formalism is discussed for ${\tilde{J}}=0$.

At first, each of the meson-meson components of the wave function
$\psi(\textbf{r})$ is written as sum of an incident wave, which is a solution
of the free Schrödinger equation, i.e. Eq. (4) with $V(\textbf{r})=0$, and an
emergent wave. This will allow us to define and determine scattering
amplitudes for energies above the lowest meson-meson threshold. In contrast to
our previous work [28], we do not restrict the scattering problem to an
incident plane wave, but allow arbitrary regular solutions of the free
Schrödinger equation, including superpositions of $\bar{M}M$ and
$\bar{M}_{s}M_{s}$ waves. We expand
$\psi(\textbf{r})=(\psi_{\bar{Q}Q}(\textbf{r}),\vec{\psi}_{\bar{M}M}(\textbf{r}),\vec{\psi}_{\bar{M}_{s}M_{s}}(\textbf{r}))$
in terms of eigenfunctions of $\tilde{\textbf{J}}^{2}$ and ${\tilde{J}}_{z}$,

$\displaystyle\psi_{QQ}(\textbf{r})$
$\displaystyle=\frac{u_{0,0}(r)}{r}Y_{0,0}(\Omega)+\sum_{\bar{J}=1}^{\infty}\sum_{{\tilde{J}}_{z}=-{\tilde{J}}}^{+{\tilde{J}}}\frac{u_{{\tilde{J}},{\tilde{J}}_{z}}(r)}{r}Y_{{\tilde{J}},{\tilde{J}}_{z}}(\Omega)$
(14) $\displaystyle\vec{\psi}_{\bar{M}_{(s)}M_{(s)}}(\textbf{r})$
$\displaystyle=\alpha_{\bar{M}_{(s)}M_{(s)},0,0}\,j_{0}(k_{(s)}r)\,\textbf{Z}_{1\rightarrow
0,0}(\Omega)+\sum_{\bar{J}=1}^{\infty}\sum_{{\tilde{J}}_{z}=-{\tilde{J}}}^{+{\tilde{J}}}\sum_{L={\tilde{J}}-1,{\tilde{J}},{\tilde{J}}+1}\alpha_{\bar{M}_{(s)}M_{(s)},{\tilde{J}},{\tilde{J}}_{z}}j_{L}(k_{(s)}r)\textbf{Z}_{L\rightarrow{\tilde{J}},{\tilde{J}}_{z}}(\Omega)$
(15) $\displaystyle+\frac{\chi_{\bar{M}_{(s)}M_{(s)},1\rightarrow
0,0}}{r}\textbf{Z}_{1\rightarrow
0,0}(\Omega)+\sum_{\bar{J}=1}^{\infty}\sum_{{\tilde{J}}_{z}=-{\tilde{J}}}^{+{\tilde{J}}}\sum_{L={\tilde{J}}-1,{\tilde{J}},{\tilde{J}}+1}\frac{\chi_{\bar{M}_{(s)}M_{(s)},L\rightarrow{\tilde{J}},{\tilde{J}}_{z}}}{r}\textbf{Z}_{L\rightarrow{\tilde{J}},{\tilde{J}}_{z}}(\Omega)$
(16)

where $Y_{{\tilde{J}},{\tilde{J}}_{z}}(\Omega)$ are the spherical harmonics
and a detailed definition of
$\textbf{Z}_{L\rightarrow{\tilde{J}},{\tilde{J}}_{z}}(\Omega)$ can be found in
our previous work [28].
$\alpha_{\bar{M}_{(s)}M_{(s)},{\tilde{J}},{\tilde{J}}_{z}}$ denotes the
expansion coefficient for an incoming $\bar{M}_{(s)}M_{(s)}$-wave with angular
momentum ${\tilde{J}}$ and $j_{L}(k_{(s)}r)$ denote spherical Bessel
functions.  The Schrödinger equation (4) can then be projected in a
straightforward way to definite ${\tilde{J}}$,

$\displaystyle\left(\frac{1}{2}\mu^{-1}\bigg{(}\partial_{r}^{2}+\frac{1}{r^{2}}L_{{\tilde{J}}}^{2}\bigg{)}+V_{{\tilde{J}}}(r)+\left(\begin{array}[]{ccccc}E_{\text{threshold}}&0&0&0&0\\\
0&2m_{M}&0&0&0\\\ 0&0&2m_{M}&0&0\\\ 0&0&0&2m_{M_{s}}&0\\\
0&0&0&0&2m_{M_{s}}\\\
\end{array}\right)-E\right)\left(\begin{array}[]{c}u_{{\tilde{J}}}(r)\\\
\chi_{\bar{M}M,{\tilde{J}}-1\rightarrow{\tilde{J}}}(r)\\\
\chi_{\bar{M}M,{\tilde{J}}+1\rightarrow{\tilde{J}}}(r)\\\
\chi_{\bar{M}_{s}M_{s},{\tilde{J}}-1\rightarrow{\tilde{J}}}(r)\\\
\chi_{\bar{M}_{s}M_{s},{\tilde{J}}+1\rightarrow{\tilde{J}}}(r)\end{array}\right)=$
(27) $\displaystyle\quad\quad=\left(\begin{array}[]{c}V_{\textrm{mix}}(r)\\\
0\\\ 0\\\ 0\\\
0\end{array}\right)\bigg{(}\alpha_{\bar{M}M,{\tilde{J}}-1}{{\tilde{J}}\over
2{\tilde{J}}+1}rj_{{\tilde{J}}-1}(kr)+\alpha_{\bar{M}M,{\tilde{J}}+1}{{\tilde{J}}+1\over
2{\tilde{J}}+1}rj_{{\tilde{J}}+1}(kr)$ (33)
$\displaystyle\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad+\alpha_{\bar{M}_{s}M_{s},{\tilde{J}}-1}{{\tilde{J}}\over
2{\tilde{J}}+1}\frac{rj_{{\tilde{J}}-1}(k_{s}r)}{\sqrt{2}}+\alpha_{\bar{M}_{s}M_{s},{\tilde{J}}+1}{{\tilde{J}}+1\over
2{\tilde{J}}+1}\frac{rj_{{\tilde{J}}+1}(k_{s}r)}{\sqrt{2}}\bigg{)}$ (34)

with
$\mu^{-1}=\textrm{diag}(1/\mu_{Q},1/\mu_{M},1/\mu_{M},1/\mu_{M_{s}},1/\mu_{M_{s}})$,
$L_{\tilde{J}}^{2}=\textrm{diag}({\tilde{J}}({\tilde{J}}+1),({\tilde{J}}-1){\tilde{J}},({\tilde{J}}+1)({\tilde{J}}+2),({\tilde{J}}-1){\tilde{J}},({\tilde{J}}+1)({\tilde{J}}+2))$
and

$\displaystyle V_{\tilde{J}}(r)$
$\displaystyle=\left(\begin{array}[]{ccccc}V_{\bar{Q}Q}&\sqrt{{{\tilde{J}}\over
2{\tilde{J}}+1}}V_{\textrm{mix}}&\sqrt{{{\tilde{J}}+1\over
2{\tilde{J}}+1}}V_{\textrm{mix}}&{1\over\sqrt{2}}\sqrt{{{\tilde{J}}\over
2{\tilde{J}}+1}}V_{\textrm{mix}}&{1\over\sqrt{2}}\sqrt{{{\tilde{J}}+1\over
2{\tilde{J}}+1}}V_{\textrm{mix}}\\\ \sqrt{{{\tilde{J}}\over
2{\tilde{J}}+1}}V_{\textrm{mix}}&0&0&0&0\\\ \sqrt{{{\tilde{J}}+1\over
2{\tilde{J}}+1}}V_{\textrm{mix}}&0&0&0&0\\\
{1\over\sqrt{2}}\sqrt{{{\tilde{J}}\over
2{\tilde{J}}+1}}V_{\textrm{mix}}&0&0&0&0\\\
{1\over\sqrt{2}}\sqrt{{{\tilde{J}}+1\over
2{\tilde{J}}+1}}V_{\textrm{mix}}&0&0&0&0\\\ \end{array}\right).$ (40)

Note that, equation (34) is degenerate with respect to ${\tilde{J}}_{z}$ thus
we dropped the index. The confining $\bar{Q}Q$ channel is represented by the
radial wave function $u_{{\tilde{J}}}(r)$ with the following boundary
conditions:

  * •

For $r\rightarrow 0$:

$\displaystyle u_{{\tilde{J}}}(r)\propto r^{{\tilde{J}}+1}.$ (41)

  * •

For $r\rightarrow\infty$:

$\displaystyle u_{{\tilde{J}}}(r)=0.$ (42)

The incident wave becomes a superposition of spherical waves represented by
Bessel functions $j_{L_{\textrm{in}}}$. These include $\bar{M}M$ waves with
$L_{\textrm{in}}={\tilde{J}}-1$ as well as $L_{\textrm{in}}={\tilde{J}}+1$ and
$\bar{M}_{s}M_{s}$ waves with $L_{\textrm{in}}={\tilde{J}}-1$ as well as
$L_{\textrm{in}}={\tilde{J}}+1$, where $L_{\textrm{in}}$ denotes orbital
angular momentum. Spherical waves with $L_{\textrm{in}}={\tilde{J}}$ are
excluded, because of parity. For example, an incident $\bar{M}M$ wave with
$L_{\textrm{in}}={\tilde{J}}-1$ translates to
$\vec{\alpha}=(\alpha_{\bar{M}M,{\tilde{J}}-1},\alpha_{\bar{M}M,{\tilde{J}}+1},\alpha_{\bar{M}_{s}M_{s},{\tilde{J}}-1},\alpha_{\bar{M}_{s}M_{s},{\tilde{J}}+1})=(1,0,0,0)$.
The momenta of these waves, $k$ and $k_{s}$, are related to the energy $E$ via

$\displaystyle E=2m_{M}+\frac{k^{2}}{2\mu_{M}}\quad,\quad
E=2m_{M_{s}}+\frac{k_{s}^{2}}{2\mu_{M_{s}}}.$ (43)

The emergent wave is described by the four radial wave functions
$\chi_{\bar{M}_{(s)}M_{(s)},L_{\textrm{out}}\rightarrow{\tilde{J}}}(r)$ with
$\bar{M}_{(s)}M_{(s)}\in\\{\bar{M}M,\bar{M}_{s}M_{s}\\}$ and
$L_{\textrm{out}}\in\\{{\tilde{J}}-1,{\tilde{J}}+1\\}$, where the subscript
indicates the meson content and the coupling of orbital angular momentum
$L_{\textrm{out}}$ and light quark spin $S_{q}=1$ to total angular momentum
${\tilde{J}}$. The boundary conditions for the emergent wave can be formulated
as follows:

  * •

For $r\rightarrow 0$:

$\displaystyle\chi_{\bar{M}_{(s)}M_{(s)},L_{\textrm{out}}\rightarrow{\tilde{J}}}\propto
r^{L_{\textrm{out}}+1}.$ (44)

  * •

For $r\rightarrow\infty$:

$\displaystyle\chi_{\bar{M}M,L_{\textrm{out}}\rightarrow{\tilde{J}}}=it_{\bar{M}_{(s)}M_{(s)},L_{\textrm{in}};\bar{M}M,L_{\textrm{out}}}rh^{(1)}_{L_{\textrm{out}}}(kr)$
(45)
$\displaystyle\chi_{\bar{M}_{s}M_{s},L_{\textrm{out}}\rightarrow{\tilde{J}}}=it_{\bar{M}_{(s)}M_{(s)},L_{\textrm{in}};\bar{M}_{s}M_{s},L_{\textrm{out}}}rh^{(1)}_{L_{\textrm{out}}}(k_{s}r),$
(46)

where

    * –

$\bar{M}_{(s)}M_{(s)}\equiv\bar{M}M$ for $\vec{\alpha}=(1,0,0,0)$ and
$\vec{\alpha}=(0,1,0,0)$, i.e. an incident $\bar{M}M$ wave with
$L_{\textrm{in}}={\tilde{J}}-1$ and $L_{\textrm{in}}={\tilde{J}}+1$,
respectively,

    * –

$\bar{M}_{(s)}M_{(s)}\equiv\bar{M}_{s}M_{s}$ for $\vec{\alpha}=(0,0,1,0)$ and
$\vec{\alpha}=(0,0,0,1)$, i.e. an incident $\bar{M}_{s}M_{s}$ wave with
$L_{\textrm{in}}={\tilde{J}}-1$ and $L_{\textrm{in}}={\tilde{J}}+1$,
respectively.

These boundary conditions define 16 quantities,
$t_{\bar{M}_{(s)}M_{(s)},L_{\textrm{in}};\bar{M}_{(s)}M_{(s)},L_{\textrm{out}}}$,
which represent scattering amplitudes and can be combined to the $4\times 4$ T
matrix

$\displaystyle\mbox{T}_{\tilde{J}}=\left(\begin{array}[]{cccc}t_{\bar{M}M,{\tilde{J}}-1;\bar{M}M,{\tilde{J}}-1}&t_{\bar{M}M,{\tilde{J}}+1;\bar{M}M,{\tilde{J}}-1}&t_{\bar{M}_{s}M_{s},{\tilde{J}}-1;\bar{M}M,{\tilde{J}}-1}&t_{\bar{M}_{s}M_{s},{\tilde{J}}+1;\bar{M}M,{\tilde{J}}-1}\\\
t_{\bar{M}M,{\tilde{J}}-1;\bar{M}M,{\tilde{J}}+1}&t_{\bar{M}M,{\tilde{J}}+1;\bar{M}M,{\tilde{J}}+1}&t_{\bar{M}_{s}M_{s},{\tilde{J}}-1;\bar{M}M,{\tilde{J}}+1}&t_{\bar{M}_{s}M_{s},{\tilde{J}}+1;\bar{M}M,{\tilde{J}}+1}\\\
t_{\bar{M}M,{\tilde{J}}-1;\bar{M}_{s}M_{s},{\tilde{J}}-1}&t_{\bar{M}M,{\tilde{J}}+1;\bar{M}_{s}M_{s},{\tilde{J}}-1}&t_{\bar{M}_{s}M_{s},{\tilde{J}}-1;\bar{M}_{s}M_{s},{\tilde{J}}-1}&t_{\bar{M}_{s}M_{s},{\tilde{J}}+1;\bar{M}_{s}M_{s},{\tilde{J}}-1}\\\
t_{\bar{M}M,{\tilde{J}}-1;\bar{M}_{s}M_{s},{\tilde{J}}+1}&t_{\bar{M}M,{\tilde{J}}+1;\bar{M}_{s}M_{s},{\tilde{J}}+1}&t_{\bar{M}_{s}M_{s},{\tilde{J}}-1;\bar{M}_{s}M_{s},{\tilde{J}}+1}&t_{\bar{M}_{s}M_{s},{\tilde{J}}+1;\bar{M}_{s}M_{s},{\tilde{J}}+1}\\\
\end{array}\right).$ (51)

The T matrix is related to the S matrix in the usual way,

$\displaystyle\mbox{S}_{\tilde{J}}=1+2i\mbox{T}_{\tilde{J}}.$ (52)

For ${\tilde{J}}=0$ one has to discard the contributions to the incident wave
with $L_{\textrm{in}}={\tilde{J}}-1$ and to the emergent wave with
$L_{\textrm{out}}={\tilde{J}}-1$. The Schrödinger equation (34) is then
reduced from five to three channels and $\mbox{T}_{\tilde{J}}$ from a $4\times
4$ to a $2\times 2$ matrix. This ${\tilde{J}}=0$ equation is extensively
discussed in Ref. [29].

In analogy to a single channel scattering problem, where the phase shift
$\delta$ is defined via $1+2i\mbox{T}=\mbox{S}=\exp(2i\delta)$, one can define
the eigenphase sum [50, 51, 52] for multi-channel scattering via

$\displaystyle\textrm{det}(\mbox{S}_{\tilde{J}})=\exp(2i\delta_{\tilde{J}}).$
(53)

The eigenphase sum $\delta_{\tilde{J}}$ is identical to a sum of phase shifts,
where each phase shift corresponds to one of the eigenvalues of the S matrix.

##  III Numerical results

###  III.1 Numerical methods to solve the coupled channel Schrödinger equation
and to determine the poles of the T matrix

In Section II.3 we defined the entries of the T matrix (51) as the a priori
unknown coefficients
$t_{\bar{M}_{(s)}M_{(s)},L_{\textrm{in}};\bar{M}_{(s)}M_{(s)},L_{\textrm{out}}}$
appearing in the $r\rightarrow\infty$ boundary conditions (45) and (46). To
determine these coefficients, one has to solve the coupled channel Schrödinger
equation (34). To cross check our results, we used two rather different
numerical methods. The first method corresponds to discretizing the radial
coordinate by a uniform grid and solving the resulting system of linear
equations using methods from standard textbooks (for details see Ref. [28]).
The second method corresponds to using an ordinary 4th order Runge-Kutta
algorithm.

To find the poles of $\mbox{T}_{\tilde{J}}$ in the complex energy plane,
chracterized by at least one of its eigenvalues approaching infinity, we
applied the Newton-Raphson method to find the roots of
$1/\textrm{det}(\mbox{T}_{\tilde{J}})$.

###  III.2 Input parameters and error analysis

In the following we present results for heavy $b$ quarks, i.e. $Q\equiv b$. We
use $m_{Q}=4.977\,\text{GeV}$ from quark models [53] and the spin-averaged
mass of the $B$ meson and the $B^{*}$ meson, i.e.
$m_{M}=(m_{B}+3m_{B^{*}})/4=5.313\,\text{GeV}$, as well as of the $B_{s}$
meson and the $B_{s}^{*}$ meson, i.e.
$m_{M_{s}}=(m_{B_{s}}+3m_{B_{s}^{*}})/4=5.403\,\text{GeV}$. The lattice data
from Ref. [30] we are using to determine the parameters of the potential
parameterizations (see Section II.2) was generated with a light quark mass
slightly below the physical strange quark mass. This is reflected by
$E_{\text{threshold}}=10.790\,\text{GeV}$, which is much closer to the spin-
averaged $B_{s}^{(*)}B_{s}^{(*)}$ threshold than to the spin-averaged
$B^{(*)}B^{(*)}$ threshold.

The uncertainties of the lattice data provided in Ref. [30] are propagated via
resampling. We generated 1000 statistically independent samples and repeated
our computations on each of them. For our results we quote asymmetric errors,
defined by the 16th and 84th percentile.

###  III.3 Eigenphase sums and poles of the T matrix

We have computed the eigenphase sum $\delta_{\tilde{J}}$ defined in Eq. (53)
for ${\tilde{J}}=0,1,2,3$ as function of the energy above the spin-averaged
$B^{(*)}B^{(*)}$ threshold at $10.627\,\text{GeV}$. We show the respective
plots in Fig. 2. For rather stable resonances, which are clearly separated in
their energies, such plots show pronounced steps of order $\pi$. The locations
of these steps then correspond to the resonance masses and the slopes are
inversely proportional to the associated decay widths. In our case, however,
where some resonances have large decay widths and their energy levels are
close, it is hard to identify them in clear and unique way.

Figure 2: Eigenphase sum $\delta_{\tilde{J}}$ as function of the energy $E$.
Narrow resonances are indicated by pronounced steps. We also show the real
parts of the positions of the poles of $\mbox{T}_{\tilde{J}}$ using green
dotted lines (the line at $\approx 10.798\,\text{GeV}$ in the lower left plot
does not correspond to the real part of a pole, but to a fit to the eigenphase
sum; see the technical discussion in Section III.3.1). A pale red dotted line
indicates the spin averaged $B_{s}^{(*)}B_{s}^{(*)}$ threshold, the light blue
solid line $E_{\text{threshold}}$. Results for energies inside the light-blue
shaded region, i.e. above $\approx 11.025\,\text{GeV}$, should not be trusted,
since we neglect decay channels containing a negative and a positive parity
heavy-light meson.

To define and to compute masses and decay widths of bottomonium resonances in
a clearer and more definite way, we have analytically continued our scattering
problem to the complex energy plane. Then we have determined the poles of the
T matrix (51) in the complex energy plane numerically. The positions of these
poles can be related to masses and decay widths according to

$\displaystyle
m=\textrm{Re}(E_{\text{pole}})\quad,\quad\Gamma=-2\textrm{Im}(E_{\text{pole}}),$
(54)

where $E_{\text{pole}}$ denotes the complex pole energy.

In Fig. 3 we show all poles of the T matrix (51) for ${\tilde{J}}=0,1,2,3$ up
to $11.2\,\text{GeV}$. Colored point clouds represent 1000 independent
computations with resampled lattice data from Ref. [30]. These point clouds
are used to determine statistical errors, which is straightforward, because
there are clear gaps between the point clouds. Bound states correspond to
poles located on the real axis below the $B^{(*)}B^{(*)}$ threshold, while
resonances correspond to poles above this threshold with a non-vanishing
negative imaginary part. The pole positions and their statistical errors are
indicated by the black crosses. The $B^{(*)}B^{(*)}$ threshold at
$10.627\,\text{GeV}$ and the $B_{s}^{(*)}B_{s}^{(*)}$ threshold at
$10.807\,\text{GeV}$ are indicated by pale red dotted lines. Our results can
only be trusted up to the threshold of one negative parity heavy-light meson
and one positive parity heavy-light meson at around $11.025\,\text{GeV}$,
because the corresponding decay channel is not included in our Schrödinger
equation (4). The region above this threshold is shaded in light-blue.

We summarize these theoretical predictions of bottomonium masses and decay
widths in Table 2, together with available experimental results. Additionally,
we show the same data in a graphical way in Fig. 4.

Figure 3: Positions of the poles of the T matrix (51) in the complex energy
plane for ${\tilde{J}}=0$ (top left), ${\tilde{J}}=1$ (top right),
${\tilde{J}}=2$ (bottom left) and ${\tilde{J}}=3$ (bottom right) representing
bound states and resonances up to $11.2\,\text{GeV}$. Colored point clouds
represent 1000 independent computations with resampled lattice data, while
black dots correspond to the mean values and error bars. The pale red dotted
lines indicate the spin averaged $B^{(*)}B^{(*)}$ and $B_{s}^{(*)}B_{s}^{(*)}$
thresholds, the light blue solid line $E_{\text{threshold}}$. Results for
energies inside the light-blue shaded region, i.e. above $\approx
11.025\,\text{GeV}$, should not be trusted, since we neglect decay channels
containing a negative and a positive parity heavy-light meson. Table 2: Masses
and decay widths of $I=0$ bottomonium with
${\tilde{J}}^{PC}=0^{++},1^{--},2^{++},3^{--}$ obtained from the poles of the
T matrix (51), where errors are purely statistical. An exception is the $n=3$
resonance for ${\tilde{J}}=2$, which was extracted by a fit to the eigenphase
sum $\delta_{2}$ (see the technical discussion in Section III.3.1). The spin-
averaged $B^{(*)}B^{(*)}$ and $B_{s}^{(*)}B_{s}^{(*)}$ thresholds are
indicated by dashed lines. Results for energies above the threshold of one
negative and another positive heavy-light meson at $\approx
11.025\,\text{GeV}$ are marked by a light-blue background and should not be
trusted, since we neglect the corresponding decay channels. (Results for
${\tilde{J}}^{PC}=0^{++}$ were already presented in Ref. [29].)  Figure 4:
Graphical summary of theoretical predictions and experimental results for
masses of $I=0$ bottomonium with
${\tilde{J}}^{PC}=0^{++},1^{--},2^{++},3^{--}$. We also show the quarkonium
and meson-meson composition as discussed in Section III.4: $\%\bar{Q}Q$ in
orange, $(\%\bar{M}M)_{{\tilde{J}}-1}+(\%\bar{M}M)_{{\tilde{J}}+1}$ in blue
and
$(\%\bar{M}_{s}M_{s})_{{\tilde{J}}-1}+(\%\bar{M}_{s}M_{s})_{{\tilde{J}}+1}$ in
green.

####  III.3.1 Technical aspects of pole finding

In addition to the physical poles of the T matrix, which correspond to
bottomonium bound states and resonances, there are also unphysical poles,
which are caused by numerical inaccuracies. Such unphysical poles can be
identified, by varying numerical parameters (e.g. the spacing of the grid, the
Runge-Kutta step size or the large, but finite value of $r$ replacing
$r\rightarrow\infty$ in the boundary conditions (45) and (46)). If a pole is
unstable with respect to these parameters, it is clearly an unphysical pole.

In cases, where a physical pole and unphysical poles are close, pole finding
might become a delicate task. To verify that we did not miss any of the
physical poles with our pole finding algorithm, we have also determined the
poles after removing the $\bar{B}^{(*)}_{s}B^{(*)}_{s}$ decay channels. In
general, the pole positions in the 3-flavor case are similar to those in the
2-flavor case and consistent with the steps observed in the eigenphase sums
shown in Fig. 2. An exception is the $n=3$ resonance in the ${\tilde{J}}=2$
sector, which or pole finding algorithm was unable to identify in the 3-flavor
case. It seems to be masked by unphysical poles. Thus, to determine the
corresponding resonance parameters, we have used the eigenphase sum
$\delta_{2}$. We have performed a 4-parameter fit with
$\alpha+\beta\text{arctan}((2/\Gamma)(E-m))$ to the data points for
$\delta_{2}$ in the region of the step at $\approx 10.798\,\text{GeV}$ (see
Fig. 5). In this way we have obtained the resonance mass
$m=10.799_{-2}^{+2}\,\text{GeV}$ and the decay width
$\Gamma=13.0_{-2.0}^{+1.8}\,\text{MeV}$, where error bars are determined by
resampling lattice data from Ref. [30] and by varying the fit range. As a
cross check we have also determined the parameters of the $n=4$ resonance in
the ${\tilde{J}}=2$ sector in the same way and find
$m=11.038_{-3}^{+5}\,\text{GeV}$ and $\Gamma=42.29_{-1.2}^{+2.0}\,\text{MeV}$.
These results are fully consistent with those obtained from our pole search
(cf. Table 2).

Figure 5: Fits of $\alpha+\beta\text{arctan}((2/\Gamma)(E-m))$ to the data
points for the eigenphase sum $\delta_{2}$ in the region of the steps at
$\approx 10.798\,\text{GeV}$ and $\approx 11.038\,\text{GeV}$.

####  III.3.2 Comparison to experimental results

For the predicted low-lying states it is straightforward to assign
experimental counterparts.

  * •

The ${\tilde{J}}=0$ states with $n=1,2,3,4$ correspond to
$\eta_{b}(1S)\equiv\Upsilon(1S)$, to $\Upsilon(2S)$, to $\Upsilon(3S)$ and to
$\Upsilon(4S)$.

  * •

The ${\tilde{J}}=1$ states with $n=1,2,3$ correspond to
$h_{b}(1P)\equiv\chi_{b0}(1P)\equiv\chi_{b1}(1P)\equiv\chi_{b2}(1P)$, to
$h_{b}(2P)\equiv\chi_{b0}(2P)\equiv\chi_{b1}(2P)\equiv\chi_{b2}(2P)$ and to
$\chi_{b1}(3P)$.

  * •

The ${\tilde{J}}=2$ state with $n=1$ corresponds to $\Upsilon(1D)$.

Our masses exhibit a pattern, which is quite similar to that found in
experiments. The largest discrepancies are observed for the lowest states,
most prominently for the ground state $\eta_{b}(1S)\equiv\Upsilon(1S)$. This,
however, does not indicate particular problems with these states, but is
rather a consequence of choosing $E_{\text{threshold}}\approx 2m_{M_{s}}$ as
reference point to introduce the energy scale. Compared to this reference
point the errors of our predicted states are of order $10\%$. To a large part
these errors can be compensated by a global multiplicative factor. A possible
reason for this error might be the scale setting in the lattice QCD
computation of Ref. [30], which is based on defining $r_{0}=0.5\,\text{fm}$,
while more recent lattice investigations indicate a smaller value for $r_{0}$
[54]. Using e.g. $r_{0}=0.45\,\text{fm}$ instead of $r_{0}=0.5\,\text{fm}$
reduces the error for $\eta_{b}(1S)\equiv\Upsilon(1S)$ by around $50\%$.
Moreover, our potentials were determined from lattice QCD data obtained at a
single coarse lattice spacing and with only two light quark flavors close to
the mass of the physical strange quark. Thus, for precise quantitative
predictions a more accurate lattice QCD computation of the relevant potentials
and mixing angles will be necessary (see also the outlook in Section IV).

The ${\tilde{J}}=0$ resonance with $n=5$ has a mass close to the experimental
result for $\Upsilon(10753)$, which was recently reported by Belle [55]. In a
previous publication we investigated the structure of this state within the
same setup and found that it is meson-meson dominated with just a small quark-
antiquark component [29] (see also Section III.4, in particular Table 3).
Thus, since it is not an ordinary quarkonium state and the heavy quark spin
can be $1^{--}$, it can be classified as a $Y$ type crypto-exotic state. Note
that we found another resonance in that energy region with quantum numbers
${\tilde{J}}=2$ and $n=3$. This state is, however, farther away from the
experimental result for $\Upsilon(10753)$ ($\approx 46\,\text{MeV}$ difference
for ${\tilde{J}}=2$, $n=3$ compared to $\approx 20\,\text{MeV}$ difference for
${\tilde{J}}=0$, $n=5$) and, thus, an identification with $\Upsilon(10753)$
seems less likely. On the other hand, since our results exhibit certain
systematic errors, as discussed in the previous paragraph, we are not in a
position to fully exclude such an identification.

The resonances $\Upsilon(10860)$ and $\Upsilon(11020)$ are typically
interpreted as $\Upsilon(5S)$ and $\Upsilon(6S)$. However, from the
experimental perspective they could as well correspond to $D$ wave states. The
${\tilde{J}}=0$ $S$ wave resonance with $n=6$ is rather close to the mass of
$\Upsilon(10860)$, whereas there is no ${\tilde{J}}=2$ $D$ wave resonance in
that energy region. Thus, our results support the interpretation of
$\Upsilon(10860)$ as $\Upsilon(5S)$. Concerning $\Upsilon(11020)$, the
${\tilde{J}}=0$ resonance with $n=7$ and the ${\tilde{J}}=2$ resonance with
$n=4$ have almost the same mass and are both close to the mass of
$\Upsilon(11020)$. Moreover, that mass is already close to the threshold of a
negative parity $B$ or $B^{\ast}$ and a positive parity $B_{0}^{\ast}$ or
$B_{1}^{\ast}$ meson, a channel we have not yet included in our approach.
Thus, we cannot decide, whether the $\Upsilon(11020)$ is indeed an $S$ wave or
rather a $D$ wave state.

###  III.4 Quarkonium and meson-meson composition

Using techniques developed in Ref. [29] we also study the structure and quark
content of bound states and resonances, to clarify, whether they are
conventional quarkonia or there are sizable $\bar{Q}Q\bar{q}q$ four-quark
components. To this end, we compute for each state the percentages of
quarkonium with $L={\tilde{J}}$ and of $\bar{M}M$ and $\bar{M}_{s}M_{s}$
meson-meson pairs with $L_{\text{out}}={\tilde{J}}-1,{\tilde{J}}+1$,

$\displaystyle\%\bar{Q}Q$
$\displaystyle=\frac{Q}{Q+M_{{\tilde{J}}-1}+M_{{\tilde{J}}+1}+M_{s,{\tilde{J}}-1}+M_{s,{\tilde{J}}+1}}$
(55) $\displaystyle(\%\bar{M}M)_{L_{\text{out}}}$
$\displaystyle=\frac{M_{L_{\text{out}}}}{Q+M_{{\tilde{J}}-1}+M_{{\tilde{J}}+1}+M_{s,{\tilde{J}}-1}+M_{s,{\tilde{J}}+1}}$
(56) $\displaystyle(\%\bar{M}_{s}M_{s})_{L_{\text{out}}}$
$\displaystyle=\frac{M_{s,L_{\text{out}}}}{Q+M_{{\tilde{J}}-1}+M_{{\tilde{J}}+1}+M_{s,{\tilde{J}}-1}+M_{s,{\tilde{J}}+1}},$
(57)

where

$\displaystyle Q$
$\displaystyle=\int_{0}^{R_{\text{max}}}\text{d}r\,\Big{|}u_{{\tilde{J}}}(r)\Big{|}^{2}$
(58) $\displaystyle M_{L_{\text{out}}}$
$\displaystyle=\int_{0}^{R_{\text{max}}}\text{d}r\,\Big{|}\chi_{\bar{M}M,L_{\text{out}}\rightarrow{\tilde{J}}}(r)\Big{|}^{2}$
(59) $\displaystyle M_{s,L_{\text{out}}}$
$\displaystyle=\int_{0}^{R_{\text{max}}}\text{d}r\,\Big{|}\chi_{\bar{M}_{s}M_{s},L_{\text{out}}\rightarrow{\tilde{J}}}(r)\Big{|}^{2}.$
(60)

$u_{{\tilde{J}}}(r)$ ,
$\chi_{\bar{M}M,L_{\text{out}}\rightarrow{\tilde{J}}}(r)$ and
$\chi_{\bar{M}_{s}M_{s},L_{\text{out}}\rightarrow{\tilde{J}}}(r)$ are the
solutions of the radial Schrödinger equation (34) for real energy
$\text{Re}(E_{\text{pole}})$, where $E_{\text{pole}}$ is the position of the
corresponding pole of the T matrix in the complex energy plane.

In case of a bound state, $Q$, $M_{L_{\text{out}}}$ and $M_{s,L_{\text{out}}}$
approach constants for
$R_{\text{max}}\raisebox{-2.15277pt}{$\,\stackrel{{\scriptstyle>}}{{\scriptstyle\sim}}\,$}2.0\,\text{fm}$
as indicated by Figure 6 to Figure 9. The corresponding asymptotic values of
$\%\bar{Q}Q$, $(\%\bar{M}M)_{L_{\text{out}}}$ and
$(\%\bar{M}_{s}M_{s})_{L_{\text{out}}}$ for bound states are collected in
Table 3.

Table 3: Percentages of quarkonium and of meson-meson pairs for $I=0$
bottomonium with ${\tilde{J}}^{PC}=0^{++},1^{--},2^{++},3^{--}$. In addition
to statistical errors we also provide systematic uncertainties, which are
discussed in Section III.4. The spin-averaged $B^{(*)}B^{(*)}$ and
$B_{s}^{(*)}B_{s}^{(*)}$ thresholds are indicated by dashed lines. Results for
energies above the threshold of one negative and another positive heavy-light
meson at $\approx 11.025\,\text{GeV}$ are marked by a light-blue background
and should not be trusted, since we neglect the corresponding decay channels.
(Results for ${\tilde{J}}^{PC}=0^{++}$ were already presented in Ref. [29].)

In case of a resonance, $M_{L_{\text{out}}}$ and $M_{s,L_{\text{out}}}$ are
linearly rising functions for large $R_{\text{max}}$, because in that region
$\chi_{\bar{M}M,L_{\text{out}}\rightarrow{\tilde{J}}}(r)$ and
$\chi_{\bar{M}_{s}M_{s},L_{\text{out}}\rightarrow{\tilde{J}}}(r)$ represent
emergent spherical waves. The corresponding slopes, however, are rather small
and result in changes of only a few percent in $\%\bar{Q}Q$,
$(\%\bar{M}M)_{L_{\text{out}}}$ and $(\%\bar{M}_{s}M_{s})_{L_{\text{out}}}$ in
the interval $1.8\,\text{fm}<R_{\text{max}}<3.0\,\text{fm}$. At such
separations $R_{\text{max}}$ the quarkonium component $u_{{\tilde{J}}}(r)$ is
already negligible and
$\chi_{\bar{M}M,L_{\text{out}}\rightarrow{\tilde{J}}}(r)$ and
$\chi_{\bar{M}_{s}M_{s},L_{\text{out}}\rightarrow{\tilde{J}}}(r)$ are almost
pure emergent spherical waves. Thus the corresponding resonance is contained
inside a sphere of radius $R_{\text{max}}$. We evaluate $\%\bar{Q}Q$,
$(\%\bar{M}M)_{L_{\text{out}}}$ and $(\%\bar{M}_{s}M_{s})_{L_{\text{out}}}$ at
the center of this interval, i.e. at $R_{\text{max}}=2.4\,\text{fm}$, but
assign asymmetric systematic uncertainties
$|\%\bar{Q}Q(R_{\text{max}}=1.8\,\text{fm})-\%\bar{Q}Q(R_{\text{max}}=2.4\,\text{fm})|$
and
$|\%\bar{Q}Q(R_{\text{max}}=3.0\,\text{fm})-\%\bar{Q}Q(R_{\text{max}}=2.4\,\text{fm})|$
to $\%\bar{Q}Q$ and in an analogous way also to
$(\%\bar{M}M)_{L_{\text{out}}}$ and to
$(\%\bar{M}_{s}M_{s})_{L_{\text{out}}}$. These results are collected in Table
3. Moreover, $\%\bar{Q}Q$, $(\%\bar{M}M)_{L_{\text{out}}}$ and
$(\%\bar{M}_{s}M_{s})_{L_{\text{out}}}$ as functions of $R_{\text{max}}$ are
shown in Figure 6 to Figure 9.

The majority of bound states have masses significantly below the
$B^{(*)}B^{(*)}$ threshold with binding energies of the order of
$100\,\text{MeV}$ or larger (${\tilde{J}}=0$, $n=1,2,3$; ${\tilde{J}}=1$,
$n=1,2$; ${\tilde{J}}=2$, $n=1,2$; ${\tilde{J}}=3$, $n=1$). These states
consist mostly of quarkonium, $\%\bar{Q}Q\approx 74\%\ldots 84\%$. The largest
quarkonium components are present for the lowest states with orbital angular
momentum $L=0$, i.e. for ${\tilde{J}}=0$ and $n=1,2$. Because of the non-
vanishing mixing angle $\theta\approx 0.35\ldots 0.40$ for
$r\raisebox{-2.15277pt}{$\,\stackrel{{\scriptstyle<}}{{\scriptstyle\sim}}\,$}1.0\,\text{fm}$
obtained by a full lattice QCD computation (see Ref. [30] and our previous
article [28], in particular Section III), larger quarkonium percentages
$\%\bar{Q}Q$ are excluded. The reason is that for two light quark flavors the
ground state potential in the $\Sigma_{g}^{+}$ sector, which is of central
importance for bound states, is a linear superposition of
$1-\sin^{2}(\theta)\approx 85\%\ldots 88\%$ quarkonium and of
$\theta^{2}\approx 12\%\ldots 15\%$ meson-meson. Thus, a very similar
composition is generated dynamically for the deeply bound states by our
Schrödinger equation (4), which was derived to be consistent with the lattice
QCD potentials from Ref. [30] and the corresponding mixing angle.

However, there is one bound state (${\tilde{J}}=1$, $n=3$) with $m\approx
10.594\,\text{GeV}$, rather close to the $B^{(*)}B^{(*)}$ threshold at
$10.627\,\text{GeV}$, where the quarkonium component is already significantly
reduced, $\%\bar{Q}Q\approx 58\%$. A reason for that could be that the
components of the radial wave functions are significantly farther extended, up
to $r\approx 1.5\,\text{fm}\ldots 2.0\,\text{fm}$ (see Figure 7). Close to the
string breaking distance $r_{\text{sb}}$ the mixing angle changes rapidly from
$\theta\approx 0.35\ldots 0.40$ to $\theta\approx\pi/2$. Thus, for
$r\raisebox{-2.15277pt}{$\,\stackrel{{\scriptstyle>}}{{\scriptstyle\sim}}\,$}r_{\text{sb}}\approx
1.25\,\text{fm}$ the ground state potential in the $\Sigma_{g}^{+}$ sector
corresponds almost exclusively to a $B^{(*)}B^{(*)}$ pair. Consequently,
$(\%\bar{M}M)_{L_{\text{out}}}$ is significantly enhanced for such a spatially
extended state compared to the more tightly bound states discussed in the
previous paragraph.

The resonances with ${\tilde{J}}=0$, $n=4$ and ${\tilde{J}}=3$, $n=2$ are
slightly above the $B^{(*)}B^{(*)}$ threshold, i.e. energetically they are
extremely close to bound states. In both cases this is reflected by a roughly
equal mix of the quarkonium component and the $B^{(*)}B^{(*)}$ component(s).
$B_{s}^{(*)}B_{s}^{(*)}$ contributions, on the other hand, are almost
negligible, because the corresponding threshold is more than $160\,\text{MeV}$
above.

For higher resonances the meson-meson components start to dominate and the
quarkonium component is somewhere between $8\%\ldots 35\%$. The corresponding
widths tend to be sizable and the resonances are rather unstable. Thus, it is
not surprising that these resonances are mostly meson-meson states.

For resonances above $\approx 11.025\,\text{GeV}$, which is the threshold of
one negative and another positive heavy-light meson, the quarkonium components
start to increase again. We interpret this, however, rather as a consequence
of our neglect of decay channels in that energy region than as a solid and
meaningful physics result. We note again that also energy levels above
$\approx 11.025\,\text{GeV}$ should not be trusted or at least be taken with
extreme caution, as e.g. discussed already in Section III.3.

Figure 6: Percentages of quarkonium and of meson-meson pairs for $I=0$
bottomonium with ${\tilde{J}}^{PC}=0^{++}$ as functions of $R_{\textrm{max}}$.

Figure 7: Percentages of quarkonium and of meson-meson pairs for $I=0$
bottomonium with ${\tilde{J}}^{PC}=1^{--}$ as functions of $R_{\textrm{max}}$.

Figure 8: Percentages of quarkonium and of meson-meson pairs for $I=0$
bottomonium with ${\tilde{J}}^{PC}=2^{++}$ as functions of $R_{\textrm{max}}$.

Figure 9: Percentages of quarkonium and of meson-meson pairs for $I=0$
bottomonium with ${\tilde{J}}^{PC}=3^{--}$ as functions of $R_{\textrm{max}}$.

##  IV Conclusions and outlook

We extended our previous work for ${\tilde{J}}=0$ [28, 29] and derived a
coupled channel Schrödinger equation for arbitrary ${\tilde{J}}$ to
investigate bottomonium bound states and resonances with $I=0$ using lattice
QCD static potentials. For each ${\tilde{J}}>0$ we take into account five
coupled channels and, thus, have developed an approach to study complicated
resonances based on lattice QCD. We solved the coupled channel Schrödinger
equation for ${\tilde{J}}=0,1,2,3$ and found multiple bound states and
resonances for all values of ${\tilde{J}}$. We also explored the structure of
these states by computing their quarkonium and meson-meson percentages.

Our results for masses of bound states and resonances are consistent with
experimentally observed states within expected errors, which arise from
resorting to several approximations. We find several bound states in the
sectors ${\tilde{J}}=0,1,2$, which all have a clear experimental counterpart.

The resonance $\Upsilon(10753)$, which was recently found by Belle [55] also
appears in our $S$ wave spectrum as ${\tilde{J}}=0$, $n=5$. See our last
publication [29], where we discuss in detail that this is a meson-meson
dominated state ,which can be classified as an $\Upsilon$ type crypto-exotic
state. We now find another state with a comparable mass and decay width in our
$D$ wave spectrum as ${\tilde{J}}=2$, $n=3$, i.e. there are two distinct
resonances very close in the spectrum. This is not surprising and can be
observed already in constituent quark models, since $D$ wave quarkonium states
have energies comparable to radially excited $S$ wave quarkonium states [56].

We can confirm the interpretation of $\Upsilon(10860)$ as $\Upsilon(5S)$.
Concerning $\Upsilon(11020)$ we find indications of corresponding resonances
in both the $S$ wave sector (${\tilde{J}}=0$, $n=7$) as well as the $D$ wave
sector (${\tilde{J}}=2$, $n=4$). We interpret these two resonances as two
states as briefly discussed in the previous paragraph. We note that it would
be very interesting, if experiments studying bottomonium could disentangle
these two states.

In what concerns a possibly existing bottomonium state close to the
$B^{(*)}B^{(*)}$ threshold, as counterpart to the $X(3872)$ charmonium state,
we do not find one. However, we find a ${\tilde{J}}=2$, $n=3$ state very close
to the $B_{s}^{(*)}B_{s}^{(*)}$ threshold, which could have similarities to
$X(3872)$ [48, 49] and its composite nature proposed by lattice QCD
computations [57] and by data analysis [58]. This state has a meson-meson
component significantly larger than its quarkonium component, $79\%$ versus
$21\%$.

Our aim for the future is to reduce systematic errors as much as possible, to
be able to reproduce the experimentally observed states not only on a
qualitative level, but rather precisely with combined statistical and
systematic errors of only a few MeV. To reduce the systematic errors, we plan
to compute the necessary potentials $V_{\bar{Q}Q}(r)$, $V_{\textrm{mix}}(r)$,
$V_{\bar{M}M,\parallel}(r)$ and $V_{\bar{M}M,\perp}(r)$ via lattice QCD
explicitly, avoiding to perform algebraic operations on string breaking
potentials from Ref. [30]. This should eliminate several sources of
uncertainty and additionally create opportunities for new and interesting
results:

  * •

With up-to-date lattice results we might also obtain a more accurate scale
setting, i.e. a more precise value for the lattice spacing $a$ or,
equivalently, the Sommer parameter $r_{0}$. Moreover, we should be in a
position to calibrate the bottom quark mass $m_{b}$ rather precisely and,
thus, avoid using a value from quark models.

  * •

There is a discrepancy between our lightest bound states and the corresponding
experimental results. We expect that this can be resolved by a precise
computation of the mixing angle $\theta(r)$, which has a strong impact on the
potential matrix (see Section II A of Ref. [28] for details).

  * •

Related to the previous item is a precise determination of the mixing
potentials between quarkonium and meson-meson pairs, providing quantitative
first principles results for quantities estimated and used in a large number
of quark models for many years [59, 60, 61, 62, 63, 64].

  * •

The details of the short range and the long range parts of the entries of the
potential matrix, if computed with sufficient precision, should be of interest
for models based on hadron-hadron interactions.

Another significant source of systematic error is the neglect of heavy spin
effects. To achieve the desired level of precision, we plan to include them
using methods developed in a related project [65] and possibly carry out a
lattice QCD computation of $1/m_{b}$ and $1/m_{b}^{2}$ corrections at least
for the confining potential $V_{\bar{Q}Q}(r)$ (see e.g. Refs. [66, 67, 68, 69,
70, 71]).

Finally, it might be interesting to include decay channels to a negative and a
positive parity heavy-light meson pair, to make solid predictions up to the
threshold of two positive parity mesons at around $11.525\,\text{GeV}$. This
would allow us to obtain information about several states, which are not yet
measurable by experiments.

###### Acknowledgements.

We acknowledge useful discussions with Gunnar Bali, Eric Braaten, Marco
Cardoso, Francesco Knechtli, Vanessa Koch and Sasa Prelovsek. P.B. and N.C.
acknowledge the support of CeFEMA under the FCT contract for R&D Units
UIDB/04540/2020. N.C. acknowledges the FCT contract SFRH/BPD/109443/2015. L.M.
acknowledges support by a Karin and Carlo Giersch Scholarship of the Giersch
foundation. M.W. acknowledges support by the Heisenberg Programme of the
Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) - project
number 399217702. Calculations on GPU servers of CeFEMA partly supported by
NVIDIA were conducted for this research. Calculations on the GOETHE-HLR and on
the FUCHS-CSC high-performance computer of the Frankfurt University were
conducted for this research. We would like to thank HPC-Hessen, funded by the
State Ministry of Higher Education, Research and the Arts, for programming
advice.

## References

  * Aaij  _et al._ [2021a] R. Aaij _et al._ (LHCb), Phys. Rev. Lett. 127, 082001 (2021a), arXiv:2103.01803 [hep-ex] .
  * Ablikim  _et al._ [2021] M. Ablikim _et al._ (BESIII), Phys. Rev. Lett. 126, 102001 (2021), arXiv:2011.07855 [hep-ex] .
  * Aaij  _et al._ [2020a] R. Aaij _et al._ (LHCb), Phys. Rev. Lett. 125, 242001 (2020a), arXiv:2009.00025 [hep-ex] .
  * Aaij  _et al._ [2020b] R. Aaij _et al._ (LHCb), Phys. Rev. D 102, 112003 (2020b), arXiv:2009.00026 [hep-ex] .
  * Rodas  _et al._ [2019] A. Rodas _et al._ (JPAC), Phys. Rev. Lett. 122, 042002 (2019), arXiv:1810.04171 [hep-ph] .
  * Aaij  _et al._ [2018] R. Aaij _et al._ (LHCb), Eur. Phys. J. C 78, 1019 (2018), arXiv:1809.07416 [hep-ex] .
  * Aaij  _et al._ [2015] R. Aaij _et al._ (LHCb), Phys. Rev. D 92, 112009 (2015), arXiv:1510.01951 [hep-ex] .
  * Ablikim  _et al._ [2015] M. Ablikim _et al._ (BESIII), Phys. Rev. D 92, 092006 (2015), arXiv:1509.01398 [hep-ex] .
  * Chilikin  _et al._ [2014] K. Chilikin _et al._ (Belle), Phys. Rev. D 90, 112009 (2014), arXiv:1408.6457 [hep-ex] .
  * Aaij  _et al._ [2014] R. Aaij _et al._ (LHCb), Phys. Rev. Lett. 112, 222002 (2014), arXiv:1404.1903 [hep-ex] .
  * Ablikim  _et al._ [2014a] M. Ablikim _et al._ (BESIII), Phys. Rev. Lett. 112, 022001 (2014a), arXiv:1310.1163 [hep-ex] .
  * Ablikim  _et al._ [2013a] M. Ablikim _et al._ (BESIII), Phys. Rev. Lett. 111, 242001 (2013a), arXiv:1309.1896 [hep-ex] .
  * Ablikim  _et al._ [2014b] M. Ablikim _et al._ (BESIII), Phys. Rev. Lett. 112, 132001 (2014b), arXiv:1308.2760 [hep-ex] .
  * Chilikin  _et al._ [2013] K. Chilikin _et al._ (Belle), Phys. Rev. D 88, 074026 (2013), arXiv:1306.4894 [hep-ex] .
  * Xiao  _et al._ [2013] T. Xiao, S. Dobbs, A. Tomaradze, and K. K. Seth, Phys. Lett. B 727, 366 (2013), arXiv:1304.3036 [hep-ex] .
  * Liu  _et al._ [2013] Z. Q. Liu _et al._ (Belle), Phys. Rev. Lett. 110, 252002 (2013), [Erratum: Phys.Rev.Lett. 111, 019901 (2013)], arXiv:1304.0121 [hep-ex] .
  * Ablikim  _et al._ [2013b] M. Ablikim _et al._ (BESIII), Phys. Rev. Lett. 110, 252001 (2013b), arXiv:1303.5949 [hep-ex] .
  * Ketzer [2012] B. Ketzer, PoS QNP2012, 025 (2012), arXiv:1208.5125 [hep-ex] .
  * Bondar  _et al._ [2012] A. Bondar _et al._ (Belle), Phys. Rev. Lett. 108, 122001 (2012), arXiv:1110.2251 [hep-ex] .
  * Mizuk  _et al._ [2009] R. Mizuk _et al._ (Belle), Phys. Rev. D 80, 031104 (2009), arXiv:0905.2869 [hep-ex] .
  * Choi  _et al._ [2008] S. K. Choi _et al._ (Belle), Phys. Rev. Lett. 100, 142001 (2008), arXiv:0708.1790 [hep-ex] .
  * Ader  _et al._ [1982] J. P. Ader, J. M. Richard, and P. Taxil, Phys. Rev. D25, 2370 (1982).
  * Ballot and Richard [1983] J. l. Ballot and J. M. Richard, Phys. Lett. B123, 449 (1983).
  * Jaffe [1977] R. L. Jaffe, Phys.Rev. D15, 267 (1977).
  * Aaij  _et al._ [2021b] R. Aaij _et al._ (LHCb), (2021b), arXiv:2109.01038 [hep-ex] .
  * Aaij  _et al._ [2021c] R. Aaij _et al._ (LHCb), (2021c), arXiv:2109.01056 [hep-ex] .
  * Lüscher [1991] M. Lüscher, Nucl. Phys. B 354, 531 (1991).
  * Bicudo  _et al._ [2020] P. Bicudo, M. Cardoso, N. Cardoso, and M. Wagner, Phys. Rev. D 101, 034503 (2020), arXiv:1910.04827 [hep-lat] .
  * Bicudo  _et al._ [2021] P. Bicudo, N. Cardoso, L. Müller, and M. Wagner, Phys. Rev. D 103, 074507 (2021), arXiv:2008.05605 [hep-lat] .
  * Bali  _et al._ [2005] G. S. Bali, H. Neff, T. Duessel, T. Lippert, and K. Schilling (SESAM), Phys. Rev. D 71, 114513 (2005), arXiv:hep-lat/0505012 .
  * Bulava  _et al._ [2019] J. Bulava, B. Hörz, F. Knechtli, V. Koch, G. Moir, C. Morningstar, and M. Peardon, Phys. Lett. B793, 493 (2019), arXiv:1902.04006 [hep-lat] .
  * Bonati and Morlacchi [2020] C. Bonati and S. Morlacchi, Phys. Rev. D 101, 094506 (2020), arXiv:2003.07244 [hep-lat] .
  * Bicudo and Wagner [2013] P. Bicudo and M. Wagner, Phys. Rev. D 87, 114511 (2013), arXiv:1209.6274 [hep-ph] .
  * Brown and Orginos [2012] Z. S. Brown and K. Orginos, Phys. Rev. D 86, 114506 (2012), arXiv:1210.1953 [hep-lat] .
  * Li  _et al._ [2020] Q. Li, M.-S. Liu, Q.-F. Lü, L.-C. Gui, and X.-H. Zhong, Eur. Phys. J. C 80, 59 (2020), arXiv:1905.10344 [hep-ph] .
  * Liang  _et al._ [2020] W.-H. Liang, N. Ikeno, and E. Oset, Phys. Lett. B 803, 135340 (2020), arXiv:1912.03053 [hep-ph] .
  * Giron and Lebed [2020] J. F. Giron and R. F. Lebed, Phys. Rev. D 102, 014036 (2020), arXiv:2005.07100 [hep-ph] .
  * Meng and Chao [2008] C. Meng and K.-T. Chao, Phys. Rev. D 77, 074003 (2008), arXiv:0712.3595 [hep-ph] .
  * Simonov and Veselov [2009] Y. A. Simonov and A. I. Veselov, Phys. Lett. B 671, 55 (2009), arXiv:0805.4499 [hep-ph] .
  * Voloshin [2012] M. B. Voloshin, Phys. Rev. D 85, 034024 (2012), arXiv:1201.1222 [hep-ph] .
  * Süngü  _et al._ [2019] J. Y. Süngü, A. Türkan, H. Dağ, and E. Veli Veliev, Adv. High Energy Phys. 2019, 8091865 (2019), arXiv:1809.07213 [hep-ph] .
  * Mizuk  _et al._ [2019a] R. Mizuk _et al._ (Belle), JHEP 10, 220 (2019a), arXiv:1905.05521 [hep-ex] .
  * Wang [2019] Z.-G. Wang, Chin. Phys. C 43, 123102 (2019), arXiv:1905.06610 [hep-ph] .
  * Ali  _et al._ [2020] A. Ali, L. Maiani, A. Y. Parkhomenko, and W. Wang, Phys. Lett. B 802, 135217 (2020), arXiv:1910.07671 [hep-ph] .
  * Tarrús Castellà [2020] J. Tarrús Castellà, AIP Conf. Proc. 2249, 020008 (2020), arXiv:1908.05179 [hep-ph] .
  * Chen  _et al._ [2020] B. Chen, A. Zhang, and J. He, Phys. Rev. D 101, 014020 (2020), arXiv:1910.06065 [hep-ph] .
  * Brambilla  _et al._ [2020] N. Brambilla, S. Eidelman, C. Hanhart, A. Nefediev, C.-P. Shen, C. E. Thomas, A. Vairo, and C.-Z. Yuan, Phys. Rept. 873, 1 (2020), arXiv:1907.07583 [hep-ex] .
  * Choi  _et al._ [2003] S. K. Choi _et al._ (Belle), Phys. Rev. Lett. 91, 262001 (2003), arXiv:hep-ex/0309032 .
  * Acosta  _et al._ [2004] D. Acosta _et al._ (CDF), Phys. Rev. Lett. 93, 072001 (2004), arXiv:hep-ex/0312021 .
  * Hazi [1979] A. U. Hazi, Phys. Rev. A 19, 920 (1979).
  * Ashton  _et al._ [1983] C. J. Ashton, M. S. Child, and J. M. Hutson, The Journal of Chemical Physics 78, 4025 (1983), https://doi.org/10.1063/1.445128 .
  * Hagino  _et al._ [2020] K. Hagino, H. Sagawa, S. Kanaya, and A. Odahara, PTEP 2020, 023D01 (2020), arXiv:1903.01634 [nucl-th] .
  * Godfrey and Isgur [1985a] S. Godfrey and N. Isgur, Phys. Rev. D 32, 189 (1985a).
  * Sommer [2014] R. Sommer, PoS LATTICE2013, 015 (2014), arXiv:1401.3270 [hep-lat] .
  * Mizuk  _et al._ [2019b] R. Mizuk _et al._ (Belle),  _13th International Workshop on Heavy Quarkonium (QWG 2019), Torino, Italy, May 13-17, 2019_ , JHEP 10, 220 (2019b), arXiv:1905.05521 [hep-ex] .
  * Godfrey and Isgur [1985b] S. Godfrey and N. Isgur, Phys. Rev. D 32, 189 (1985b).
  * Padmanath  _et al._ [2015] M. Padmanath, C. B. Lang, and S. Prelovsek, Phys. Rev. D 92, 034501 (2015), arXiv:1503.03257 [hep-lat] .
  * Kang and Oller [2017] X.-W. Kang and J. A. Oller, Eur. Phys. J. C 77, 399 (2017), arXiv:1612.08420 [hep-ph] .
  * Micu [1969] L. Micu, Nucl. Phys. B 10, 521 (1969).
  * Le Yaouanc  _et al._ [1973] A. Le Yaouanc, L. Oliver, O. Pene, and J. C. Raynal, Phys. Rev. D 8, 2223 (1973).
  * Kokoski and Isgur [1987] R. Kokoski and N. Isgur, Phys. Rev. D 35, 907 (1987).
  * van Beveren  _et al._ [1986] E. van Beveren, T. A. Rijken, K. Metzger, C. Dullemond, G. Rupp, and J. E. Ribeiro, Z. Phys. C 30, 615 (1986), arXiv:0710.4067 [hep-ph] .
  * Bicudo and Ribeiro [1990] P. J. d. A. Bicudo and J. E. F. T. Ribeiro, Phys. Rev. D 42, 1635 (1990).
  * Bruschini and González [2020] R. Bruschini and P. González, Phys. Rev. D 102, 074002 (2020), arXiv:2007.07693 [hep-ph] .
  * Bicudo  _et al._ [2017] P. Bicudo, J. Scheunert, and M. Wagner, Phys. Rev. D 95, 034502 (2017), arXiv:1612.02758 [hep-lat] .
  * Bali  _et al._ [1997] G. S. Bali, K. Schilling, and A. Wachter, Phys. Rev. D 56, 2566 (1997), arXiv:hep-lat/9703019 .
  * Brambilla  _et al._ [2001] N. Brambilla, A. Pineda, J. Soto, and A. Vairo, Phys. Rev. D 63, 014023 (2001), arXiv:hep-ph/0002250 .
  * Pineda and Vairo [2001] A. Pineda and A. Vairo, Phys. Rev. D 63, 054007 (2001), [Erratum: Phys.Rev.D 64, 039902 (2001)], arXiv:hep-ph/0009145 .
  * Brambilla  _et al._ [2005] N. Brambilla, A. Pineda, J. Soto, and A. Vairo, Rev. Mod. Phys. 77, 1423 (2005), arXiv:hep-ph/0410047 .
  * Koma  _et al._ [2006] Y. Koma, M. Koma, and H. Wittig, Phys. Rev. Lett. 97, 122003 (2006), arXiv:hep-lat/0607009 .
  * Koma and Koma [2012] Y. Koma and M. Koma, PoS LATTICE2012, 140 (2012), arXiv:1211.6795 [hep-lat] .