Tetraquark resonances computed with static lattice QCD potentials and scattering theory

We study tetraquark resonances with lattice QCD potentials computed for two static quarks and two dynamical quarks, the Born-Oppenheimer approximation and the emergent wave method of scattering theory. As a proof of concept we focus on systems with isospin $I = 0$, but consider different relative angular momenta $l$ of the heavy $b$ quarks. We compute the phase shifts and search for $\mbox{S}$ and $\mbox{T}$ matrix poles in the second Riemann sheet. We predict a new tetraquark resonance for $l = 1$, decaying into two $B$ mesons, with quantum numbers $I(J^P) = 0(1^-)$, mass $m = 10576_{-4}^{+4} \, \textrm{MeV}$ and decay width $\Gamma = 112_{-103}^{+90} \, \textrm{MeV}$.


Introduction
A long standing problem in QCD is to understand exotic hadrons. In this work we specialize in tetraquark systems with two heavy antiquarksbb and two lighter quarks qq, where q ∈ {u, d, s, c}. The existence of bound states has been extensively investigated in the recent past by combining static lattice QCD potentials and the Born-Oppenheimer approximation. A stable udbb tetraquark with quantum numbers I(J P ) = 0(1 + ) has been predicted [1][2][3][4][5][6][7][8][9] and been confirmed by similar computations using four quarks of finite mass [10]. Here we extend our investigation by including a new technique from scattering theory, the emergent wave method [11], and search for possibly existing tetraquark resonances (cf. also [12] for more details).

Lattice QCD potentials of two static antiquarksQQ in the presence of two lighter quarks qq
In a first step we have computed potentials V(r) of two static antiquarksQQ in the presence of two lighter quarks qq, where q ∈ {u, d, s, c}, using lattice QCD [2,4]. There are both attractive and repulsive channels. Most promising with respect to the existence of stable tetraquarks or tetraquark resonances are light quarks q ∈ {u, d} together with (I = 0, j = 0) or (I = 1, j = 1), where I denotes isospin and j light total angular momentum. The corresponding potentials are not only attractive, but also rather wide and deep [7].
Speaker, e-mail: bicudo@tecnico.ulisboa.pt  The lattice QCD results for the potentials can be parametrized by a screened Coulomb potential, inspired by one-gluon exchange at smallQQ separations r and a screening of the Coulomb potential by the two B mesons at large r (cf. Figure 1). Clearly, the (I = 0, j = 0) potential is more attractive than the (I = 1, j = 1) potential as can be seen in Figure 2. Numerical results for the parameters α and d are collected in Table 1. The potential parametrization is then used in the Schrödinger equation for the relative coordinate of the heavy antiquarksbb ≡QQ to search for either for bound states (cf. [5,[7][8][9]) or for resonances (cf. sections 3 and 4 ).

The emergent wave method
In this section, we summarize the emergent wave method, which is suited to compute phase shifts and resonances. More details can be found in Ref. [11].

Emergent and incident wavefunctions
We consider the Schrödinger equation used for studying bound states: The first step in the emergent wave method is to split the wave function of the Schrödinger Eq. into two parts, where Ψ 0 is the incident wave, a solution of the free Schrödinger equation, H 0 Ψ 0 = EΨ 0 , and X is the emergent wave. Inserting this in Eq.
(2), we obtain: For any energy E we calculate the emergent wave X by providing the corresponding Ψ 0 and fixing the appropriate boundary conditions. From the asymptotic behaviour of the emergent wave X we then determine the phase shifts δ l , the S matrix and the T matrix. Continuing to complex energies E ∈ C we find the poles of the S matrix and the T matrix in the complex plane. We identify a resonance with a pole of S in the second Riemann sheet at m − iΓ/2, where m is the mass and Γ is the resonance decay width.

Partial wave decomposition
The two heavy antiquarksbb at zero total angular momentum are described by the Hamiltonian: with reduced mass µ = M/2, where M = 5 280 MeV is the mass of the B meson from the PDG [13]. For simplicity we omit the additive constant 2M in Eq. (5), i.e. all resulting energy eigenvalues are energy differences with respect to 2M. We consider an incident plane wave Ψ 0 = e ik·r , which can be expressed as a sum of spherical waves, where j l are spherical Bessel functions, P l are Legendre polynomials and the relation between energy and momentum is k = 2µE. For a spherically symmetric potential V(r) as in Eq. (1) and an incident wave Ψ 0 = e ik·r , the emergent wave X can also be expanded in terms of Legendre polynomials P l , Inserting Eq. (6) and Eq. (7) into Eq. (4) leads to a set of ordinary differential equations for χ l :

Solving the differential equations for the emergent wave
The potentials V(r), Eq. (1), are exponentially screened, i.e. V(r) ≈ 0 for r ≥ R, where R d. For large separations r ≥ R the emergent wave is, hence, a superposition of outgoing spherical waves, i.e.
where h (1) l 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 (8). The corresponding boundary conditions are the following: • At r = 0: χ l (r) ∝ r l+1 .
The boundary condition for r ≥ R fixes t l as a function of E.
We solve it numerically, with two different numerical techniques 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.

Phase shifts and S and T matrix poles
The quantity t l is a T matrix eigenvalue (c.f. a standard textbook on quantum mechanics, e.g. [14]). For instance, at large distances r ≥ R, the radial wavefunction is  From t l we can calculate the phase shift δ l and also read off the corresponding S matrix eigenvalue s l , 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 ∈ C. Thus, our numerical method can as well be applied to solve the differential Eq. (8) for complex E ∈ C. We find the S and T matrix poles by scanning the complex plane (Re(E), 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.

Phase shifts
We first consider the udbb potential corresponding to isospin I = 0 and light spin j = 0 (cf. Section 3), since this potential is most attractive. We compute t l and via Eq. (11) the phase shift δ l for real energy E and angular momenta l = 0, 1, 2, . . .. We do not find a fast increase of the phase shift δ l as a function of the energy E which would clearly indicate a resonance (cf. Figure 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 attractive, i.e. we increase the parameter α. We keep the parameter d = 0.45 fm fixed here to preserve the scale of the potential. The corresponding results for the phase shift δ 1 are shown in Figure 4. Indeed, for α ≈ 0.65 we find clear resonances. For α = 0.72, we find a bound state, since the phase shift δ 1 starts at π and decreases monotonically to 0, when increasing the energy E. However, it is not clear from this observation, for which values of α a resonance exists or not.

Resonances as poles of the S and T matrices
To clearly identify a resonance, 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, α = 0.34 and d = 0.45 fm. We show this pole in Figure 5 by plotting t 1 as a function of the complex energy E. To understand the dependence of the resonance pole on the shape of the potential, we again scan different values of the parameter α 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 α in the complex plane (Re(E), Im(E)) in Figure 6. Indeed, starting with α = 0.21 we find a pole. This confirms our prediction of a resonance for angular momentum l = 1 and physical values of the parameters, α = 0.34 and d = 0.45 fm. In what concerns angular momenta l 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).

Statistical and systematic error analysis
Finally we perform a detailed statistical and systematic error analysis of the pole of t 1 and the corresponding values (Re(E), Im(E)) for angular momentum l = 1. We use the same analysis method as for our previous study of the bound state for l = 0, cf. [7]. To parametrize the lattice QCD data for the potentials, V lat (r), discussed in Section 3, we perform uncorrelated χ 2 minimizing fits with the ansatz of Eq. (1). To this end we minimize the expression ∆V lat (r) 2 (12) with respect to the parameters α, d and V 0 defined in Eq. (1) and in Refs. [2,4,15]. In Eq. (12), ∆V lat (r) denotes the corresponding statistical errors. To quantify systematic errors, we perform a large  Figure 6. Trajectory of the pole of the eigenvalue t 1 of the T matrix in the complex plane (Re(E), Im(E)), corresponding to a variation of parameter α. We also illustrate with a cloud of diamond points the systematic error [7] . number of fits. We vary the range of temporal separations t min ≤ t ≤ t max of the correlation function where V lat (r) is read off as well as the range of spatialbb separations r min ≤ r ≤ r max considered in the χ 2 minimizing fits to determine the parameters α, d and V 0 .
To also include statistical errors, we compute the jackknife errors of the medians of Re(E) and 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 Re(E) = 17 +4 −4 MeV and a decay width Γ = −2Im(E) = 112 +90 −103 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 + Re(E) = 10 576 +4 −4 MeV.

Summary and outlook
We 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 to search for udbb resonances. First we computed the scattering phase shifts of a BB meson pair. Then we performed the analytic continuation of the S matrix and the T matrix to the second Riemann sheet, where we have searched for poles ∈ C. From these results we have predicted a novel resonance with quantum numbers I(J P ) = 0(1 − ). Performing a careful statistical and systematic error analysis, we found a resonance mass m = 10 576 +4 −4 MeV and a decay width Γ = 112 +90 −103 MeV. For more details, please see our recent publication [12].
In the future we plan to address the experimentally observed quarkonia tetraquarks, including bb or cc heavy quarks [15], with our method.