Bottomonium resonances with I = 0 from lattice QCD correlation functions with static and light quarks ( 1 )

We discuss, how to study I = 0 quarkonium resonances decaying into pairs of heavy-light mesons using static potentials from lattice QCD. These static potentials can be obtained from a set of correlation functions containing both static and light quarks. As a proof of concept we focus on bottomonium with relative orbital angular momentum L = 0 of the b̄b pair corresponding to J = 0−+ and J = 1−−. We use static potentials from an existing lattice QCD string breaking study and compute phase shifts and T matrix poles for the lightest heavy-light meson-meson decay channel. We discuss our results in the context of corresponding experimental results, in particular for Υ(10860) and Υ(11020).

In this work we continue to use lattice QCD potentials and the Born-Oppenheimer approximation and focus on quarkonium bound states and resonances, which might be exotic, i.e. states containing a heavy quark and a heavy antiquark and possibly an additional light quarkantiquark pair. To study for example the experimentally observed Z b tetraquark resonances, it is necessary to extend the techniques introduced in [28] from a single channel to a coupled channel Schrödinger equation. One has to consider a confined quarkonium channelQQ and at least one scattering channelM M with two heavy-light mesons M =Qq andM =qQ. In the least complicated case of a single scattering channel the potential in the Schrödinger equation is a 2 × 2 matrix of the form as we will derive in detail in section II. It is important to note that the off-diagonal terms V mix (r) couple the quarkonium channel and the meson-meson channel. The consequence is that quarkonium bound states only exist below theM M threshold, whereas above this threshold all quarkonium states are resonances. Once we have set up an appropriate coupled channel Schrödinger equation we proceed as in Ref. [28] and compute phase shifts and T matrix poles, where the latter provide masses of bound states below threshold as well as resonance masses and decay widths above threshold. The main advantage of our approach is that it is in principle straightforward to consider several decay channels, since this is done in the framework of quantum mechanics in the form of a coupled channel Schrödiger equation. Each decay channel, however, requires the lattice QCD computation of static potentials with specific quantum numbers. Even though this might be time consuming and challenging, these static potentials are not only important for the approach we propose in this paper, but they will be of interest for theoretical hadron physics in general (as an example see the detailed discussion of static potentials for thebbud case in Ref. [9]). As mentioned above, it would be very interesting to study the experimentally observed Z b tetraquark resonances with I = 1. The lattice QCD computation of the potential matrix (1) for I = 1 is, however, very difficult (see e.g. Refs. [29,30]). Therefore, we decided to first explore the simpler I = 0 case, where corresponding high-quality lattice QCD potentials are provided by the string breaking computation of Ref. [31]. This allows us to study the established quarkonium states η b (1S), Υ(1S), Υ(2S), Υ(3S) and Υ(4S), possibly also Υ(10860) and Υ(11020) and to predict additional quarkonium resonances not yet observed experimentally.
This paper is organized as follows. In section II we detail the theoretical basics of our approach. We start by discussing quantum numbers of quarkonium and of pairs of heavy-light mesons for I = 0. Then we show, how to set up a corresponding coupled channel Schrödiger equation in a consistent way, and explain, how the potential matrix is related to static potentials from QCD, which can be computed with lattice QCD. We also discuss the boundary conditions of the wave function, which are appropriate for a coupled channel scattering problem. Moreover, we specialize the coupled channel Schrödiger equation for the specific case of relative orbital angular momentum L = 0 for the two heavy quarks, which will be the starting point for all numerical results presented later in the paper. In section III we extract the static potentials we need, i.e. the elements of the matrix (1), from lattice QCD data from Ref. [31]. In section IV we discuss numerical methods to solve the coupled channel Schrödiger equation. We use these methods in section V to predict quarkonium bound states and resonances for I = 0 and compare to existing experimental results. Finally, in section VI we conclude and present an outlook.

II. THEORETICAL BASICS OF STUDYING QUARKONIUM RESONANCES USING LATTICE QCD POTENTIALS
In this section we present the theoretical basics of our approach to study quarkonium resonances for isospin I = 0 and various J P C in the Born-Oppenheimer approximation with lattice QCD potentials. We start by analyzing quantum numbers of quarkonium systems and of corresponding decay channels of two heavy-light mesons. Then we derive a coupled channel Schrödinger equation containing quarkonium and two-meson channels. We also discuss, how the potentials appearing in the Schrödinger equation are related to static potentials, which can be computed using lattice QCD. Finally we formulate the boundary conditions for meson-meson scattering and specialize the Schrödinger equation to a specific sector by performing a partial wave decomposition.
Notice that we ignore decays of quarkonium to a lighter quarkonium and a light I = 0 meson, e.g. a σ or an η meson. Such a decay is suppressed by the OZI rule [32][33][34], when compared to the decay to a pair of heavy-light mesons. This is consistent with experimental observations, where the dominant hadronic decay is the decay to a pair of heavy-light mesons. It might be possible to also study these OZI suppressed decays using lattice QCD potentials and the Born-Oppenheimer approximation, but this seems even more technical and difficult than the decays addressed here and, thus, we leave them for future research.
A. Quantum numbers ofQQ (quarkonium) and of M M (two heavy-light mesons) We consider systems with a heavy quark-antiquark pairQQ and either no light quarks (i.e. quarkonium) or another light quark-antiquark pairqq with isospin I = 0 (two heavy-light mesons M =Qq andM =qQ for largē QQ separation). The quantum numbers of these systems are denoted in the following way: • J P C : total angular momentum, parity and charge conjugation of theQQ or theQQqq system.
• S P C Q : spin ofQQ and corresponding parity and charge conjugation.
• J P C : total angular momentum excluding the heavy spins of Q andQ and corresponding parity and charge conjugation (for quarkonium, i.e.QQ without light quarks, J P C coincides with the relative orbital angular momentum L P C of the two heavy quarks).
For the heavy quarks Q andQ we use the following approximations: (1) Heavy quark spins are conserved quantities.
Consequently, the energy levels of theQQ and thē QQqq systems as well as their decays and and resonance parameters do not depend on the spins of the heavy quarks Q andQ, i.e. are independent of S P C Q .
(2) Two of the four components of the Dirac spinors of the heavy quarks Q andQ vanish. Thus, one can write Q = P − Q andQ =QP + , where P ± = (1 ± γ 0 )/2 are the projectors to the large and small components of the non-relativistic limit.
These approximations become exact for static quarks and should still yield reasonably accurate results for b quarks, possibly even for c quarks.
Since J P C as well as S P C Q are conserved, J P C is also conserved. For each value of J P C the corresponding coupled channel Schrödinger equation is different (see section II D 4 for J P C = 0 ++ ). Thus J P C is of central im-portance throughout this work, similar as J P C for systems without heavy quarks, while both J P C and S P C Q are less relevant.

1.QQ (quarkonium bound states and resonances)
As discussed above, for quarkonium, i.e. aQQ pair without light quarks, J P C coincides with the relative orbital angular momentum L P C of the two quarks. Thus, possible values are J P C = 0 ++ , 1 −− , 2 ++ . . .
The coupling of the two heavy spinors can be written according to where Γ Q is a 4 × 4 matrix. It is easy to show, that there are only four linearly independent choices for Γ Q such that the right hand side of Eq. (2) does not vanish: Γ Q = P + γ 5 corresponding to S P C Q = 0 −+ and Γ Q = P + γ j (j = 1, 2, 3) corresponding to S P C Q = 1 −− (see also Table I, where the coupling of two heavy spinors is summarized). Within the approximations discussed at the beginning of section II A quarkonium energy levels are independent of the heavy spins and, thus, independent of Γ Q . QQ (two heavy spinors) Qq andqQ (one heavy and one light spinor) S P C Q and J P C = L P C can be coupled in the usual way to definite total angular momentum J P C . For J = L = 0, 1, 2 all possibilities are listed in Table II. Table II.QQ (quarkonium): possibilities to couple the separately conserved S P C Q and J P C to definite J P C .

2.M M (two heavy-light mesons, the decay channels of quarkonium resonances)
It is convenient to write the spin coupling of the four quarks forming the two heavy-light mesons M andM as (A, B, C and D are spin indices), i.e. such that the heavy spins are coupled with Γ Q and the light spins are coupled with Γ q . To study quarkonium resonances, which can decay into two heavy-light mesons with a coupled channel Schrödinger equation, the quarkonium channel and the two-meson decay channels must have identical S P C Q and J P C (and identical corresponding z components of S Q J P C L P C S P C q Γq type of M andM 0 ++ 0 ++ 0 ++ 1 one P = − and one P = + meson two P = + mesons 2 ++ 0 ++ 1 one P = − and one P = + meson 1 ++ γjγ5 one P = − and one P = + meson 3 −− 1 −− P+γj two P = − mesons P−γj two P = + mesons . . . . . . . . . . . . Table III.M M (two heavy-light mesons): possibilies to couple relative orbital angular momentum L P C and light spin S P C q to given J P C . and J). This implies that the heavy quark spins of M andM have to be coupled with the same Γ Q as in the quarkonium case (2). There are, however, several possibilities to couple the relative orbital angular momentum of the two mesons L P C and the light spin S P C q to given J P C . The algebra is straightforward, when using Table I. For J = 0, 1, 2 all possibilities are listed in Table III together with the corresponding Γ q .
A heavy-light meson M orM with the heavy and the light quark in an S wave can have total angular momen-tum J M = 0 or J M = 1 and parity P = − or P = + . The  negative parity mesons have similar mass and the positive  parity mesons have similar mass However, the positive parity mesons are roughly 400 MeV . . . 500 MeV heavier (cf. e.g. [35,36]). Since including lighter decay channels in a coupled channel Schrödinger equation seems more important than including heavier decay channels, it is a necessary step to analyze, which types of heavy-light mesons correspond to the spin couplings listed in Table III, in particular, which Γ q corresponds to two P = − mesons. This can be done using the Fierz identity, Since the resulting energy levels and resonance parameters will not depend on the heavy spins, we consider from now on exclusively the technically simpler S P C Q = 0 −+ , i.e. Γ A = Γ Q = P + γ 5 . For Γ B = Γ q there are four possibilities as can be seen from Table III: i.e. a linear combination of two negative parity heavy-light mesons, which can be read off by comparing to Table I. Notice that this is the lightest decay channel and, thus, of central importance. Later, when setting up a coupled channel Schrödinger equation, we will include this channel, but neglect the following three channels, which are heavier.
• Γ q = 1: i.e. a linear combination of a negative and a positive parity heavy-light meson.
• Γ q = γ j γ 5 : i.e. a linear combination of a negative and a positive parity heavy-light meson.
• Γ q = P − γ j : i.e. a linear combination of two positive parity heavy-light mesons.

B.QQ andM M coupled channel Schrödinger equation
Setting up the coupled channel Schrödinger equation is independent of the heavy spins. Therefore, as already stated in the previous subsection, we consider the technically simpler S P C Q = 0 −+ , i.e. Γ Q = P + γ 5 . In addition to a quarkonium channel characterized by J P C = L P C , we consider the lightest two-meson decay channel, which contains two negative parity mesons M andM (see Eq. (5)). We ignore the heavier two-meson decay channels, which include one or even two positive parity mesons (see Eqs. (6), (7) and (8)), because their threshold energy is higher by more than 400 MeV or 800 MeV, respectively. The corresponding light spin is S P C q = 1 −− and Γ q = P + γ j (see Table III). Thus, the wave function of the coupled channel Schrödinger equation has four components, ψ(r) = (ψQ Q (r), ψM M (r)). The first component ψQ Q (r) represents theQQ quarko-nium channel, while the remaining three components ψM M (r) correspond to the three spin orientations of the S q = 1 triplet of theM M two-meson channel. r is the relative coordinate of the two heavy quarksQQ, which is for ψM M (r) equivalent to the separation of the two heavy-light mesonsM M .
The spin algebra is [S j , S k ] = i jkl S l . For the three components of ψM M (r) we choose as generators i.e. the generators of rotations around the three Cartesian axes. The eigenvectors of S z are v 0 = (0, 0, 1) with eigenvalue 0 and v ± = (1, ±i, 0)/ √ 2 with eigenvalues ±1. The coupled channel Schrödinger equation reads where are the reduced heavy quark and heavy-light meson masses and L = r × p is the orbital angular momentum operator. The potential V (r) is also a 4 × 4 matrix, which can be written as This particular structure is derived and discussed in the following section II C, where we treat the heavy quarks Q andQ in the static limit and relate V (r) to static potentials from QCD. In section III we explain, how to compute the four functions VQ Q (r), VM M, (r), VM M,⊥ (r) and V mix (r) on the right hand side of Eq. (11) using lattice QCD. The non-zero off-diagonal elements in the first column and the first row proportional to V mix (r) lead to mixing of the quarkonium channel and the two-meson channels and, thus, to quarkonium resonances. Note that in a previous paper [23] we have used similar techniques to derive a coupled channel Schrödinger equation for an I(J P ) = 0(1 + )QQqq tetraquark system, to explore the effect of the heavy quark spins. In this subsection we treat the heavy quarks Q andQ as static quarks. This allows to relate the potential matrix V (r) appearing in the coupled channel Schrödinger equation (10), i.e. the four potentials VQ Q (r), VM M, (r), VM M,⊥ (r) and V mix (r) (see Eq. (11)) to static potentials from QCD, which can be computed using lattice QCD. Moreover, we explain, why V (r) has the particular structure given in Eq. (11).
The positions of static quarks are frozen, for Q and Q w.l.o.g. at +r/2 and −r/2, respectively, i.e. their separation is r = |r|. Thus, rotational symmetry, parity and charge conjugation are broken. Remaining symmetry transformations are rotations around theQQ separation axis, parity combined with charge conjugation (operator P • C) and spatial reflection along an axis perpendicular to theQQ separation axis (corresponding operator denoted as P x ). States are labeled by quantum numbers Λ η , where theQQ spin is not included, because static spins are conserved quantities (for a detailed discussion of this notation, which is also used for homonuclear diatomic molecules and for excited flux tubes, see Refs. [31,[37][38][39]): • Λ = 0, 1, 2, . . . = Σ, Π, ∆, . . . is the absolute value of total angular momentum with respect to theQQ separation axis.
• η = +, − = g, u is the eigenvalue with respect to the operator P • C.
• = +, − is the eigenvalue with respect to the operator P x (sectors with Λ = 0, which differ only in the quantum number , have degenerate spectra; therefore, it is common to list only for Λ = Σ and to omit for Λ = Π, ∆, . . .).
Since the three quantum numbers Λ η do not include thē QQ spin, they play a similar role as J P C , when the heavy quark positions are not frozen. A quarkonium systemQQ with static quarks and the gluonic flux tube in the ground state has quantum numbers Λ η = Σ + g , i.e. the flux tube is invariant under rotations around and reflections parallel and perpendicular to theQQ separation axis. For theM M channels with quark contentQQqq and static quarks Q andQ the situation is more complicated, because of the spins of the two light quarks, which are coupled according toqΓ q q (see Eq. (3)). Since we exclusively consider channels containing two negative parity heavy-light mesons, there are only three independent possibilities, Γ q ∈ {P + γ 1 , P + γ 2 , P + γ 3 } (see Table III). The linear combination Γ q = e r P + γ = (e r ) j P + γ j corresponds to light spin S q = 1 parallel to theQQ separation axis (here and in the following e r , e ϑ and e ϕ denote the orthogonal basis vectors of spherical coordinates). The quantum numbers are Λ η = Σ + g , i.e. are identical to those of the quarkonium system. The remaining two combinations, Γ q = e ϑ P + γ = (e ϑ ) j P + γ j and Γ q = e ϕ P + γ = (e ϕ ) j P + γ j , correspond to light spin S q = 1 perpendicular to theQQ separation axis with quantum numbers Λ η = Π + g and Λ η = Π − g , respectively (for definiteness we have choosen as direction of reflection for P x the ϕ direction; note, however, that final results are independent of this choice).
The computation of static potentials with quantum numbers Λ η with lattice QCD can be done as explained in detail in Ref. [31]. For Λ η = Σ + g there is mixing between a quarkonium-like static potential and a two-meson-like static potential. To compute both potentials, two creation operators are needed, where the first operator is of quarkonium type (U (−r/2; +r/2) is a straight spatial parallel transporter connecting −r/2 and +r/2, typically a product of smeared spatial links), while the second operator is of two-meson type. From a Fig. 1 and C M (t) denotes the correlation function of the J P = 0 − static-light meson (see e.g. Refs. [36,40]), one can extract the energy eigenvalues V Σ + g 0 (r) and V Σ + g 1 (r). These eigenvalues correspond to the two lowest energy eigenstates in the Σ + g sector, |0; Σ + g and |1; Σ + g , and are normalized with respect to 2m M . The Π + g sector as well as the Π − g sector do not include quarkonium states. Thus, for each of them a single creation operator is sufficient, differing from the operator (13) only in Γ q , From the corresponding normalized correlation functions one can extract the energy eigenvalues V Π + g 0 (r) and V Π − g 0 (r) of the lowest energy eigenstate in each of the sectors, |0; Π + g and |0; Π − g . Note that the spectra in these sectors are The previously defined energy eigenvalues and eigenstates fulfill where H denotes the Hamiltonian.
If the gluonic parallel transporter U and the light quark fields u and d appearing in the creation operators (12) and (13) where θ(r) is the mixing angle (for details see Ref. [31], section 5.A). For separations r somewhat below the string breaking distance r sb ≈ 1.1 fm the lowest energy eigenstate is predominantly a quarkonium state, while the first excitation is a two-meson state, i.e. θ(r) ≈ 0.3 . . . 0.4. For r somewhat above r sb the situation is reversed, i.e. θ(r) ≈ π/2. For a detailed discussion, of how to compute the mixing angle from the 2 × 2 correlation matrix see Ref. [31]. One can obtain the analog of Eq. (16) for a basis including |QQ and |M M instead of |0; Σ + g and |1; Σ + g by using Eq. (17), where we have defined |M M ± ⊥ = |0; Π ± g and V mix (r) = cos(θ(r)) sin(θ(r)) V To express the 4 × 4 potential matrix V (r) appearing in the Schrödinger equation (10) in terms of static potentials in QCD, recall that the first row and column of this matrix corresponds to aQQ quarkonium channel, while the remaining three rows and columns correspond to the spin-1 triplet of a two-mesonM M channel. Thus, where |M M j , j = 1, 2, 3 are linear combinations of two-meson states |M M , |M M + ⊥ and |M M − ⊥ , which have light spin in j direction. These linear combinations are given by Combining Eqs. (18), (23) and (24) leads to which is identical to Eq. (11). Thereby, we have derived the structure of the 4 × 4 potential matrix appearing in the Schrödinger Eq. (10). Moreover, via Eqs. (19) to (20) we have related the four potentials VQ Q (r), VM M, (r), V mix (r) and VM M,⊥ (r) to the static potentials V Σ + g 0 (r), V Σ + g 1 (r) and V Πg 0 (r) and the mixing angle θ(r), which can be computed using lattice QCD.
D. Boundary conditions for MM scattering and partial wave decomposition Now we consider specific boundary conditions appropriate to describe scattering of two heavy-light mesons M andM , an incident plane wave and an emergent spherical wave for largeQQ separations r (for previous work on similar systems see Refs. [2,28]). Moreover, we do a partial wave decomposition and specialize the Schrödinger equation (10) to definite J P C = 0 ++ , i.e. vanishing total angular momentum excluding the heavy spins. In other words, we formulate the Schrödinger equation (10) specifically for quarkonium bound states and resonances with J P C = L P C = 0 ++ (for the relation of J P C to the common J P C as e.g. used by the Particle Data Group [41] see Table II). This reduces the partial differential equation (10) to a system of two coupled ordinary differential equations in the radial coordinate r = |r|, which is much simpler to solve numerically.

Basis functions for the partial wave decomposition
In standard textbooks on quantum mechanics scattering theory is typically discussed for a spin-0 system and a single channel, i.e. the corresponding Schrödinger equation has only one component and the partial wave decomposition is done in terms of spherical harmonics, which are eigenfunctions of L 2 and L z . For our particular problem an equivalent decomposition is technically more complicated, because the Schrödinger equation (10) has two channels, aQQ quarkonium channel with spin 0 (upper component of the wave function) and aM M twomeson channel with spin 1 (lower three components of the wave function), i.e. ψ(r) = (ψQ Q (r), ψM M (r)). The partial wave decomposition has to be done in terms of four-component eigenfunctions of J 2 and J z , which are conserved quantities and replace the non-conserved L 2 and L z . An orthonormal and complete set of eigenfunctions with respect to to the solid angle Ω and the four components of ψ(r) is the following: i.e. proportional to spherical harmonics.
• Eigenfunctions of J 2 and J z , which are non-zero in the lower three components, can be constructed by Clebsch-Gordan coupling of spherical harmonics Y L,Lz (Ω) (representing relative orbital angular momentum L) and the three spin-1 components of the light quarks. They are given by where L = 1 for J = 0 and J − 1 ≤ L ≤ J + 1 for J ≥ 1. Z L→ J, Jz (Ω) is defined as follows: (j = x, y, z replaces the index J z = −1, 0, +1 in the usual way, i.e. z ≡ 0 and ∓(x ± iy) ≡ ±1).
-J ≥ 2: The corresponding functions Z L→ J, Jz (Ω) can be constructed in a straightforward way. Since these functions are not needed explicitly in this work, we do not provide equations.
Any three-component function G(r) can be written as an expansion in Z L→ J, Jz (Ω), where the expansion coefficients g L→ J, Jz (r) are functions of r = |r|. Note that all Z L→ J, Jz (Ω) have parity (−1) L+1 . Thus for parity even functions G(r) all coefficients g L→ J, Jz (r) with even L are zero, while for parity odd functions G(r) all coefficients g L→ J, Jz (r) with odd L are zero. In section III, where we determine the potentials VQ Q (r), VM M, (r), VM M,⊥ (r) and V mix (r) by parameterizing lattice QCD results from Ref. [31], we find VM M, (r) → 0, VM M,⊥ (r) → 0 and V mix (r) → 0 for QQ separations r → ∞ (see Eqs. (61), (59) and (60)). This is expected, because the potentials are normalized by an additive constant in such a way that a value of 0 corresponds to theM M threshold. Consequently, the 4 × 4 potential matrix (11) reduces to i.e. the Schrödinger Eq. (11) decouples into four independent partial differential equations, one for each of the four components of the wave function ψ(r) = (ψQ Q (r), ψM M (r)).
The potential in the quarkonium equation, VQ Q (r), is linear for r → ∞, i.e. confining (see Eq. (58)). Thus the boundary conditions for ψQ Q (r) is The where t L→ J, J /k are the scattering amplitudes. We have included the prefactor − √ 4π in front of t 1→0,0 , because using probability conservation one can show This in turn allows to define the corresponding scattering phase δ 1→0,0 as which closely resembles textbook conventions for standard scattering of spinless particles. Eq. (36) is equivalent to which is the optical theorem.

Partial wave decomposition
The partial wave decomposition of ψQ Q (r) is an ordinary expansion in spherical harmonics, with functions u J, Jz (r) as coefficients (see Eq. (26)). It is convenient to include the prefactor √ 4πiA z in front of u 0,0 (r), to avoid unnecessary factors in the coupled channel Schrödinger equation for J P C = 0 ++ , which we will derive in section II D 4. The boundary conditions are where the latter follows from Eq. (34). We write ψM M (r) as a sum of the incident wave Ae +ikz and an emergent wave X(r), ψM M (r) = Ae +ikz + X(r).
The partial wave decomposition of ψM M (r) follows Eq.
with functions a L→ J, Jz (r) and χ L→ J, Jz (r) as coefficients.
Again it is convenient to include the prefactor √ 4πiA z in front of χ 1→0,0 (r). For the incident wave the coefficients a L→ J, Jz (r) can be calculated in a straightforward way, e.g. a 1→0,0 (r) = dΩ (Z 1→0,0 (Ω)) * Ae ikz = √ 4πiA z j 1 (kr) (j 1 denotes a spherical Bessel function of the first kind), which is particularly relevant in the following. The boundary conditions are These two equations can also be expressed in matrix form, Eqs. (48) and (49) or equivalently Eq. (50) are corner stone equations of this work. In section V we will solve them to obtain numerical results for bottomonium bound states and resonances. Specializing the coupled channel Schrödinger equation (10) to J P C = 1 −− or higher J can be done in the same way. We leave that for future publications.

III. UTILIZING THE LATTICE QCD STATIC
POTENTIALS FROM REF. [31] To determine the potentials VQ Q (r), VM M, (r), VM M,⊥ (r) and V mix (r), one has to compute the following correlation functions with lattice QCD, as discussed in section II C: • A normalized 2 × 2 correlation matrix using the creation operators (12) and (13).
• A normalized correlation function using either the creation operator (14) or the creation operator (15).
Such computations are quite challenging and technicalities are discussed in detail e.g. in [31,42]. In the future we plan to perform such computations. Here we follow a different strategy and reuse the existing lattice QCD results for static potentials from Ref. [31], to determine VQ Q (r), VM M, (r), VM M,⊥ (r) and V mix (r) within certain approximations. In Ref.
[31] a 2 × 2 correlation matrix was computed at light u and d quark mass corresponding to a pion mass m π ≈ 654 MeV and lattice spacing a ≈ 1/(2.37 GeV) ≈ 0.083 fm using the creation operators O [31] QQ = Q (−r/2)e r P + γU (−r/2; +r/2)Q(+r/2) O [31] M M = Q (−r/2)P + γ 5 u(−r/2) ū(+r/2)γ 5 P − Q(+r/2) + (u → d) (see Eqs. (11) and (15) in Ref. [31]; for convenience we have expressed these operators in the same notation used in previous sections of this work and we have inserted projectors P + using Q = P − Q andQ =QP + ). Using the Fierz identity (4) these operators can be rewritten according to Four of the five operators appearing on the right hand sides are defined in Eqs. (12) to (15), where Γ Q has to be replaced as indicated in square brackets. O Note that the two-meson creation operator O [31] M M used in Ref. [31] does not only probe the Σ + g sector, i.e. the sector of the ordinary static potential, but also the three sectors Σ − u , Π + g and Π − g . To parameterize the three independent elements of the normalized 2 × 2 correlation matrix for large temporal separations t, we follow the arguments of Ref. [31], It is assumed that for large t two energy eigenstates are sufficient to describe the contributions from the Σ + g sector (which are a mixtures of aQQ quarkonium state and aM M two-meson state), while for each of the sectors Σ − u , Π + g and Π − g only a single energy eigenstate is needed (aM M two-meson state). V  [31], where it is now taken into account that O [31] M M probes several Λ η sectors. The correlation matrix data of Ref. [31] is not publicly available, but a parametrization of this data is given (see. Eqs. (68) to (70) and Table I in Ref. [31]). This parametrization allows us to resample the correlators from Ref. [31] and to perform χ 2 minimizing fits of our parametrization (55) to (57) to the resampled data. For stable fits the number of parameters is, however, too large. Thus, we assume V Inserting the fit results for V Σ + g 0 (r), V Σ + g 1 (r) and θ(r) into Eqs. (19) to (21) leads to the potentials VQ Q (r), VM M, (r) and V mix (r), which appear in the coupled channel Schrödinger equation (Eqs. (10) and (11) or specifically for J P C = 0 ++ Eq. (50)). These potentials are also collected in Table IV, columns 7  represent the interaction in the quarkonium and the twomeson channel and the mixing between these channels: • V Σ + g 0 (r) for separations r < ∼ 1 fm is a standard quantity computed in lattice QCD and commonly referred to as "the static potential". It has a negative curvature for small r and is almost linear for larger r. Quite often it is parameterized via V Σ + g 0 (r) = const−α/r+σr (see e.g. Ref. [43]). With this parameterization it is straightforward to determine the common hadronic scale r 0 defined via [44]. Fitting the parameterization in the range 3a ≤ r ≤ 12a yields r 0 = 6.004(41) a, which is in excellent agreement with the original result quoted in Ref. [31]), r 0 = 6.009(53) a.
• VQ Q (r) is the potential between a static quark and a static antiquark. In contrast to V Σ + g 0 (r) this potential is linear also for separations r > ∼ 1 fm, i.e. separations larger than the string breaking distance. At smaller r there is a sizable r-dependent mixing of the two lowest energy eigenstates (see Fig. 2, lower plot), and VQ Q (r) is an r-dependent linear combi-  (19)). Thus one should not expect a simple parameterization similar to the previously mentioned const − α/r + σr. The mixing manifests itself by a clearly visible bump around r ≈ 3a and needs to be taken into account, when parameterizing VQ Q (r).
• VM M, (r) is the potential between two static-light mesons. For larger r the residual strong force between the two mesons is expected to vanish and the potential should approach two times the static-light meson mass. This expectation is in excellent agreement to what we observe for r > ∼ 8a. For smaller r the statistical errors are larger and the potential might either be constant or slightly repulsive. Note that V Σ + g 1 (r) is somewhat larger than VM M, (r). Again, this is a consequence of the non-vanishing mixing angle.
To solve the Schrödinger equation (10), which we do in the following for J P C = 0 ++ (see Eq. (50)), it is necessary to have continuous parameterizations of the potentials VQ Q (r), VM M, (r) and V mix (r). We tested several parameterizations and found the following most suitable: with 11 parameters, E 0 , α, σ, cQ Q,j , λQ Q,j , c mix,j and λ mix,j (where j = 1, 2). These parameters can be determined by χ 2 minimizing fits in a stable way and the corresponding χ 2 /dof indicate reasonable fits. Results are collected in Table V. and the fitted parameterizations are shown in Fig. 3. Moreover, Eq. (22) and our above assumption V Πg 0 (r) = 0 lead to VM M,⊥ (r) = 0. (61)

A. S and T matrix poles in the complex energy plane and their relation to quarkonium resonances
The quantity t 1→0,0 appearing in the r → ∞ boundary condition (47) of the radial wave function χ 1→0,0 (r) of the coupled channel Schrödinger equation (50) is an eigenvalue of the T matrix. From t 1→0,0 we can read off the corresponding S matrix eigenvalue, s 1→0,0 = 1 + 2it 1→0,0 = e 2iδ1→0,0 .
Moreover, both the S matrix and the T matrix are analytical in the complex plane. They are well-defined for potential parameter in units of a in units of GeV  complex energies E. The poles of the S and the T matrix, i.e. the poles of t L→ J, Jz , which are in the second Riemann sheet of the complex energy plane with a negative imaginary part, correspond to quarkonium resonances. For a pole at complex energy E the resonance energy and the decay width are For more details see e.g. our recent work [28].
B. Numerical methods to determine t1→0,0 and to find poles Since we restrict our numerical calculations in this work to the sector J P C = 0 ++ , i.e. the Schrödinger equation (50), it is convenient to use the simplified notation u(r) ≡ u 0,0 (r) and χ(r) ≡ χ 1→0,0 (r).
The boundary conditions of the solutions (u(r), χ(r)) of the Schrödinger equation (50) can be read of from Eqs. (40), (41), (46) and (47) and are Note that the boundary condition (67) depends on t 1→0,0 . For a given value of the energy E this boundary condition is only fulfilled for a specific corresponding value of t 1→0,0 . In other words the boundary condition (67) fixes t 1→0,0 as a function of E. Thus, our numerical goals in the following are to compute t 1→0,0 for given values of the energy E and to find the poles of t 1→0,0 in the complex energy plane.
The coupled channel Schrödinger equation (50) is then a system of 2(N − 1) linear equations, for the 2(N − 1) unknowns u 1 , . . . , u N −1 and χ 1 , . . . , χ N −1 representing the radial wave functions. Inserting the boundary conditions leads to To determine t 1→0,0 for a given value of the energy E is now straightforward. We solve the linear system (71) and insert χ N −1 in Eq. (68).
In case we are just interested to find the positions of poles of t 1→0,0 in the complex energy plane, we use a more efficient method. From x = M −1 (E)b, where all components of b are finite, it is obvious that M −1 (E) must have at least one infinite eigenvalue, if |t 1→0,0 | → ∞ or equivalently |χ N −1 | → ∞. Thus, M (E) must have at least one zero mode, which implies det(M (E)) = 0. Thus, we determine poles of t 1→0,0 by finding the roots of det(M (E)). To this end we apply the Newton-Raphson method for analytic functions of a single complex variable. Note that det(M (E)) is typically a large number beyond machine double precision. Thus, we rescale M (E) by an appropriate factor, after a LU decomposition, but before multiplying the diagonal elements to obtain the determinant. The roots of the determinant are, of course independent, of such a rescaling, i.e. the final results for the positions of the poles of t 1→0,0 in the complex energy plane are unaffected.
Instead of rewriting the coupled channel Schrödinger equation (50) as a large system of linear equations (71), one can also solve it and determine t 1→0,0 using Runge-Kutta methods. We cross checked and verified our numerical results by implementing a 4th order Runge-Kutta solver.

A. Choice of parameters and error analysis
In the following we focus on b quarks and bottomonium, where the heavy quark approximations discussed in section II A should be justified. For m M , which is the energy reference of our system, we use the spinaveraged mass of the B meson and the B * meson, i.e. m M = (m B + 3m B * )/4 = 5.313 GeV [41]. µ Q = m Q /2 in the kinetic term of the coupled channel Schrödinger equation is the reduced mass of the b quark. In the quark model perspective the B meson is composed of a b quark/antiquark and a light antiquark/quark u or d. Thus the B meson mass is heavier than the heavy quark mass, where the difference is of the order of the light constituent quark mass m l , i.e. m M = m Q + m l . Since results are only weakly dependent on m Q (see e.g. previous work following a similar approach [9,43]), we use for simplicity m Q = 4.977 GeV from quark models [45].
For the upper boundary of the r axis we use R = 15.0/GeV ≈ 2.96 fm ≈ 35.6 a, where VQ Q (r) is quite large and both VM M, (r) and V mix (r) are essentially vanishing (see e.g. Figure 3). For the 1-dimensional lattice discretizing the interval [0, R] we use N = 600 sites corresponding to the spacing d = R/N = 0.025/GeV ≈ 0.005 fm. We verified the independence of our results from these parameters for R ≥ 15.0/GeV and N ≥ 600 by performing identical computations with several different R and N .
We propagate the uncertainties provided in Ref. [31], TABLE I by resampling. We generate 1000 statistically independent samples and repeat all computations on each of the samples. This is a computer time consuming task, for which we employ GPUs utilizing the algebra package of CUDA [46] and the cuSOLVER library [47]. For t 1→0,0 , the corresponding phase shift δ 1→0,0 , energy lev-els and positions of poles of t 1→0,0 in the complex energy plane we quote asymmetric errors, which are defined via the 16th and 84th percentile of the 1000 samples, respectively.

B. The bottomonium spectrum from single channel Schrödinger equations
In this subsection we are interested in a qualitative understanding, how results from single channel Schrödinger equations compare to experimental results. Thus we just use the mean values of the parameters of the potentials VQ Q (r), VM M, (r) and V mix (r) (Eqs. (58) to (60)), but ignore their uncertainties. It will also be interesting to compare these single channel results to more realistic results from the coupled channel Schrödinger equation (50), which will be discussed in sections V C and V D.
In a first step we compute the quarkonium spectrum, setting V mix (r) = 0, i.e. decoupling the quarkonium channel from the meson-meson channel. The corresponding Schrödinger equation is then the upper component of Eq.
There is an infinite number of bound states, because the potential VQ Q (r) is confining. The energies of the lightest states are listed in Table VI. Three states are below the B ( * )B( * ) threshold at 2m M , which is marked by a horizontal line, all other states are above. In particular the mass of the lowest state (n = 1) is significantly larger than the corresponding masses from experiment (compare columns "from VQ Q (r)" and "from experiment" of Table VI; for a full summary of the experimental results on bottomonium see Table VII). This sizable discrepancy is expected, because the mixing angle is non-zero, where the wave function is large, i.e. 0.3 < ∼ θ(r) < ∼ 0.5 for separations r < ∼ 1.1 fm (see Fig. 2). Consequently, the potential VQ Q (r) is a mixture of the ground state potential and the first exited potential, which leads to unphysically large bottomonium masses.
More realistic estimates for these masses are obtained, when replacing VQ Q (r) in Eq. (72) by the ground state potential V Σ + g 0 (r), i.e. by solving the Schrödinger equation This is a standard approach appearing frequently in the literature (for recent work see e.g. [43]). The resulting bottomonium masses are now smaller and compare much better to experimental results (see Table VI). From Table VII one can also read off that the decay widths of bottomonium states are indeed quite small below the BB threshold, except for η b (1S), which decays electromagnetically via the U (1) anomaly. Thus, our approach to neglect the OZI suppressed decays of excited  name  Table VII. Bottomonium states with isopspin I = 0 according to the Review of Particle Physics [41]. We also list the quantum numbers J P C conserved in the limit of infinite b quark mass ( J = 0, 1, 2 corresponds to S, P, D in the meson name; the quantum numbers J P C = 1 −− of Υ(10860) and Υ(11020) are consistent with J P C = 0 ++ , which is the sector we focus on in this work). B ( * )B( * ) thresholds are marked by horizontal lines.
bottomonium states to lighter bottomonium and a light I = 0 meson should be a good approximation.
C. t1→0,0 and the phase shift δ1→0,0 for real energies We proceed by computing the scattering amplitude t 1→0,0 and the phase shift δ 1→0,0 for real energies E above the B ( * )B( * ) threshold at 10.627 GeV. In contrast to the previous subsection, we now include the meson-meson channel, i.e. we consider the coupled channel Schrödinger equation (50). The available experimental data for bottomonium goes up to ≈ 11 GeV (see Table VII). Two resonances consistent with J P C = 0 ++ were already observed, Υ(10860) and Υ(11020). We perform our computations up to 11.6 GeV.
In Fig. 4 we show the real and the imaginary part of t 1→0,0 as functions of the energy E. If there would be "simple resonances", Im(t 1→0,0 ) would exhibit clear peaks, top to bottom, since according to Eq. (38) it is identical to |t 1→0,0 | 2 , i.e. proportional to the absolute square of the partial wave scattering amplitude. In such a case we could determine the decay widths from the peaks at half height. However, the system we investigate is more complicated, e.g. the resonances seem to mutually impact each other and there is a large background.  Notice Eq. (38) is equivalent to the equation of circle in the complex plane centered at i/2 with radius 1/2. In Fig. 5 we show the corresponding Argand diagram. We get the expected circle, which is a good test that we are complying with the probability conservation expressed in the optical theorem as well as of the numerical precision of our results. In Fig. 6 we show the phase shift δ 1→0,0 , which can be obtained from t 1→0,0 via Eq. (37), as function of the energy. In this plot the resonances can be identified more clearly, since each of them corresponds to a "jump" of the order of π. The steepness of each jump is inversely proportional to the the corresponding decay widths. There are three clear resonances close to 11.1 GeV, 11.3 GeV and 11.5 GeV. Moreover, there seems to be another wider and less clear resonance around 10.85 GeV. To determine resonance energies and decay widths precisely, we consider the analytic continuation of our scattering problem to the complex energy plane. There we search for the poles of t 1→0,0 using the Newton-Raphson method as discussed in section V A. The positions of the poles E are related to the resonance masses and decay widths according to Eq. (63), i.e. m = Re(E) and Γ = −2Im(E).
In Fig. 7 we show the positions of the poles of t 1→0,0 in the complex energy plane for all bound states and resonances below 11.6 GeV. For each bound state and each resonance there is a differently colored point cloud representing the 1000 resampled sets of parameters of the potentials, which we use to determine statistical errors (see section V A for details). For the bound states the poles are located on the real axis below the B ( * )B( * ) threshold. For the resonances the positions of the poles follow curved bands, where the imaginary parts range from almost vanishing values to finite values comparable to those observed in experiments (see Table VII), i.e. of the order of tens of MeV. There are clear gaps between the point clouds representing different bottomonium bound states and resonances, which allows a straightforward error analysis. The corresponding mean values and errors are indicated by the black circles and crosses. These results are also summarized in Table VIII We also find several resonances (i.e. n ≥ 5), where the majority has not yet been observed experimentally. Our resonance mass for n = 6 is quite similar to the experimental result for Υ(10860), which might indicate that Υ(10860) should be interpreted as Υ(5S). For Υ(11020), on the other hand, there is no perfect match among our theoretical results. The closest resonance we find (n = 7) is almost 100 MeV heavier. An explanation could be that Υ(11020) is not a state with J P C = 0 ++ , but with higher J. It will be interesting to explore this further in the future, by deriving and solving coupled channel Schrödiger equations also for J ≥ 1.
Note that we also predict a resonance (n = 5) below the resonance masses of the two experimental J P C = 0 ++ candidates Υ(10860) and Υ(11020), not far away from the B ( * )B( * ) threshold. This resonance and the resonance with n = 6, which is a candidate for Υ(10860), with masses close to 10.8 GeV and 10.9 GeV, respectively, are illustrated in Fig. 8, which is a 3D plot of the absolute value and the phase of t 1→0,0 in the complex energy plane. As expected, the phase performs a full 2π revolution around each of the corresponding poles. Note that a clear identification and separation of those two resonances is only possible from a pole analysis in the complex energy plane, but not from Re(t 1→0,0 ), Im(t 1→0,0 ) and δ 1→0,0 at real energies (see Fig. 4 and Fig. 6). For example, as discussed at the end of section V C, Fig. 6 suggests that there might be a wide resonance around 10.85 GeV, but it is impossible to see from that figure that there are actually two resonances in that energy re-gion.
A closer inspection reveals that the lower of the two resonances (n = 5) has a fully dynamical origin. This can be seen from Fig. 9, where we study both Im(t 1→0,0 ) and δ 1→0,0 for mixing potentials cV mix (r) with c ∈ {1/ √ 10, 1/ √ 2, 1}. In particular in the left plot, where we show Im(t 1→0,0 ) as a function of real energy E, one can see that the peak at around 10.9 GeV becomes more pronounced for decreasing c. At the same time the very wide peak below 10.8 GeV also transforms into a sharp and clear peak and moves to significantly smaller energies slightly above the B ( * )B( * ) threshold. At small mixing, c = 1/ √ 10, however, there are only three bound states, not four as for c = 1. In other words, when decreasing the mixing potential, the bound state with n = 4 becomes a clear resonance close to the threshold, while the wide resonance with n = 5 disappears.
In what concerns the imaginary parts of the positions of the poles, we obtain below the B ( * )B( * ) threshold exactly zero as expected. Above the B ( * )B( * ) threshold only the lightest resonances (n = 5 and n = 6) have decay widths similar to the experimental J P C = 0 ++ candidates Υ(10860) and Υ(11020). All higher resonances (n ≥ 7) have significantly smaller widths. The reason for this is probably that we consider only the coupling of quarkonium to the lightest meson-meson channel. To obtain larger widths, we would need to include all excited meson-meson channels up to the respective resonance masses. Clearly, this goes beyond the scope of the present work.
We stress that the errors on our theoretical results quoted in Table VIII and shown in the figures are purely statistical. These results, however, were obtained by resorting to certain approximations. First of all, the coupled channel Schrödinger equation (10) was derived in the static limit (see the discussion in section II A), while b quarks have a large, but finite mass. Similarly, we use lattice QCD potentials V Σ + g 0 (r) and V Σ + g 1 (r) and a mixing angle θ(r), which were computed in the static limit and at unphysically heavy u and d quark mass corresponding to m π ≈ 654 MeV and a single lattice spacing a ≈ 1/(2.37 GeV) ≈ 0.083 fm. Moreover, we assumed V Πg 0 (r) = 0. As already mentioned, our coupled channel Schrödinger equation contains only the lightest decay channel to two negative parity heavy-light mesons, while the next ones containing a negative and a positive parity meson are around 400 MeV . . . 500 MeV above (see Eqs. (6) and (7)). Thus we expect that the neglect of this channel has little effect on our results up to the corresponding threshold around 11.025 GeV . . . 11.125 GeV, i.e. for n ≤ 6, while higher resonances with n ≥ 7 might be strongly affected. Finally we separated the treatment of heavy and light degrees of freedom using the Born-Oppenheimer approximation. For the obtained masses we crudely estimate systematic errors to be around 50 MeV, which is the maximum discrepancy of experimental results and our theoretical predictions for  Table VIII. Masses and decay widths for J P C = 0 ++ bottomonium from the coupled channel Schrödinger equation (50) (column "from poles of t1→0,0"). For comparison we also list corresponding single channel results (column "from V Σ + g 0 (r)") and experimental results (column "from experiment"). Relevant B ( * )B( * ) thresholds are marked by horizontal lines. Errors on our theoretical results are purely statistical. the three bound states Υ(1S), Υ(2S) and Υ(3S).

VI. CONCLUSIONS
We proposed and derived a formalism to study quarkonium bound states and resonances with I = 0 based on static potentials from QCD, which can be computed with lattice QCD. We applied the Born-Oppenheimer approximation by inserting these potentials in a specifically derived coupled channel Schrödinger equation for the dynamics of heavy quarks. This equation, which contains a quarkonium and a heavy-light meson meson channel, allows to predict masses of bound states and resonances as well as decay widths. For the resonances we apply scattering theory, which enables to compute phase shifts arg(t1→0,0) and eigenvalues of the T matrix in the complex energy plane.
Within our framework we studied bottomonium states with I = 0 up to 11.6 GeV focusing on the J P C = 0 ++ channel, which corresponds to L = 0 for thebb pair in the quarkonium channel and L = 1 for theB ( * ) B ( * ) pair in the meson-meson channel. Even though we resort to several approximations (see the discussion at the end of section V D), we find reasonable agreement with the experimentally observed bottomonium spectrum. There are four bound states, which can clearly be identified with η b ≡ Υ(1S), Υ(2S), Υ(3S), Υ(4S). We also obtain a resonance around 10.870 GeV, which matches Υ(10860) rather well, suggesting that Υ(10860) could be interpreted as Υ(5S). For Υ(11020), on the other hand, we do not find a close-by resonance, which might be an indication that Υ(11020) is not an S wave state. Moreover, we predict a new, dynamically generated resonance close the theB ( * ) B ( * ) threshold with mass ≈ 10.790 GeV and decay width ≈ 51 MeV.
A straightforward next step will be to study bottomonium with J ≥ 1. The corresponding coupled channel Schrödiger equations will have at least a 3 × 3 matrix structure. For instance J P C = 1 −− corresponds to L = 1 for thebb pair in the quarkonium channel and L = 0 or L = 2 for theB ( * ) B ( * ) pair in the meson-meson channels. Thus, one can study a possibly existing X b meson, the counterpart of the famous X c (3872) [48].
Another direction for the future could be to include the decay channels to a negative and a positive parity heavy-light meson. This would allow to make more realistic predictions for resonances with n ≥ 7 up to the threshold of two positive parity heavy-light mesons at around 11.525 GeV. The corresponding static potentials, however, have not yet been computed with lattice QCD. Moreover, our current determination of the potentials VQ Q (r), VM M, (r), VM M,⊥ (r) and V mix (r) from the lattice QCD results of Ref. [31] requires certain assumptions. Thus, we plan to perform a dedicated lattice QCD computation of all those static potentials, possibly also with u and d quark mass closer to the physical value.
Finally it would be worthwhile to include the effects of the heavy spins, either on the level of the coupled channel Schrödinger equation as in Ref. [23] or even in a direct way, when computing the static potentials with lattice QCD (see e.g. Refs. [49][50][51]).