$u d \bar{b} \bar{b}$ tetraquark resonances with lattice QCD potentials and the Born-Oppenheimer approximation

TL;DR

This paper predicts a new tetraquark resonance with specific quantum numbers and decay properties using lattice QCD potentials and the Born-Oppenheimer approximation, expanding understanding of exotic hadron states.

Contribution

It introduces a novel lattice QCD and Born-Oppenheimer approach to identify and characterize a tetraquark resonance with specific quantum numbers and decay width.

Findings

Predicted a tetraquark resonance at 10576 MeV with a decay width of about 112 MeV.

Identified a resonance with quantum numbers I(JP)=0(1-).

Extended previous bound state predictions to include resonant states with angular momentum l=1.

Abstract

We study tetraquark resonances with lattice QCD potentials computed for a static bbar bbar pair in the presence of two lighter quarks u d, the Born-Oppenheimer approximation and the emergent wave method. As a proof of concept we focus on the system with isospin I = 0, but consider different relative angular momenta l of the heavy quarks bbar bbar. For l=0 a bound state has already been predicted with quantum numbers I(JP) = 0(1+). Exploring various angular momenta we now compute the phase shifts and search for S and T matrix poles in the second Riemann sheet. We predict a tetraquark resonance for l =1, decaying into two B mesons, with quantum numbers I(JP) = 0(1-), mass m = 10 \, 576^{+4}_{-4} MeV} and decay width Gamma = 112^{+90}_{-103} MeV.

Figures

Tables

Table 1. Table 1: Parameters α 𝛼 \alpha and d 𝑑 d of the potential of Eq. ( 3 ) for two static antiquarks Q ¯ ​ Q ¯ ¯ 𝑄 ¯ 𝑄 \bar{Q}\bar{Q} , in the presence of two light quarks q ​ q 𝑞 𝑞 qq with quantum numbers I 𝐼 I and j 𝑗 j , as determined in Bicudo  et al. ( 2016 ) .

$I$	$j$	$\alpha$	$d$ in fm
$0$	$0$	${0.34}_{-0.03}^{+0.03}$	${0.45}_{-0.10}^{+0.12}$
$1$	$1$	${0.29}_{-0.06}^{+0.05}$	${0.16}_{-0.02}^{+0.05}$

Equations

$$\mathcal{O}[f,f^{\prime },\Gamma ,\tilde{\Gamma }]=(\mathcal{C}\Gamma {)}_{AB}(\mathcal{C}\tilde{\Gamma }{)}_{CD}$$

$$\displaystyle=\Big{(}\bar{Q}_{C}^{a}(\mathbf{r}_{1})\psi_{A}^{(f)a}(\mathbf{r}_{1})\Big{)}\Big{(}\bar{Q}_{D}^{b}(\mathbf{r}_{2})\psi_{B}^{(f^{\prime})b}(\mathbf{r}_{2})\Big{)}\ ,$$

$$⟨\Omega |\mathcal{O}[f,f^{\prime },\Gamma ,\tilde{\Gamma }{]}^{†}(t)\mathcal{O}[f,f^{\prime },\Gamma ,\tilde{\Gamma }](0)|\Omega ⟩.$$

$$V(r)=-\frac{\alpha }{r}{e}^{-r^2/d^2}+V_0.$$

$$\Big{(}H_{0}+V(r)\Big{)}\Psi=E\Psi\ .$$

$$\Psi ={\Psi }_0+X,$$

$$H_0{\Psi }_0=E{\Psi }_0,$$

$$\Big{(}H_{0}+V(r)-E\Big{)}X=-V(r)\Psi_{0}\ .$$

$$H=H_0+V(r)=-\frac{{\hslash }^2}{2\mu }△+V(r)$$

$${\Psi }_0=e^{i\mathbf{k}\cdot \mathbf{r}}=\sum _l(2l+1)i^l{j}_l(kr)P_l(\hat{\mathbf{k}}\cdot \hat{\mathbf{r}}),$$

$$X=\sum _l(2l+1)i^l\frac{{\chi }_l(r)}{kr}{P}_l(\hat{\mathbf{k}}\cdot \hat{\mathbf{r}}).$$

$$\displaystyle\bigg{(}-\frac{\hbar^{2}}{2\mu}\frac{d^{2}}{dr^{2}}+\frac{l(l+1)}{2\mu r^{2}}+V(r)-E\bigg{)}\chi_{l}(r)=$$

$$=-V(r)kr{j}_l(kr).$$

$$\frac{{\chi }_l(r)}{kr}=i{t}_l{h}_l^{(1)}(kr),$$

$$s_l\equiv 1+2i{t}_l=e^{2i{\delta }_l}.$$

$${\chi }^2=\sum _{r=r_{\text{min}},...,r_{\text{max}}}{\left(\frac{V(r)-V^{\text{lat}}(r)}{\Delta V^{\text{lat}}(r)}\right)}^2$$

Full text

I Introduction

A long standing problem in particle physics is to understand exotic hadrons, i.e. hadrons which have a structure more complicated than a quark-antiquark pair or a triplet of quarks Jaffe (1977). The problem of identifying exotic hadrons, say tetraquarks, pentaquarks, hexaquarks, hybrids or glueballs – expected since the onset of QCD – turned out to be much harder than initially expected Bicudo and Cardoso (2016). The observed candidates are resonances high in the spectrum, not only difficult to observe, but also technical to address in quark or hadron models. They possibly require the development of new techniques, potentially relevant to other areas of physics, to be studied theoretically from first principles, e.g. with lattice QCD Wu et al. (2014); Morningstar et al. (2016).

Our main motivation is to investigate tetraquarks by combining lattice QCD and quantum mechanics techniques. We specialize in systems with two heavy antiquarks, which are expected to form bound states, when sufficiently heavy Ader et al. (1982); Ballot and Richard (1983); Heller and Tjon (1987); Carlson et al. (1988); Lipkin (1986); Brink and Stancu (1998); Gelman and Nussinov (2003); Vijande et al. (2004); Janc and Rosina (2004); Cohen and Hohler (2006); Vijande et al. (2007). The starting point are potentials of two static antiquarks in the presence of two light quarks, which can be computed with state of the art lattice QCD techniques (cf. e.g. Detmold et al. (2007); Wagner (2010); Bali and Hetzenegger (2010); Wagner (2011); Brown and Orginos (2012); Bicudo et al. (2016)). If the masses of the two heavy quarks are much larger than the scale of QCD, which is the case for two $\bar{b}$ quarks, their dynamics can then be described by a quantum mechanical Hamiltonian with the aforementioned lattice QCD potentials. This two-step approach is the Born-Oppenheimer approximation Born and Oppenheimer (1927). Using this approach, a $ud\bar{b}\bar{b}$ tetraquark bound state with quantum numbers $I(J^{P})=0(1^{+})$ has recently been predicted Bicudo and Wagner (2013); Brown and Orginos (2012); Bicudo et al. (2015, 2016, 2017) and confirmed by a lattice QCD computation with four quarks of finite mass Francis et al. (2016). So far, however, resonances have not been studied in this framework.

Notice there are two classes of double-heavy tetraquarks. The tetraquarks with one heavy quark and one heavy antiquark including the $Z_{c}$ and $Z_{b}$ are easier to detect experimentally. Their observation at Belle Bondar et al. (2012); Liu et al. (2013); Chilikin et al. (2014), Cleo-C Xiao et al. (2013), BESIII Ablikim et al. (2013a, 2014a, b, 2014b, 2014c) and LHCb Aaij et al. (2014) collaborations turned tetraquarks into a main highlight of particle physics in recent years. But since they have more coupled channels we opt here to study tetraquarks with two heavy antiquarks (or quarks), which are theoretically simpler. This “theoretical simplicity” is convenient for a first study of resonances with lattice QCD potentials. Moreover, with the recent observation at LHCb of hadronic systems with two heavy quarks Maciuła et al. (2016); Aaij et al. (2017) we expect this second class of tetraquarks to be observed in the near future.

In this work we extent the previous Born-Oppenheimer studies with lattice QCD potentials, reviewed in Section II. We utilize the emergent wave method, a technique from scattering theory detailed in Section III, to compute phase shifts, S and T matrix poles in the second Riemann sheet and the corresponding resonance masses and decay widths. For the first time, we apply this technique with lattice QCD potentials, and our results are presented in Section IV. We conclude in Section V.

II Lattice QCD potentials of two static antiquarks in the presence of two light quarks and prediction of a stable $ud\bar{b}\bar{b}$ tetraquark

In preceding papers we have computed potentials $V(r)$ of two static antiquarks $\bar{Q}\bar{Q}$ in the presence of two light quarks $qq$ using lattice QCD. The computations have been carried out for many different quantum numbers including light flavor combinations $qq$ with $q\in\{u,d,s,c\}$, parity $P$ and total angular momentum of the light quarks and gluons $j$ (cf. e.g. Bicudo et al. (2015, 2016)). There are both attractive and repulsive channels. Most promising with respect to the existence of tetraquark bound states or resonances are attractive potentials with light quarks $q\in\{u,d\}$, since they are rather wide and deep. There are two such potentials, with quantum numbers $(I=0,j=0)$ and $(I=1,j=1)$, where $I$ denotes isospin.

We have used creation operators

[TABLE]

where $r=|\mathbf{r}_{2}-\mathbf{r}_{1}|$, $a,b$ denote color and $A,B,C,D$ spin indices and $\psi^{(f)}\psi^{(f^{\prime})}=ud-du$ for $I=0$ and $\psi^{(f)}\psi^{(f^{\prime})}\in\{uu,ud+du,dd\}$ for $I=1$. For the $(I=0,j=0)$ potential $\Gamma=(1+\gamma_{0})\gamma_{5}$, while for the $(I=1,j=1)$ potential $\Gamma=(1+\gamma_{0})\gamma_{j}$ ($j=1,2,3$). Since the potentials are independent of the static quark spins, one can choose arbitrarily $\tilde{\Gamma}\in\{(1-\gamma_{0})\gamma_{5},(1-\gamma_{0})\gamma_{j}\}$. As usual in lattice QCD hadron spectroscopy we have extracted the potentials from the asymptotic exponential decay in the temporal separation $t$ of correlation functions

[TABLE]

Example plots for lattice spacing $a\approx 0.079\,\textrm{fm}$ and $u/d$ quark masses corresponding to a pion mass $m_{\pi}\approx 340\,\textrm{MeV}$ are shown in Fig. 1.

Since it is known that the existence of a stable tetraquark as well as its binding energy exhibits a sizable dependence on the light quark mass Bicudo et al. (2015), we have performed computations of the potentials for three different $u/d$ quark masses corresponding to $m_{\pi}\in\{340\,\textrm{MeV},480\,\textrm{MeV},650\,\textrm{MeV}\}$. Then we have used these results to extrapolate to the physical $u/d$ quark mass corresponding to $m_{\pi}=140\,\textrm{MeV}$. Moreover, we have crudely estimated systematic errors due to the finite lattice spacing $a\approx 0.079\,\textrm{fm}$ by performing the computations with two different Wilson twisted mass lattice QCD discretizations. We have found that discretization errors are negligible compared to the current statistical uncertainties (for more details cf. Bicudo et al. (2016)). Similarly, effects due to the finite spatial volume of the lattice are expected to be negligible as well.

To search for bound states and resonances we parameterize the potentials by a screened Coulomb potential,

[TABLE]

This ansatz is inspired by one-gluon exchange at small $\bar{Q}\bar{Q}$ separations $r$ and a screening of the Coulomb potential due to the formation of two $B$ mesons at large $r$, as illustrated in Fig. 2. The ansatz, even though phenomenologically motivated, is consistent with our lattice QCD results, which are based on first principles, i.e. a fit of (3) to the lattice QCD data yields an acceptable \chi^{2}/\textrm{dof}\raisebox{-2.15277pt}{,\stackrel{{\scriptstyle<}}{{\scriptstyle\sim}},}1. Vice versa, parameterizing the lattice QCD data by using ansätze different from (3) leads e.g. to similar results for masses of tetraquark bound states. The values of the two parameters $\alpha$ and $d$ as determined in Bicudo et al. (2016) are listed in Table 1. Clearly, the $(I=0,j=0)$ potential is more attractive than the $(I=1,j=1)$ potential. Note that there is also an uncertainty associated with the lattice spacing, $a=0.079(3)\,\textrm{fm}$ (cf. Baron et al. (2010)), which is not included in the parameter $d$ in Table 1. We investigate the effect of this uncertainty at the end of our analysis in section IV.

Finally we have applied the Born-Oppenheimer approximation, where Eq. (3) is used as a potential for two heavy antiquarks, i.e. $\bar{b}\bar{b}$, in the presence of two light quarks $ud$ or for two heavy-light mesons, i.e. $B^{(\ast)}B^{(\ast)}$. Solving the Schrödinger equation for the $(I=0,j=0)$ potential and angular momentum $l=0$ of the two $\bar{b}$ quarks a bound state has been predicted with binding energy $90^{+43}_{-36}\textrm{MeV}$ and quantum numbers $I(J^{P})=0(1^{+})$ Bicudo et al. (2016).

The use of the Born-Oppenheimer approximation entails a systematic error from quantizing the $\bar{b}\bar{b}$ system with the kinetic energy only. The kinetic energy naturally emerges in the next to leading term in a non-relativistic series expansion. However, the spin dependent terms of the potential are of the same non-relativistic expansion order of the kinetic energy and so far we have not taken them into account directly. Nevertheless, in Ref. Bicudo et al. (2017) the spin effects have been estimated and they have little effect on the binding energy of the tetraquark. Finally, very recent computations in lattice QCD with non-relativistic bottom quarks, which account for both the kinetic and spin effects, confirm our previous results for the binding energy, obtained with the Born-Oppenheimer approximation Junnarkar et al. (2017). Thus we expect that the use of the Born-Oppenheimer approximation is adequate for our study.

III The emergent wave method

We now summarize the emergent wave method, explained in detail for instance in Ref. Bicudo and Cardoso (2016), which is suited to study phase shifts and resonances. Let us consider the same Schrödinger equation utilized in the bound state study,

[TABLE]

The first step is to split the wave function into two parts,

[TABLE]

where $\Psi_{0}$ is the incident wave, a solution of the free Schrödinger equation,

[TABLE]

and $X$ is the emergent wave. Substituting Eq. (5) into Eq. (4) and using Eq. (6) we obtain

[TABLE]

For any energy $E$ we can use this equation to calculate the emergent wave $X$ by providing the corresponding $\Psi_{0}$ and fixing the appropriate boundary conditions. From the asymptotic behaviour of $X$ we then determine the phase shifts, the S matrix and the T matrix.

The problem can be continued to complex energies in a straightforward way and we can, therefore, find the poles of the S matrix and the T matrix in the complex plane. We identify a resonance with a pole, when located in the second Riemann sheet at $m-i\Gamma/2$, where $m$ is the mass and $\Gamma$ is the decay width of the resonance.

III.1 Partial wave decomposition

The Hamiltonian describing the two heavy antiquarks $\bar{b}\bar{b}$ at vanishing total momentum, i.e. in the rest frame of the system, is

[TABLE]

with reduced mass $\mu=M/2$, where $M=5\,280\,\textrm{MeV}$ is the mass of the $B$ meson from the PDG Olive et al. (2014). For simplicity we omit the additive constant $2M$ in Eq. (8), i.e. all resulting energy eigenvalues are energy differences with respect to $2M$. We consider an incident plane wave $\Psi_{0}=e^{i\mathbf{k}\cdot\mathbf{r}}$, which can be expressed as a sum of spherical waves,

[TABLE]

where $j_{l}$ are spherical Bessel functions, $P_{l}$ are Legendre polynomials and the relation between energy and momentum is $\hbar k=\sqrt{2\mu E}$. For a spherically symmetric potential $V(r)$ as in Eq. (3) and an incident wave $\Psi_{0}=e^{i\mathbf{k}\cdot\mathbf{r}}$ the emergent wave $X$ can also be expanded in terms of Legendre polynomials $P_{l}$,

[TABLE]

Inserting Eq. (9) and Eq. (10) into Eq. (7) leads to a set of ordinary differential equations for $\chi_{l}$,

[TABLE]

III.2 Solving the differential equations for the emergent wave

The potentials $V(r)$, Eq. (3), are exponentially screened, i.e. $V(r)\approx 0$ for $r\geq R$, where $R\gg d$. For large separations $r\geq R$ the emergent wave is, hence, a superposition of outgoing spherical waves, i.e.

[TABLE]

where $h_{l}^{(1)}$ are the spherical Hankel functions of first kind.

Our aim is now to compute the complex prefactors $t_{l}$, which will eventually lead to the phase shifts. To this end we solve the ordinary differential equation (11). The corresponding boundary conditions are the following:

• •

At $r=0$: $\chi_{l}(r)\propto r^{l+1}$.

• •

For $r\geq R$: Eq. (12).

Note that the boundary condition for $r\geq R$ depends on $t_{l}$. For a given value of the energy $E$ this boundary condition is only fulfilled for a specific corresponding value of $t_{l}$. In other words the boundary condition for $r\geq R$ fixes $t_{l}$ as a function of $E$.

The numerical solution of the differential Eq. (11) is rather straightforward. To check our results and to exclude any numerical artefacts we implemented two different approaches: (1) a fine uniform discretization of the interval $[0,R]$, which reduces the differential equation to a large set of linear equations, which can be solved rather efficiently, since the corresponding matrix is tridiagonal; (2) a standard 4-th order Runge-Kutta shooting method.

III.3 Phase shifts and $\mathbf{S}$ and $\mathbf{T}$ matrix poles

The quantity $t_{l}$ is a T matrix eigenvalue (cf. standard textbooks on quantum mechanics and scattering, e.g. Merzbacher (1998)). From $t_{l}$ we can calculate the phase shift $\delta_{l}$ and also read off the corresponding S matrix eigenvalue $s_{l}$ 111 At large distances $r\geq R$, the radial wavefunction is $kr[j_{l}(kr)+i\,t_{l}h_{l}^{(1)}(kr)]=(kr/2)[h_{l}^{(2)}(kr)+e^{2i\delta_{l}}h_{l}^{(1)}(kr)]$. ,

[TABLE]

Moreover, note that both the S matrix and the T matrix are analytical in the complex plane. They are well-defined for complex energies $E$. Thus, our numerical method can as well be applied to solve the differential Eq. (11) for complex $E$. We find the S and T matrix poles by scanning the complex plane $(\textrm{Re}(E),\textrm{Im}(E))$ and applying Newton’s method to find the roots of $1/t_{l}(E)$. The poles of the S and the T matrix correspond to complex energies of resonances. Note the resonance poles must be in the second Riemann sheet with a negative imaginary part both for the energy $E$ and the momentum $k$.

IV Results for phase shifts, $\mathbf{S}$ matrix and $\mathbf{T}$ matrix poles and resonances

We first consider the more attractive $ud\bar{b}\bar{b}$ potential corresponding to isospin $I=0$ and light spin $j=0$ (cf. Sec. II). We compute $t_{l}$ and via Eq. (13) the phase shift $\delta_{l}$ for real energy $E$ and angular momenta $l=0,1,2,\ldots$ A very clear signal for a resonance would be a fast increase of the phase shift $\delta_{l}$ as a function of $E$ from [math] to $\approx\pi$, almost like a step function. However, we do not find such a pronounced increase (cf. Fig. 3). Thus, we must search more thoroughly for possibly existing resonances.

Starting with angular momentum $l=1$ we first search for clear resonance signals by making the potential more and more attractive. We increase the parameter $\alpha$, while keeping the parameter $d=0.45\,\textrm{fm}$ fixed, to preserve the scale of the potential. The corresponding results for the phase shift $\delta_{1}$ are shown in Fig. 4. Indeed, for \alpha\raisebox{-2.15277pt}{,\stackrel{{\scriptstyle>}}{{\scriptstyle\sim}},}0.65 we find clear resonances with $\delta_{1}$ increasing from [math] to $\approx\pi$. Then, for $\alpha=0.72$, we find a bound state, since the phase shift $\delta_{1}$ starts at $\pi$ and decreases monotonically to [math], when increasing the energy $E$. However, from these phase shifts it is not clear, for which values of $\alpha$ a resonance exists or not, i.e. it is not possible to say, whether there is a resonance for e.g. $\alpha\approx 0.50$ or even for the physical $\alpha=0.34$.

Thus, we search directly for poles of the T matrix eigenvalues $t_{l}$. With this technique we clearly find a pole for angular momentum $l=1$ and physical values of the parameters, $\alpha=0.34$ and $d=0.45\,\textrm{fm}$. We show this pole in Fig. 5 by plotting $t_{1}$ as a function of the complex energy $E$. The pole is clearly visible as a sharp peak.

To understand the dependence of the resonance pole on the shape of the potential, we again scan different values of the parameter $\alpha$ and determine each time the pole of the eigenvalue $t_{1}$ of the T matrix. We show the trajectory of the pole corresponding to a variation of $\alpha$ in the complex plane $(\textrm{Re}(E),\textrm{Im}(E))$ in Fig. 6. Indeed, starting with $\alpha=0.21$ we find a pole. This confirms our prediction of a resonance for angular momentum $l=1$ and physical values of the parameters, $\alpha=0.34$ and $d=0.45\,\textrm{fm}$.

Finally we perform a detailed statistical and systematic error analysis of the pole of $t_{1}$ and the corresponding values $(\textrm{Re}(E),\textrm{Im}(E))$. We use the same analysis method as for our previous study of the bound state for $l=0$, cf. Bicudo et al. (2015). To parameterize the lattice QCD data for the potentials, $V^{\textrm{lat}}(r)$, discussed in Section II, we perform uncorrelated $\chi^{2}$ minimizing fits with the ansatz (3). To this end we minimize the expression

[TABLE]

with respect to the parameters $\alpha$, $d$ and $V_{0}$ ($\Delta V^{\textrm{lat}}(r)$ denote the corresponding statistical errors). To quantify systematic errors, we perform a large number of fits, where we vary the following parameters:

• •

The range of temporal separations $t_{\textrm{min}}\leq t\leq t_{\textrm{max}}$ of the correlation function (2), where $V^{\textrm{lat}}(r)$ is read off, according to

• –

$t_{\textrm{max}}-t_{\textrm{min}}\geq a$,

• –

$4a\leq t_{\textrm{min}}$, $t_{\textrm{max}}\leq 9a$

($a\approx 0.079\,\textrm{fm}$ is the lattice spacing).

• •

The range of spatial $\bar{b}\bar{b}$ separations $r_{\textrm{min}}\leq r\leq r_{\textrm{max}}$ considered in the $\chi^{2}$ minimizing fits to determine the parameters $\alpha$, $d$ and $V_{0}$ according to

• –

$r_{\textrm{min}}\in\{2a,3a\}$,

• –

$r_{\textrm{max}}\in\{8a,9a,10a\}$.

We obtain a large number of different, but similar potential parameterizations $V(r)$ characterized by sets of values for $\alpha$, $d$ and $V_{0}$. For each potential parameterization we determine the position of the pole of $t_{1}$, i.e. $(\textrm{Re}(E),\textrm{Im}(E))$ as discussed above and shown as a cloud of blue points in Fig. 6. For both $\textrm{Re}(E)$ and $\textrm{Im}(E)$ we construct a distribution by considering all corresponding results weighted by $\exp(-\chi^{2}/\textrm{dof})$ with $\chi^{2}$ from Eq. (14). The central values of $\textrm{Re}(E)$ and $\textrm{Im}(E)$ are then defined as the medians of the corresponding distributions and the lower/upper systematic uncertainties are given by the differences of the 16th/84th percentiles to the medians. To also include statistical errors, we compute the jackknife errors of the medians of $\textrm{Re}(E)$ and $\textrm{Im}(E)$ and add them in quadrature to the corresponding systematic uncertainties. With our combined statistical and systematic error analysis we find a resonance energy $\textrm{Re}(E)=17^{+4}_{-4}\,\textrm{MeV}$ and a decay width $\Gamma=-2\textrm{Im}(E)=112^{+90}_{-103}\,\textrm{MeV}$. Using the Pauli principle and considering the symmetry of the quarks with respect to colour, flavour, spin and their spatial wave function one can determine the quantum numbers of the resonance, which are $I(J^{P})=0(1^{-})$. The resonance will decay into two $B$ mesons and, hence, its mass is $m=2M+\textrm{Re}(E)=10\,576^{+4}_{-4}\,\textrm{MeV}$.

Note that there is also an uncertainty associated with the lattice spacing, $a=0.079(3)\,\textrm{fm}$ (cf. Ref. Baron et al. (2010) for details), which has not been taken into account so far. We have investigated the impact of this uncertainty on our final results for the resonance energy $\textrm{Re}(E)$ and the decay width $\Gamma$. We have found that both quantities exhibit only a mild dependence on the lattice spacing $a$ and the propagation of the uncertainty of $a$ has a negligible effect on the results for $\textrm{Re}(E)$ and $\Gamma$ quoted above within the current combined systematic and statistical errors.

In what concerns angular momenta $l\neq 1$, we find no clear signal for a resonance pole (except for the bound state pole for $l=0$). We also find no poles for any $l$ in the less attractive case of $(I=1,j=1)$.

V Conclusions and outlook

As a case study for the investigation of resonances above the $BB$ meson pair threshold, we have explored the $ud\bar{b}\bar{b}$ four-quark system. We have utilized lattice QCD potentials computed for two static antiquarks in the presence of two light quarks, the Born-Oppenheimer approximation and the emergent wave method for the $BB$ system. First we have computed scattering phase shifts. Then we have performed the analytic continuation of the S matrix and the T matrix to the second Riemann sheet and have searched for poles as signals of resonances.

From these results we have predicted a new resonance, with quantum numbers $I(J^{P})=0(1^{-})$. Performing a careful statistical and systematic error analysis has led to a resonance mass $m=10\,576^{+4}_{-4}\,\textrm{MeV}$ and a decay width $\Gamma=112^{+90}_{-103}\,\textrm{MeV}$.

Acknowledgements.

We acknowledge useful conversations with K. Cichy. P.B. acknowledges the support of CFTP (grant FCT UID/FIS/00777/2013) and is thankful for hospitality at the Institute of Theoretical Physics of Johann Wolfgang Goethe-University Frankfurt am Main. M.C. acknowledges the support of CFTP and the FCT contract SFRH/BPD/73140/2010. M.W. acknowledges support by the Emmy Noether Programme of the DFG (German Research Foundation), grant WA 3000/1-1. This work was supported in part by the Helmholtz International Center for FAIR within the framework of the LOEWE program launched by the State of Hesse. Calculations on the LOEWE-CSC and on the 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.