Zb tetraquark channel from lattice QCD and Born-Oppenheimer approximation

Two $Z_b$ hadrons with exotic quark structure $\bar bb\bar du$ were discovered by Belle experiment. We present a lattice QCD study of the $\bar bb\bar du$ system in the approximation of static b quarks, where the total spin of heavy quarks is fixed to one. The energies of eigenstates are determined as a function of separation $r$ between $b$ and $\bar b$. The lower eigenstates are related to a bottomonium and a pion. The eigenstate dominated by $B\bar B^*$ has energy significantly below $m_B+m_{B^*}$, which points to sizable attraction for small $r$. The attractive potential $V(r)$ between $B$ and $\bar B^*$ is extracted assuming that this eigenstate is related exclusively to $B\bar B^*$. The Schr\"odinger equation for $B\bar B^*$ within the extracted potential leads to a virtual bound state, whose mass depends on the parametrization of the lattice potential. For certain parametrizations, we find a virtual bound state slightly below $B\bar B^*$ threshold and a narrow peak in the $B\bar B^*$ rate above threshold - these features could be related to $Z_b(10610)$ in experiment. We surprisingly find also a deep bound state within undertaken approximations.

Two Z b hadrons with exotic quark structurebbdu were discovered by Belle experiment. We present a lattice QCD study of thebbdu system in the approximation of static b quarks, where the total spin of heavy quarks is fixed to one. The energies of eigenstates are determined as a function of separation r between b andb. The lower eigenstates are related to a bottomonium and a pion. The eigenstate dominated by BB * has energy significantly below mB + mB * , which points to sizable attraction for small r. The attractive potential V (r) between B andB * is extracted assuming that this eigenstate is related exclusively to BB * . The Schrödinger equation for BB * within the extracted potential leads to a virtual bound state, whose mass depends on the parametrization of the lattice potential. For certain parametrizations, we find a virtual bound state slightly below BB * threshold and a narrow peak in the BB * rate above threshold -these features could be related to Z b (10610) in experiment. We surprisingly find also a deep bound state within undertaken approximations.
We explore this channel within first-principle lattice QCD. The only preliminary lattice study of this channel has been reported in [18,19] and is reviewed below. No other lattice results are available since this channel presents a severe challenge. Scattering matrix would have to be determined using the Lüscher method for at least 7 coupled channels listed in the previous paragraph. Poles of scattering matrix would render possible Z b states. Following this path seems too challenging at present.
In the present study we consider the Born-Oppenheimer approximation [20], inspired by the study of this system in [18,19]. It is applied in molecular physics since ions are much heavier than other degrees of freedom. It is valuable also for the Z b systembbdu, where b andb represent heavy degrees of freedom (h), while light quarks and gluons are light degrees of freedom (l), see for example [21,22]. The simplification comes from the fact that heavy degrees of freedom have large mass and therefore small velocity and kinetic energy.
In the first step we treat b andb as static at fixed distance r (Figure 1a) and the main purpose is to determine eigen-energies E n (r) of this system. This energy represents the total energy without the kinetic and rest energies of the b andb, so E n (r) is related to the potential V (r) felt by the heavy degrees of freedom. In the second step, we study the motion of the heavy degrees of free-dom (with physical masses) under the influence of the extracted potential V (r). Solutions of the Schrödinger equation render information on possible (virtual) bound states Z b , resonances and cross-sections.
The low-lying eigenstates of the system in Fig. 1a and quantum numbers (2) are related to two-hadron states in (1) The eigen-energy E n (r) related to BB * in Fig. 1b is of most interest since Z b lie near BB * threshold. Bottomonium-pion states Υ(r)π( p) represent ground state at small r. Here Υ(r) denotes spin-one bottomonium whereb and b are separated by r. Pion can have zero or non-zero momentum p = n 2π L since total momenta of light degrees of freedom is not conserved in the presence of static quarks, i.e. pion momentum can change when it scatters on infinitely heavy Υ. Our task is to extract energies of all these eigenstates E n (r) as a function of separation r. The only previous lattice study of this system [18] presents preliminary results based on Fock components BB * and Υπ(0); the presence Υπ( p = 0) was mentioned in [19], but not included in the simulation.
Quantum numbers and operators: We consider Z 0 b that has quantum numbers I = 1, I 3 = 0, J P C = 1 +− and J z = 0 in experiment. The list of conserved quantum numbers is slightly different in the systems with two static particles. We study the system in Fig. 1a with quantum numbers where the neutral system is considered so that Cconjugation can be applied ( Fig. 1  ) is conserved. The quantum number corresponds to the reflection over the yz plane. P refers to inversion with respect to mid-point between b andb and C is charge conjugation, where only their product is conserved.
The spin of the infinitely heavy quark can not flip via interaction with gluons, so spin S h ofbb is conserved. We choose to study the system with S h = 1, which can decay to Υ, while it can not decay to η b and h b since these carry S h = 0. Note that the physical Z b and BB * with finite m b can be a linear combination of S h = 1 as well as S h = 0, and we study only S h = 1 component here. We have in mind this component, which includes BB * ,BB * , B * B * (O 1 in Eq. 3), when we refer to "BB * " throughout this paper.
The eigen-energies E n of the system in Fig. 1a are determined from the correlation functions O i (t)O † j (0) . We employ 6 operators O i that create/annihilate the system with quantum numbers (2) and resemble Fock components (1) in Figs. 1 (b-d) Here , momenta is given in units of 2π/L, capital (small) letters represent Dirac (color) indices, color singlets are denoted by [..] and U is a product of gauge links between 0 and r. First line in O 1 decouples spin indices of light and heavy quarks, so that J l z and S h (z) (2) are more transparent [18], while the second line is obtained via Fierz transformation. O 2 is obtained from O 1 by replacing all q(x) with ∇ 2 q(x). O 4,5 have pion momenta in z direction due to J l z = 0 and have two terms to ensure C · P = −1. The Υb 1 is not a decay mode for finite m b where C and P are separately conserved, but it is has quantum numbers (2) for m b → ∞. The pairqq indicates combination uu −dd with I = 1 and I 3 = 0. All light quarks q(x) are smeared around position x using full distillation [23] with radius about 0.3 fm, while heavy quarks are point-like.
We verified there are no other two-hadron states in addition to (1) with quantum numbers (2) and with noninteracting energies (4) Lattice details: Simulation is performed on an ensemble with dynamical Wilson-clover u/d quarks, m π 266(5) MeV, a 0.1239 (13) fm and 280 configurations [24,25]. We choose an ensemble with small N L = 16 and L 2 fm so that Υπ(p z ) with p z > 2 2π L appear at E > m B + m B * + 0.2 GeV above our interest; larger L would require further operators like O 4,5 with higher p.

Calculation of eigen-energies and overlaps:
are evaluated using the full distillation method [23] andbb annihilation Wick contraction is omitted as in practically all quarkonium lattice studies. C ij are averaged over 8 3 or 16 3 space positions ofb, while sub-matrix for O 3−6 is averaged over all source time slices to increase accuracy. Eigen-energies E n and overlaps O i |n are extracted from 6 × 6 matrices C ij (t) = n O i |n e −Ent n|O † j using the widely used GEVP variational approach [26][27][28].
Eigen-energies ofbbdu system as a function of r: The main result of our study are the eigen-energies of thē bbdu system ( Fig. 1a) with static b andb separated by r. They are shown by points in Figure 2. The colors of points indicate which Fock-component (1) dominates an eigenstate, as determined from normalized overlaps of an eigenstate |n to operators O i . Normalized overlapZ n i ≡ O i |n / max m O i |m is normalized so that its maximal value for given O i across all eigenstates is equal to one.
The dashed lines in Fig. 2 provide related noninteracting (n.i.) energies E n of two-hadron states (1)  to m B + m B * for r > 0.5 fm, but it has significantly lower energy for r [0.1, 0.4] fm (red circles in Fig. 2). This indicates sizable strong attraction between B and B * in this system -something that might be related to the existence of Z b tetraquarks. This is the most important and robust result of this lattice study. Other eigenstates are dominated by Υπ( p) and Υb 1 . Their energies E lie close to non-interacting energies E n.i. (4) given by dot-dashed lines, so E E n.i. . We point out that we can not claim nonzero energy shifts E − E n.i. for Υπ and Υb 1 states (although Fig. 2 shows small deviations from zero in some cases) since the statistical and systematic errors are not small enough.
Towards masses of Z b states within certain approximations: Eigen-energies ofbbdu system in Fig. 1a indicate that eigenstate dominated by the BB * has significantly lower energy than m B + m B * at small separation r between static b andb. This suggests a possible existence of exotic hadrons (related to Z b ) and peaks in cross-section near BB * threshold. Such physical observables require the study of motion of heavy degrees of freedom based on energies E n (r) according to Born-Oppenheimer approach. The precise prediction of such observables is not possible at present since lattice eigen-energies are not known for r < a. In addition, the accurate study would require coupled-channel treatment of all Fock components (1) through coupled-channel Schrödinger equation, which  is a challenging task left for the future (this was recently elaborated in [29] for conventionalbb with I = 0). We apply two simplifying approximations in order to shed light on the possible existence of Z b based on energies in Figure 2. The first assumption is that the eigenstate indicated by red circles in Fig. 2 is related exclusively to BB * Fock component and does not contain other Fock components in (1). This is supported by our lattice results to a very good approximation, since this eigenstate couples almost exclusively to O BB * and has much smaller coupling to O Υπ and O Υb1 : the normalized overlap of this state to O Υπ,Υb1 isZ 3−6 ≤ 0.02 for r ≤ 3, while overlap to O BB * isZ 1,2 O(1). In the reminder we explore physics implications of this eigenenergy E BB * (r).
The energy E BB * (r) represents the total energy without the kinetic energy of heavy degrees of freedom. The difference V (r) = E BB * (r) − m B − m B * therefore represents the potential felt by the heavy degrees of freedom, in this case between B andB * mesons. The extracted potential is plotted in Fig. 3. The potential shows sizable attraction for r = [0.1, 0.4] fm and is compatible with zero for r ≥ 0.6 fm within sizable errors. Lattice study that would probe whether one-pion exchange dominates at large r would need higher accuracy.
The problem is that potential V (r) is not determined from lattice for r < a, it might be affected by discretization effects at r a and the analytic form of r-dependence is not known apriori. This brings us to the second simplifying approximation where we assume a certain form of the regular potential V reg. (r) that has no singularity at r → 0. The fits of lattice potential for various choices of power F and Here W = E tot −m B −m B * is the energy with respect to BB * threshold. The B and B * can couple to Z b channel with J P = 1 + in partial waves l = 0, 2. Below we extract (virtual) bound states and scattering rates for l = 0, while l = 2 is not discussed since V (r) + l(l+1) 2µr 2 > 0 is repulsive for all r. Wave functions of the Schrödnger equation render the phase shift δ l=0 (W ) and BB * scattering matrix S(W ) = e 2iδ0(W ) . Resonances above threshold do not occur for purely attractive s-wave potentials since there is no barrier to keep the state metastable, while (virtual) bound states below threshold may be present. Bound state (virtual bound state) corresponds to the pole of S(W ) for real W < 0 and imaginary momenta k = i|k| (k = −i|k|) of B andB * in the center of momentum frame.
We find a virtual bound state below threshold and its location is shown by diamonds in Fig. 4; this pole is present when the parameter F in V (5) is F<1.9. If this pole is close below threshold, it enhances the BB * crosssection above threshold. For example, we find a virtual bound state with mass M slightly below threshold for the values of parameters F = 1.3, A = 1.139(50), d = 1.615(71) in V (5). This state is responsible for a peak in the BB * rate N BB * ∝ kσ ∝ sin 2 δ 0 (W )/k above threshold, shown in Fig. 5 for the central value of parameters. The shape of the peak resembles the Z b (10610) peak in the BB * rate observed by Belle (Fig. 2 of [3]).
The significantly attractive BB * potential (Fig. 3) and the resulting virtual bound state (diamonds in Fig. 4) could be related to the existence of Z b in experiment. The reliable relation between both will be possible only when simplifications employed here will be overcome in the future simulations. We note that Z b (10610) was found as a virtual bound state slightly below threshold also by the re-analysis of the experimental data [4] when the coupling to bottomonium light-meson channels was turned off [4] (the position of the pole is only slightly shifted when this small coupling is taken into account).
Surprisingly, the strongly attractive potential V (r) (5) leads also to a deep bound state at M − m B − m B * = −411 ± 20 MeV. Such a state was never reported by experiments. If our approach can be trusted for such deep bound states and if such a bound state exists, it could be searched for in Z b → Υ(1S)π + decays. The invariant mass distribution observed by Belle is indeed not flat (Fig. 4a of [2]) and it would be valuable to explore if some structure becomes significant at better statistics.
The exotic Z b resonances were observed only by Belle, so their confirmation by another experiment would be highly welcome. LHCb could try to search for it in inclusive final state BB * .

Conclusions:
We presented a lattice QCD study of a channel with quark structurebbdu, where Belle observed two exotic Z b hadrons. We find significantly attractive potential V (r) between B andB * at small r when the total spin of the heavy quarks is equal to one. Dynamics of BB * system within the extracted V (r) leads to a virtual bound state, whose mass depends on the parametrization of V . Certain parametrizations give a virtual bound state slightly below BB * threshold and a narrow peak in BB * rate just above threshold, resembling Z b in experiment.
For quantitative comparison to experiment, future lattice studies need to explore how the dynamics of BB * is influenced by the coupling to Υπ channels, and by the component where the total spin of the heavy quarks is equal to zero. Derivation of the appropriate analytic form for V (r) would be very valuable.

S1: SYMMETRIES AND OPERATORS
Here we provide more details on the transformation properties of the investigated system in Fig. 1(a)  The quantum number = −1 is related to the reflection of light-degrees of freedom over yz plane, which is a product of rotation R x,π by π around x and inversion I with respect to the midpoint between 0 and r. The light-quark partq b (1 − γ t )γ 5 q a of O 1 is invariant under rotations and has P = −1, therefore = −1. The pion with momenta in z-direction within O 4 transforms as π p= ez Rx,π −→ π − ez I −→ −π ez , while the straight gauge link is invariant under this reflection, so = −1.
The Dirac structure for the heavy quark part in all operators isbγ z P + b, which ensures S h = 1 and S h z = 0. The C · P = −1 is related to the product of chargeconjugation and inversion with respect to the mid-point between 0 and r. Both refer to the transformation of the light-degrees of freedom as well as the transformation of the static colour sources 1 . This is most conveniently accomplished by the usual transformation rules ψ C → Cψ T and ψ P → γ t ψ for both ψ = q and b, where this operation does not affect the heavy quark spin, while C = iγ 2 γ t . The operator O 1 has C ·P = −1 since 1 If the color of the static source was not transformed under the charge-conjugation, the color-singletb a q aqb b b would transform under C-conjugation tob a Cq aT q bT Cb b , which is not gauge invariant.
where P exchanges positions 0 and r, γ t CΓ T Cγ t =Γ forΓ = γ z P + , γ t CΓ T Cγ t = −Γ for Γ = P − γ 5 , and dummy indices a ↔ b can be exchanged in the last expression. The linear combination π p= ez + π p=− ez The second line in operator O 1 (3) is obtained via the Fierz rearrangement and further simplifies since static heavy quarks appear in the combination P + b andbP − .

S2: EFFECTIVE ENERGIES AND OVERLAPS
Effective energies E ef f n of the system in Fig. 1a are shown in Fig. S1(a) for separation r/a = 2 and all eigenstates n = 1, ..6. Effective energies are obtained from correlation matrices C ij (t) via variational approach C(t)u n (t) = λ n (t)C(t 0 )u n (t), where effective energies are given by the eigenvalues E ef f n (t) ≡ ln[λ n (t)/λ n (t + 1)]. Reference time t 0 = 2 is used and agreement for t 0 = 3, 4 is verified. Effective energies render eigen-energies E n in the plateau region, indicated in the plots.
The overlaps O i |n of each eigenstate n to employed operators O i (3) are shown in terms of the normalized overlapsZ n i in Fig. S1(b). HereZ n i ≡ O i |n / max m O i |m is normalized so that its maximal value for given O i across all eigenstates is equal to one. Overlaps of the eigenstate dominated by BB * (red circles in Fig. 2) are in Fig. S2.
The lattice potential V (r) between B andB * from Fig.  3a is tabulated in Table S1.  S1: Potential V (r) between B andB * extracted from our simulation and plotted in Fig. 3a. Potential V (r) for separations r/a > 4 is equal to zero within statistical and systematic errors.  Here we consider the potential between B andB * analytically, where b andb are separated by a very small distance r r B , such that r is much smaller than average distance r B between b andq in B ( * ) meson (i.e. average radius r B of a static B ( * ) meson). We address the question whether this potential has a singular form V 1/r (r) = K r for r → 0 and determine prefactor K, while we omit all sub-leading contributions that are finite at r → 0. Among all pairs of the four quarksbbqq, only interaction between b andb at very small r could give potential proportional to 1/r. All other pairs are at average distance of the order of O(r B ), which is finite for r → 0; these pairs do not lead to infinite potential for r → 0 and we therefore omit their contribution to V 1/r . The task is therefore to determine potential between b andb within a pair of color-singlet B ( * ) mesons The color structure with indices a and b matches with employed operators O BB * , while other indices will not be relevant below. In order to determine the potential, |BB * is expressed in terms of color singlets and octets Thebb singlet within BB * |BB * renders singlet potential V 0 (r), all eight octets render octet potential V 8 (r) 2q λ A q → 1 are properly normalized to one. The resulting BB * potential at very small r is therefore The singlet and octet contributions cancel in the case of one-gluon exchange, i.e at the order O(α s ). The lowest non-zero contribution can be obtained from the perturbative calculation of both potentials in [30] and comes at O(α 3 s ) with δa 2 = −189.2. Employing the value of α s 0.31 obtained from the fit of the singletbb static potential in our simulation, we arrive at V 1/r (r) −0.0051/r. The masses of (virtual) bound states in Fig. 4 of the main article were based on the fits of the lattice potential in the range r/a = [1,4] and form of potential V (r) = V reg. (r) + V 1/r (r) in Eq. (5): these masses are shown again in Fig. S3(a) for completeness. The lattice potential at r/a = 1 can be prone to lattice discretization errors, therefore we investigate in sensitivity of the results on excluding this point from the fit. The masses based on the fits in the range r/a = [2,4] are shown in Fig. S3(b). They also support the presence of a virtual bound state for values of parameter F < 2 and a deep bound state 350 − 400 below threshold, so the results qualitatively agree in both cases.
The part of the potential V 1/r (r) that is singular as r → 0 was determined perturbatively (S5) for very small separations r. It is equal to zero at the one-gluon exchange level and the lowest non-zero contribution comes at the order O(α 3 s ). The sensitivity of the results on including or excluding this part of the potential is explored in Fig. S3. The masses based solely on the regular potential V reg. in Figs. S3(c,d) agree within errors with masses based on V reg. + V 1/r in Figs. S3(a,b). This agreement is a consequence of the suppression in V 1/r . are also compared.