Construction of bb u ¯ d ¯ tetraquark states on lattice with NRQCD bottom and HISQ up/down quarks
aa r X i v : . [ h e p - l a t ] A ug Construction of bb ¯ u ¯ d tetraquark states on latticewith NRQCD bottom and HISQ up/down quarks Protick Mohanta and Subhasish Basak
School of Physical Sciences, National Institute of ScienceEducation and Research, HBNI, Odisha 752050, India (Dated: August 26, 2020)
Abstract
We construct bb ¯ u ¯ d states on lattice using NRQCD action for bottom and HISQ action for the lightup/down quarks. The NRQCD-HISQ tetraquark operators are constructed for “bound” [ bb ][¯ u ¯ d ] and“molecular” [ b ¯ u ][ b ¯ d ] states. Corresponding to these different operators, two different appropriatelytuned light quark masses are needed to obtain the desired spectra. We explain this requirementof different m u/d in the light of relativised quark model involving Hartree-Fock calculation. Themass spectra of double bottom tetraquark states are obtained on MILC N f = 2 + 1 Asqtad latticesat three different lattice spacings. Variational analysis has been carried out to obtain the relativecontribution of “bound” and “molecular” states to the energy eigenstates. . INTRODUCTION The multiquark hadronic states other than the mesons and baryons are relatively newentrants particularly in the heavy quark sector. The signature of some of such states con-taining four or more quarks and/or antiquarks have been observed in experiments [1–6].Such states are characterized by J P C quantum numbers that cannot be arrived at fromthe quark model. (Such four quark states are popularly referred to as tetraquark, whichis used to denote either bound or often both bound and mesonic molecular states. In thispaper we use the term tetraquark in the latter sense.) However, heavy hadronic tetraquarkstates QQ ¯ l ¯ l and their stability in the infinite quark mass limit had been studied in [7, 8]which raised the possibility of existence of heavy four quark bound states below the Q ¯ l − Q ¯ l threshold. Of late, the observations of Z − (4430 , + ) of minimal quark content being c ¯ cd ¯ u [3]have been reported along with the 1 + states like Z b (10610) and Z ′ b (10650), having minimalquark content of four quarks (containing a b ¯ b pair) that are a few MeV above the thresholdsof B ⋆ ¯ B (10604 .
6) and B ⋆ ¯ B ⋆ (10650 .
2) [1, 2]. The proximity of Z b , Z ′ b to the B ⋆ ¯ B ⋆ thresholdvalues perhaps suggest molecular, instead of bound, nature of the states.Around the same time, lattice QCD has been employed to investigate the bound and/ormolecular nature of the heavy tetraquark states, not only to understand the above experi-mentally observed states but also to identify other possible bound tetraquark states in both0 + and 1 + channels. In the charm sector, some early lattice studies involve T cc and T cs tetraquark states [9], cc ¯ c ¯ c [10], X (3872) and Y (4140) [11] and more recently D ⋆s (2317) [12].The bottom sector received intense attention where, instead of B ∗ ¯ B or B ∗ ¯ B ∗ , relativelysimpler BB, BB ∗ systems are studied. The lattice investigations this far involve four bot-tom bb ¯ b ¯ b [13] and two bottom tetraquark states ¯ b ¯ bl l , where l , l ∈ c, s, u, d , [14–17]. Animportant observation of these lattice studies is that the possibility of the existence of bb ¯ l ¯ l tetraquark bound states increases with decreasing light quark masses, while they becomeless bound with decreasing heavy (anti)quark mass.Besides the usual lattice simulations, the heavy tetraquark systems have also been studiedusing QCD potential and Born-Oppenheimer approximation [18–21]. The main idea in thismethod is to investigate tetraquark states with two heavy (anti)quarks, which was ¯ b ¯ b in thestudy, and two lighter quarks using quantum mechanical Hamiltonian containing screenedCoulomb potential. This approach has been used to explain our two different choices of light2 /d quark masses for different classes of tetraquark states.In this work our goal is to construct tetraquark states, having quark content bb ¯ l ¯ l in1 + both below and above B − B ∗ threshold, by a combination of lattice operators andtuning quark masses based on quantum mechanical potential calculation. For the b quark,we employed nonrelativistic QCD formulation [22, 23], as is the usual practice, and HISQaction [24] for l , l = u/d . Here we also explore through variational/GEVP analysis howthe trial states created by our operators contribute to the energy eigenstates.First, we briefly review the salient features and parameters of both NRQCD and HISQactions along with the steps involve in combining the relativistic u/d HISQ propagators withthe NRQCD b quark propagators in the section II. We have considered two different kind ofoperators – the local heavy diquark and light antidiquark (often referred as “good diquark”configuration) and molecular meson-meson, we described these constructions in the sectionIII. We collect our spectrum results in section IV that contains subsections on quark masstuning (IV A), Hartree-Fock calculation of two light quark in the presence of a heavy quark(IV B), tetraquark spectra (IV C) and GEVP analysis (IV D). Finally we summarized ourresults in section V. II. QUARK ACTIONS
Lattice QCD simulations with quarks require quark mass to be am l ≪
1, where a is thelattice spacing. In the units of the lattice spacings presently available, the b quark mass isnot small i.e. am b ≮
1. As is generally believed, the typical velocity of a b quark inside ahadron is nonrelativistic v ∼ . b quarks on lattice. We have used O ( v ) NRQCD action[23], where the Hamiltonian is H = H + δH , where H is the leading O ( v ) term, the O ( v )and O ( v ) terms are in δH with coefficients c through c , H = − ˜∆ m b − a n (∆ ) m b (1) δH = − c (∆ ) m b + c ig m b (cid:16) ~ ∆ ± · ~E − ~E · ~ ∆ ± (cid:17) − c g m b ~σ · (cid:16) ˜ ~ ∆ ± × ˜ ~E − ˜ ~E × ˜ ~ ∆ ± (cid:17) − c g m b ~σ · ˜ ~B − c g m b n ∆ , ~σ · ~B o − c g m b n ∆ , ~σ · (cid:16) ~ ∆ ± × ~E − ~E × ~ ∆ ± (cid:17)o − c ig m b ~σ · ~E × ~E (2)3here ∆ ± and ∆ are discretized symmetric covariant derivative and lattice Laplacian re-spectively. Both the derivatives are O ( a ) improved as are the chromoelectric ~E and chromo-magnetic ~B fields. The b quark propagator is generated by time evolution of the Hamiltonian H , G ( ~x, t + 1; 0 ,
0) = (cid:18) − aH n (cid:19) n (cid:18) − aδH (cid:19) U ( ~x, t ) † (cid:18) − aδH (cid:19) (cid:18) − aH n (cid:19) n G ( ~x, t ; 0 , G ( ~x, t ; 0 ,
0) = t < δ ~x, for t = 0The tree level value of all the coefficients c , c , c , c , c , c and c is 1. Here n is the factorintroduced to ensure numerical stability at small am b , where n > / m b [22].In NRQCD, the rest mass term does not appear either in the equation (1) or in (2),and therefore, hadron masses cannot be determined from their energies at zero momentumdirectly from the exponential fall-off of their correlation functions. Instead, we calculate thekinetic mass M k of heavy-heavy mesons from its energy-momentum relation, which to O ( p )is [25], E ( p ) = E (0) + q p + M k − M k ⇒ E ( p ) = E (0) + E (0) M k p . (4)where M k is the kinetic mass of the meson which is calculated considering E ( p ) at differentvalues of lattice momenta ~p = 2 ~nπ/L . The b quark mass is tuned from the spin average ofkinetic masses of Υ and η b , and matching them with the experimental spin average value, M b ¯ b = 3 M Υ + M η b M b ¯ b is tuned to, however, is not 9443 MeV that is obtainedfrom spin averaging Υ (9460 MeV) and η b (9391 MeV) experimental masses, but to anappropriately adjusted value of 9450 MeV [26], which we denote as M modphys in the equation(6) below. The hadron mass is then obtained from M latt = E latt + n b (cid:0) M modphys − E η b latt (cid:1) (6)where E latt is the lattice zero momentum energy in MeV, n b is the number of b -quarks inthe bottom hadron. 4he u/d light quarks comfortably satisfy the criteria am l ≪ u/d quarks, which is given in [24], S = X x ¯ q ( x ) (cid:0) γ µ D HISQ µ + m (cid:1) q ( x ) where, D HISQ µ = ∆ µ ( W ) − a ǫ ) ∆ µ ( x ) . (7)Because HISQ action reduces O ( α s a ) discretization error found in Asqtad action, it is wellsuited for u/d (and s ) quarks. The parameter ǫ in the coefficient of Naik term can beappropriately tuned to use the action for c quarks, which we do not have here. For u/d (and s ) quarks, the ǫ = 0.HISQ action is diagonal in spin space, and therefore, the corresponding quark propagatorsdo not have any spin structure. The full 4 × x ) = Y µ =1 ( γ µ ) x µ = γ x γ x γ x γ x . (8) III. TETRAQUARK OPERATORS
In the present paper, we have considered two kinds of tetraquark operators – the lo-cal heavy diquark and light antidiquark and molecular meson-meson. The b quark, beingnonrelativistic, is expressed in terms of two component field ψ h . We convert it to a fourcomponent spinor Q having vanishing two lower components, Q ≡ ψ h (9)which help us to combine b field and relativistic four component light quark fields in the usualway. The heavy-light meson operator, that we will make use of in the operator construction,is written as O hl ( x ) = ¯ Q ( x ) Γ l ( x ) (10)where l ( x ) stands for the light quark fields, ¯ Q = Q † γ and depending on pseudoscalar andvector mesons Γ = γ and γ i respectively.Because of the vanishing lower components, the states with Q can only be projectedto the positive parity states. The local double bottom tetraquark operators that we can5onstruct for bb ¯ l ¯ l system are, O M ≡ O B ∗ B = (cid:2) ¯ l ( x ) γ i Q ( x ) (cid:3) (cid:2) ¯ l ( x ) γ Q ( x ) (cid:3) , (11) O M ≡ O B ∗ B ∗ = ǫ ijk (cid:2) ¯ l ( x ) γ j Q ( x ) (cid:3) (cid:2) ¯ l ( x ) γ k Q ( x ) (cid:3) , (12) O D ≡ O Q ∗ ˜ π = (cid:2) Q T ( x ) Cγ i Q ( x ) (cid:3) (cid:2) ¯ l ( x ) Cγ ¯ l T ( x ) (cid:3) (13)where l = l and l , l ∈ u, d . The naming convention above is borrowed from reference [15]but the exact construction of the operators is different. In literature the operators in (11)and (12) are often referred to as “molecular”. The quark fields within the square bracketsare color contracted. The diquark-antidiquark 1 + four quark state bb ¯ l ¯ l with l = l in (13)can actually be defined in two ways [28], O Q ∗ ˜ π = (cid:2) Q a T Cγ i Q b (cid:3) (cid:2) ¯ l a Cγ ¯ l b T − ¯ l b Cγ ¯ l a T (cid:3) O Q ˜ π ∗ = (cid:2) Q a T Cγ Q b (cid:3) (cid:2) ¯ l a Cγ i ¯ l b T + ¯ l b Cγ i ¯ l a T (cid:3) (14)where a, b are color indices. The subscripts Q ∗ and ˜ π in the operator O Q ∗ ˜ π are in ¯3 c and3 c respectively, while Q and ˜ π ∗ in the operator O Q ˜ π ∗ are in 6 c and ¯6 c . But both O Q ∗ ˜ π and O Q ˜ π ∗ correspond to the 1 + state. Of these the O Q ∗ ˜ π is our desired “bound” tetraquarkoperator because one-gluon-exchange interaction is attractive for a heavy quark pair in ¯3 c diquark configuration [8] and spin dependent attraction exists for light quark pairs in “gooddiquark” configuration characterized by color ¯3 c , spin J = 0 and isospin I = 0 or 1 / O Q ∗ ˜ π contribute identically in the final correlator, hence we consider onlythe first term in the calculation. The generic form of the temporal correlation among theoperators at zero momentum is, C XY ( t ) = X x (cid:10) [ O X ( x , t )] [ O Y ( , † (cid:11) , (15)where X, Y can be any of
D, M , M in equations (11), (12) and (13). For example, theexplicit forms of the zero momentum correlators, including cross-correlator, when X and Y M = B ∗ B and D = Q ∗ ˜ π , are C M M ( t ) = X ~x Tr h γ M † ( x, γ γ i G ( x, γ i i × Tr h M † ( x, G ( x, i (16) C D D ( t ) = X ~x Tr h(cid:0) G ad ( x, (cid:1) T γ i γ γ G bc ( x, γ γ γ i i × Tr (cid:20) γ γ M † da ( x, γ γ (cid:16) γ M † cb ( x, γ (cid:17) T (cid:21) (17) C D M ( t ) = X ~x Tr h G ad ( x, γ i γ M † da ( x, γ γ γ γ i × Tr h γ γ i γ G bc ( x, γ γ M † cb ( x, γ i (18)Above in the equations (17) and (18), traces and transposes are taken over the spin indices,while in the equation (16) the traces are taken over both the spin and color indices. Here G ( x,
0) denotes the heavy quark propagators while the M ( x,
0) are the light quark propaga-tors. In terms of coding, we have to keep in mind that NRQCD and MILC library suit usesdifferent representation of gamma matrices. Therefore the heavy quark propagator G ( x, S γ
MILC µ S † = γ NR µ where, S = 1 √ σ y σ y − σ y σ y . (19) IV. NUMERICAL STUDIES
We calculated the double bottom tetraquark spectra using the publicly available N f =2 + 1 Asqtad gauge configurations generated by MILC collaboration. Details about theselattices can be found in [31]. It uses Symanzik-improved L¨uscher-Weisz action for the gluonsand Asqtad action [32, 33] for the sea quarks. The lattices we choose have a fixed ratio of am l /am s = 1 / A. Quark mass tuning
For bb ¯ u ¯ d mass calculation, we need nonperturbative tuning of both m b and m u/d . Withthe help of equation (5), the tuning of m b has been carried out by calculating the spin7ABLE I: MILC N f = 2 + 1 Asqtad configurations used in this work. The gauge couplingis β and the lattice spacing is a . The u/d and s sea quark masses are m l and m s respectively and the lattice size is L × T . The N cfg is number of configurations used inthis work. β = 10 /g a (fm) am l am s L × T N cfg ×
48 6006.76 0.12 0.0100 0.0500 20 ×
64 6007.09 0.09 0.0062 0.0310 28 ×
96 300 average Υ and η b kinetic masses and comparing the same with the spin average and suitablyadjusted experimental Υ and η b masses as discussed in section II. The tuned bare am b quarkmasses for lattices used in this work are given in Table II.TABLE II: Tuned b and u/d quark bare masses for lattices used in this work. For u/d -quark mass, we mention the particle states used to tune. Quark Tuning 16 ×
48 20 ×
64 28 × am b Υ − η b am u/d Λ b (5620) 0.105 0.083 0.064 am u/d B (5280) 0.155 0.118 0.087 But the tuning of am u/d is rather tricky and the way to achieve this numerically has beendiscussed in the study of bottom baryon spectra in [30]. In the present paper, we try tounderstand the light quark tuning in more details with the help of relativised quark model[34, 35] and Hartree-Fock calculation. The basic idea is that m u/d has to be tuned to twodifferent values corresponding to two different construction of the pairs [¯ u ¯ d ] and [ b ¯ u ]. In theoperator O D ≡ [ bb ][¯ u ¯ d ], the antidiquark part formed with two light u/d quarks is same asin Λ b ≡ ( u T Cγ d ) b , and therefore, we use experimental Λ b mass 5620 MeV to tune the bare am u/d . For the operators O M /M , the diquark part is formed between heavy quark and lightantiquark [ b ¯ u ] which is same as in the B -meson (¯ bγ (5 ,k ) u ) or Σ b ≡ ( Q T Cγ u ) u . In such case8e tend to use B -meson mass 5279 MeV to tune the am u/d . The result of u/d quark masstuning that are made use of in this work is given in the Table II.Before we calculate the spectra of the D, M , M tetraquark states, in the followingsubsection IV B we try to understand the diquark dependent different tuning of the light u/d quark masses. For this we consider Schr¨odinger Hamiltonian for ( a ) hydrogen-likesystem, namely B meson with an ¯ u antiquark in the potential of a static b quark and ( b )helium-like system, which is Λ b baryon with ¯ u, ¯ d quarks in the same b quark field. B. Hartree-Fock calculation of tetraquark states
In order to gain a qualitative understanding of two different tunings of m u/d , we considerthe light antiquark and light-light diquark in the potential of heavy, nearly static colorsource, the b quark(s). This picture is akin to hydrogen and helium-like quantum mechanicalsystems. For the molecular tetraquark states, the basic assumption is that the light antiquarkwave functions do not have significant overlap with each other and they are effectively inthe potential of their respective heavy b quarks [19] i.e. a two B -meson like system. But forthe diquark-antidiquark tetraquark state [ bb ][¯ u ¯ d ] where antidiquark component is similar tothe Λ b light-light diquark, we tune the u/d quark mass using the Λ b baryon. The situationis depicted schematically in Figs. 1 and 2. The relevant interpolating operators for Λ b and B meson are fairly standard but for HISQ light quarks a whole array of bottom baryonoperators, including Λ b can be found in [30]. b b ¯ u ¯ d B − B FIG. 1: Molecular tetraquark state viewed as bound state of two B mesons, which issimilar to two hydrogen atoms forming a hydrogen molecule.9he relativised quark model [34, 35] helps us to numerically calculate the masses of B mesonand Λ b baryon using the light (anti)quark mass as parameter. The molecular tetraquarkstate can be visualized as two B meson molecule as shown in the Fig. 1. Then for each B meson, the light u/d antiquark is taken to be in the field of “static” b quark and we solvethe problem by considering the radial part of the Schr¨odinger equation numerically usingsuitably modified Herman-Skillman code [36]. − m u/d d U ( r ) dr + V ( r ) U ( r ) = EU ( r ) (20)Here U ( r ) = rψ ( r ) and the potential V ( r ) is given by V ( r ) = − α r + βr (21)The B meson mass M B is, therefore, determined from the energy eigenvalue E , M B = m b + m u/d + E (22)where m b = 4 .
18 GeV (MS) is the mass of the bottom quark, the α = π/
16 [20] and β = 0 . [34]. For M B = 5 .
279 GeV, the light quark mass obtained is m u/d ≈ .
227 GeV. bu d Λ b (a) u/d quarks in Λ b baryons form a 3 c diquark in presence of a b quark. b b ¯ u ¯ d [ bb ][¯ u ¯ d ] (b) Like Λ b , two b quarks form a (nearly)static nucleus surrounded by ¯ u, ¯ d cloud. FIG. 2: Schematic diagram of helium-like Λ b and [ bb ][¯ u ¯ d ] tetraquark state used forHartree-Fock treatment.For Λ b baryon, we used Hartree-Fock method [37, 38] to solve the helium-like Hamiltonian, H = − m u/d ∇ − α r + βr − m u/d ∇ − α r + βr − α ′ r + β ′ r r is the relative distance between two light quarks “orbiting” the two heavy quarksand their interaction potential is the last two terms in the equation (23) with coefficient α ′ and β ′ . For the Hartree-Fock calculation of the energy E , we take β ′ = β and α ′ = 0 . HF = 1 √ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) χ ( x ) χ ( x ) χ ( x ) χ ( x ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) , (24)where x i ≡ ( ~r, s ) collectively denotes the space and spin indices, χ i ( ~r, s ) = φ i s ( ~r ) S ( s ) with φ ( ~r ) being the 1S state. Therefore, the expectation value of the the Hamiltonian can bewritten as h Ψ HF | H | Ψ HF i = h T i + Z ρ ( ~r ) V ext ( ~r ) d~r − Z ′ Z Z ρ ( ~r ) ρ ( ~r ) | ~r − ~r | d~r d~r + B ′ Z Z ρ ( ~r ) ρ ( ~r ) | ~r − ~r | d~r d~r + Z ′ X i,j,s Z Z φ ⋆i s ( ~r ) φ ⋆j s ( ~r ) φ i s ( ~r ) φ j s ( ~r ) | ~r − ~r | d~r d~r − B ′ X i,j,s Z Z φ ⋆i s ( ~r ) φ ⋆j s ( ~r ) φ i s ( ~r ) φ j s ( ~r ) | ~r − ~r | d~r d~r (25)where, we have used h T i = X i,s (cid:28) φ i s ( ~r ) (cid:12)(cid:12)(cid:12)(cid:12) − m u/d ∇ (cid:12)(cid:12)(cid:12)(cid:12) φ i s ( ~r ) (cid:29) ρ ( ~r ) = X i,s | φ i s ( ~r ) | , V ext ( ~r ) = − α r + βr Z ′ = 2 α ′ B ′ = β ′ . In contrast to the helium atom, the presence of linear r -terms in the Hamiltonian leadsto additional exchange-energy terms in the calculation. With these linear r -terms in, theHartree-Fock equation becomes E φ i s ( ~r ) = (cid:20) − m u/d ∇ + V ext ( ~r ) − Z ′ Z ρ ( ~r ) | ~r − ~r | d~r + B ′ Z ρ ( ~r ) | ~r − ~r | d~r (cid:21) φ i s ( ~r ) − B ′ X j,s Z φ ⋆j s ( ~r ) φ i s ( ~r ) φ j s ( ~r ) | ~r − ~r | d~r + Z ′ X j,s Z φ ⋆j s ( ~r ) φ i s ( ~r ) φ j s ( ~r ) | ~r − ~r | d~r (26)11e solve for E in equation (26) iteratively and, eventually, the Λ b mass is calculated from M Λ b = m b + 2 m u/d + E (27)The PDG value of Λ b (5620) is obtained by setting the m u/d to 0.157 GeV.TABLE III: Comparison of m u/d obtained from various lattices with quark massparameters in the equations (20) and (26). Lattice B meson: m u/d = 227 MeV Λ b baryon: m u/d = 157 MeV am u/d m u/d (MeV) am u/d m u/d (MeV)16 ×
48 0.155 204 0.105 13820 ×
64 0.118 194 0.083 13728 ×
96 0.087 191 0.064 143
In Table III, we compare the nonperturbatively tuned m u/d on our lattices with those ob-tained by solving the equations (20) and (26). The bare lattice light quark masses cannotbe directly compared to the parameter m u/d in these equations mainly because of the useof renormalized b quark mass (in MS scheme) in the Hartree-Fock calculation. Therefore,the m u/d ’s in the above calculation return a sort of “renormalized constituent” quark mass.Nonetheless it is obvious that we need two different m u/d for two different systems, namely B and Λ b . So by comparing the two sets, we simply wish to point out that the lattice tuned m u/d ’s are in same order of magnitude as Schr¨odinger equation based quark model but havea difference of 10 – 15%. This helps us to understand the possible physics behind two dif-ferent tuning of light quark mass in determining the masses of single bottom hadron(s) anddouble bottom tetraquark states. C. bb ¯ u ¯ d spectrum A plot of variation of bb ¯ u ¯ d mass with various am u/d , including the Λ b and B tuned valuesis shown in Fig. 3. Here we make a naive comparison of our data with the earlier quarkmodel, lattice calculations and the PDG values, and it shows an interesting trend.12 M a ss i n M e V a m q D η s Λ b B B s [ bb ][¯ u ¯ d ] this work [ b ¯ u ][ b ¯ d ] this workPDG Z ′ b PDG Z b Ei hten et al.Leskove et al.Junnarkar et al.Fran is et al.Karliner et al.
FIG. 3: Variation of bb ¯ u ¯ d mass at various am u/d in 16 ×
48 lattice. Λ b -tuned tetraquarkstates almost overlap with many of the quark model and lattice calculations, namelyEichten et al. [8], Leskovec at al. [17], Junnarkar et al. [16], Francis et al. [14, 15],Karliner at al. [39]. The B -tuned states instead coincide with Z b , Z ′ b PDG results [40].Firstly, PDG Z b , Z ′ b and the lattice results are clustered around two different masses. Ourdata at B meson tuning point coincides with the PDG Z b (10610) and Z ′ b (10650) statesaligning with the idea that they decay mostly into ¯ BB ∗ and ¯ B ∗ B ∗ respectively, possiblyindicating molecular nature of the state. However, our tetraquark state with Λ b tuningoverlaps mostly with other lattice results indicating the possibility of capturing a boundtetraquark state [ bb ][¯ l ¯ l ] much like the b [ ll ] state of Λ b . The effective masses of the states,obtained from O D and O M , when compared with the B − B ∗ threshold, we find D to exhibita shallow bound state while M is just marginally above. The majority of the lattice results[14–17] are found to be below this threshold as is obvious from the Fig. 3.To this end, in Fig. 4 we plot the effective masses of these two states obtained at differentlattice spacings together with the lattice thresholds for easy comparison. The colored bandsrepresent fitted am eff values. The superscripts Λ b and B denote the light quark tuning.13 .01.21.41.61.82.02.22.42.62.8 0 10 20 30 40 50 E (cid:27) e t i v e m a ss a m e ff t −→ × O Λ b D × O BM × O Λ b D × O BM × O Λ b D × O BM FIG. 4: Effective mass plot of the states of the operators O D and O M calculated on 16 × ×
64 and 28 ×
96 lattices. Dashed lines are B − B ∗ thresholds for different lattices. Foreasy viewing, the effective masses and thresholds on 20 ×
64 (purple colored) are multipliedby a common factor of 0.85, while that of 28 ×
96 (green colored) by 0.70.In the Table IV, we present our results of the tetraquark states corresponding to theoperators given in the expressions (11, 12, 13). We use two-exponential uncorrelated fit tothe correlation functions, the fitting range being chosen by looking at the positions of whatwe consider plateau in the effective mass plots. In the columns showing various lattices, wepresent the masses both in lattice unit a E latt and physical unit M latt in MeV, the notationsbeing introduced in equation (6). The errors quoted are statistical, calculated assuming thelattice configurations of different lattice spacings are statistically uncorrelated. The secondcolumn shows the tuning used for the corresponding states. In the last column we providethe masses averaged over all the lattice ensembles.From the Fig. 4 and Table IV it is clear that the trial state generated by our O D operatoris below B − B ∗ threshold which possibly indicates a bound state. On the other hand, thestates for O M and O M are just above it. We tabulate the difference of the masses from thetheir respective thresholds ∆ M D/M /M = M D/M /M − M B − M B ∗ in the Table V. In thistable, we calculated the following correlator ratio to determine the mass differences which14ABLE IV: Masses of tetraquark states for different am u/d tuning in lattice unit aE latt and M latt in MeV. We also include the B and B ∗ states that are used for threshold calculation. Operators Tuning 16 ×
48 20 ×
64 28 ×
96 Average aE latt M latt aE latt M latt aE latt M latt (MeV) O D = [ bb ][¯ u ¯ d ] Λ b O M = [ b ¯ u ][ b ¯ d ] B O M = ǫ ijk [ b ¯ u ] j [ b ¯ d ] k B B = b γ ¯ u B B ∗ = b γ k ¯ u B ∗ gives us an estimate of the binding energy, C X − B − B ∗ ( t ) = C X ( t ) C B ( t ) × C B ∗ ( t ) ∼ e − ( M X − M B − M B ∗ ) t (28)TABLE V: Mass differences of “bound” D and “molecular” M , M states from B − B ∗ threshold. The X subscript denotes any of the D, M , M . Operators Lattices a ∆ M X ∆ M X in MeV ∆ M X (MeV) O D × − . − − × − . − − × − . − − − − O M ×
48 0 . ×
64 0 . ×
96 0 . O M ×
48 0 . ×
64 0 . ×
96 0 .
15n the last column, we calculate our lattice average of ∆ M X in MeV and compare withsome of the previous lattice results. To our knowledge, the binding energies of the M , M states have been calculated in the framework of chiral quark model [41] for B − ¯ B ∗ and B ∗ − ¯ B ∗ states but there are no lattice results. But the binding energies for the first excitedstates, along with the ground states, obtained on different lattice ensembles are given in[17]. Though their tuning of light quark mass is very different compared to ours, still wecan use their result as a reference.Our binding energy for the bound tetraquark state D = [ bb ][¯ u ¯ d ] lies somewhere in themiddle of the previously quoted lattice results. Our statistical errors of the molecular states M and M are rather large but still they tentatively indicate non-bound molecular natureof the states. We will revisit the binding energy calculation for the molecular state(s) aftervariational analysis of the O M × O M correlation matrix.As we know, on lattice the operators for states having the same quantum numbers canmix and, therefore, a GEVP analysis can help resolve the issue of mutual overlap of variousstates on the energy eigenstates. In this work, rather than the energies of the eigenstates,we are more interested to learn the overlap of our trial states, namely D, M and M on thefirst few energy eigenstates, where | i is the ground state and | i , | i etc. are the excitedstates. D. Variational analysis
For the 2-bottom tetraquark system with quantum number 1 + , we consider the threelocal operators – “good” diquark O D , molecular O M and vector meson kind O M as definedabove in the expressions (11 – 13) – to capture the ground state ( | i , E ) and possibly thefirst excited state ( | i , E ). Apart from the construction itself, as we have discussed before,one basic difference in these operators is the use of different m u/d for simulation of thetetraquark states – Λ b tuning for O D whereas B tuning for O M /M operators.As is generally understood, these operators are expected to have overlap with the desiredground and excited states of the tetraquark system of our interest. The variational analysiscan be performed to determine the eigenvalues and the eigenvectors from the mixed statesformed by lattice operators. This is typically achieved by constructing a correlation matrix16nvolving the lattice operators O X , C XY ( t ) = D O X ( t ) O † Y (0) E = ∞ X n =0 h Ω |O X ( t ) | n ih n |O † Y (0) | Ω i e − E n t (29)where X, Y can either be all or any two combinations of
D, M , M in the expressions (11 –13). The terms h n |O † X | Ω i are the coefficients of expansion of the trial states O † X | Ω i , where | Ω i is the vacuum state, and written in terms of the energy eigenstates | n i as, O † X | Ω i = X n | n ih n |O † X | Ω i ≡ X n Z nX | n i (30)Presently, we are interested in expressing the energy eigenstates in terms of the trial statesto understand the contribution of each to the former. If we confine ourselves to the first fewenergy eigenstates, we can write | m i = X X v Xm O † X | Ω i ⇒ h l | m i = δ lm ≈ X X v Xm Z lX (31)The v Xm are equivalent to the eigenvector components obtained by solving a GEVP w.r.t asuitably chosen reference time t [12], C ( t ) v m ( t, t ) = λ m ( t ) C ( t ) v m ( t, t ) . (32)The eigenvalues λ m ( t ) are directly related to the energy of the m -th state, ground and thefirst few excited states, of our system through the relation λ m ( t ) = A m e −E m ( t − t ) (33)The component of eigenvectors v m ( t, t ) gives information about the relative overlap of thethree local operators to the m -th eigenstate. The eigenvectors v m ’s are normalized to 1.To determine the parameter t , we solve the GEVP and found that the ground and excitedstate energies are almost independent for t = 3 , , , B -tuned am u/d for all the operators O X but the resultsare similar with Λ b tuning and, hence, not shown. We chose t = 5 for our calculations.17 .002.20 3 5 7 9 E O M × O M E O M × O M E n e r g y e i g e n s t a t e s E E O D × O M E O D × O M t →E O D × O M t →E O D × O M × × × FIG. 5: Variation of ground and excited state energies E m of the equation (33) with t ,obtained by solving 2 × O M × O M , O D × O M and O D × O M . In this plot we used B -tuned am u/d for all the operators.TABLE VI: Comparing masses of Λ b -tuned D (from Table IV) with lowest two energyeigenstates in B -tuned O M × O M GEVP analysis.
Operators Energy 16 ×
48 20 ×
64 28 ×
96 ∆ M E-states a M a ∆ M a M a ∆ M a M a ∆ M (MeV) O D E − . − . − . O M × O M E E The GEVP analysis has been carried out in two steps because of differences in the tuningof am u/d for the “molecular” states M , M and “good” diquark state D . In the Table VI wehave shown our GEVP results of O M × O M correlation matrix and compare the differencesof energy level(s) from the threshold with that of D taken from Table IV. However, dueto large errors in calculating ∆ M in the first excitation E , we can only reliably quote thelowest energy E . Here the E and E are w.r.t. the O M × O M correlation matrix. The18round state energy E is significantly closer to the threshold than either of M or M . Nextwe look at the contribution of M and M in constructing the lowest energy eigenstate | i . | i , E = 2 .
063 (10) (16 × v M −→ v X v M −→ | i , E = 1 .
959 (12) (20 × v M −→ v X v M −→ | i , E = 1 .
888 (7) (28 × v M −→ v X v M −→ FIG. 6: The histogram plots of the normalized components ( v M , v M ) which define theenergy eigenstate | i = v M | M i + v M | M i .In the Fig. 6, we plot the histogram of the components of the normalized eigenvectors v =( v M , v M ) corresponding to the lowest energy E for all three lattices. Assuming that thecoefficients v M , M approximately remain the same on all time slices and for all the individualgauge configurations of an ensemble, the histogram figures are obtained by plotting the M , M components of normalized eigenvector v for all time points and individual gauge19onfigurations. As is expected, all three lattices return identical histogram of the coefficientsand hence, in the subsequent histogram plots we will show only the results from 28 × v M shows a peak around 0.9 indicating the lowest energy state | i receives dominant contribution from M trial state. We recall here that O M correspondsto the B − B ∗ molecular state as defined in the equation (11).However, the first excitation | i , for which our data is rather noisy to reliably estimate∆ M , the | M i and | M i states appear to have comparable contribution and are broadlydistributed over different time slices and vary significantly over configurations. This isevident from the histogram plot in the Fig. 7. This may have a bearing with the fact thatabove the threshold, the Z b tetraquark can couple to multiple decay channels resulting in abroad spectrum. | i , E = 1 .
906 (5) (28 × v M −→ v X v M −→ FIG. 7: The histogram plot of v M and v M that define the energy eigenstate | i = v M | M i + v M | M i .Including O D along with the O M and O M to form a 3 × b or B tuned am u/d in all three trial states. An important issue here is to interpretwhat Λ b tuning meant for B − B ∗ meson system or, conversely, what B tuning is meant forΛ b like system. Certainly, a B -tuned bound D state above threshold is not well-defined andwe find it has statistically small and varying overlap with the energy eigenstates much likein Fig. 7. On the other hand, Λ b tuned molecular states can possibly have finite overlap tothe eigenstates below threshold. However, we always expect dominance of D in | i becauseof the difference in construction of wave functions of the M and M .20 | i , E = 1 .
784 (12) (28 × v D −→ v M −→ v X v M −→ | i , λ = 1 .
816 (8) (28 × v D −→ v M −→ v X v M −→ | i , λ = 1 .
820 (22) (28 × v D −→ v M −→ v X v M −→ FIG. 8: Histogram plots of the normalized eigenvector components v X , v X and v X of 3 × X = D, M , M , on 28 ×
96 lattice.The histogram of the eigenvector components of 3 × ×
96 lattices. The lowest energy eigenstate | i is clearly dominated by D showing peak around 0.8, although it receives sizeable overlap from both M and M peaking around 0.45. But overlap of D on | i is rather small and it is mostly molecular M E is below the threshold. Our data for | i is too noisy toextract much information. Based on this Λ b tuned 3 × −
186 (15) MeV.
V. SUMMARY
In this work we have attempted to study two possible bb ¯ u ¯ d tetraquark states – one whichis bound and where most other lattice results are centered and, the other which is reported inPDG as Z b and Z ′ b . The experimentally observed states are believed to contain a b ¯ b pair butours is bb pair which is considered as theoretically simpler. However, the possible molecularnature of Z b and Z ′ b suggests that our molecular states should have similar masses. For thebottom quarks we have used NRQCD action while HISQ action for the u/d quarks. ThisNRQCD-HISQ combination has been employed earlier in [26] for bottom meson and recentlyin [30] for bottom baryons. We have constructed the three lattice 1 + trial states: a bound D containing “good diquark” 3 c configuration and two meson-meson molecular M , M withthe expectation that they will contribute to the states above B − B ∗ threshold. There arenot many lattice results on the states above threshold possibly because of the complicationthat they can couple to multiple decay channels besides B ∗ ¯ B and B ∗ ¯ B ∗ . Our motivationhere is to obtain a tentative estimate of the M , M states above threshold and their relativeoverlap with the bound state below the threshold.An important component of the present investigation is the tuning of the light u/d quarkmass. Depending on the wave function of the operators we need two different tuning of the u/d mass. For the operators made of heavy-light [ b ¯ l ] mesonic wave functions we find it isnecessary to tune am l to match b ¯ l meson observed mass. Similarly, for light-light diquark[ l l ], where l and l may or may not be equal, in presence of one or more heavy b quarks the am l i is tuned with Λ b . We applied this approach with fair success in bottom baryon [30] andpresently with double bottom tetraquark we attempted the same. In order to understandand explain these two different tuning, we solve the quantum mechanical Hamiltonians of B -meson system, where a single light quark is in the potential of a static bottom quark, andthe Λ b -baryon system, where the two light quarks are in the same field of the static b quark.In this problem b mass is the experimental mass and the light quark mass is treated as aparameter which is tuned to reproduce the experimental B and Λ b masses. We find that22he meson and baryon systems are solved for two different light quark masses which justifiesour need for two different tunings. However, the actual numbers from these two sets of lightquark masses, one from solving the Schr¨odinger equation and the other lattice tuned, cannotbe compared directly due to two different b masses used in these two instances.Once tuned, we find the spectra of the lattice states D, M and M . Naive calculation of hO D O † D i spectrum yields a bound state −
167 (19) MeV measured from the B − B ∗ threshold.On lattice, operators for states having the same quantum numbers can mix and, therefore, itis natural to construct correlation matrices to solve the generalized eigenvalue problem in or-der to obtain the first few lowest lying energies. Besides, the components of the eigenvectorsprovide the relative contribution of the trial states, corresponding to the lattice operators,to the energy eigenstates. They are the coefficient of expansion of the eigenstates whenexpressed in terms of trial states as shown in the equation (31). The GEVP analysis revealsthat tetraquark molecular state just above the threshold by only 17(14) MeV is dominatedby M lattice state while the lowest lying bound state recieves dominant contribution from D along with significantly large contribution from both M and M . From 3 × b tunedcorrelation matrix, we get our final binding energy number for bb ¯ u ¯ d tetraquark system tobe −
186 (15) MeV, where the error is statistical.
VI. ACKNOWLEDGEMENT
The numerical part of this work, involving generation of heavy quark propagators, hasbeen performed at HPC facility in “Kalinga” cluster at NISER funded by Dept. of AtomicEnergy (DAE), Govt. of India. The construction of the correlators and other analysis partof this paper has been carried out in the “Proton” cluster funded by DST-SERB projectnumber SR/S2/HEP-0025/2010. The authors acknowledge useful discussions with RabeetSingh (Banaras Hindu University, India) on Hartree-Fock calculation. One of the authors(PM) thanks DAE for financial support. 23
EFERENCES [1] A.E. Bondar, A. Garmash, A.I. Milstein, R. Mizuk and M.B. Voloshin, Phys. Rev. D ,054010 (2011).[2] A. Bondar et al. (Belle Collaboration), Phys. Rev. Lett. , 122001 (2012).[3] R. Aaij et al. (LHCb Collaboration), Phys. Rev. Lett , 222002 (2014).[4] R.F. Lebed, R.E. Mitchell and E.S. Swanson, Prog. Part. Nucl. Phys. , 143 (2017).[5] A. Esposito, A. Pilloni and A.D. Polosa, Phys. Rep. , 1 (2017).[6] S.L. Olsen, T. Skwarnicki and D. Zieminska, Rev. Mod. Phys. , 015003 (2018).[7] A.V. Manohar and M.B. Wise, Nucl. Phys. B399 , 17 (1993).[8] E.J. Eichten and C. Quigg, Phys. Rev. Lett. , 202002 (2017).[9] Y. Ikeda, B. Charron, S. Aoki, T. Doi, T. Hatsuda, T. Inoue, N. Ishii, K. Murano, H. Nemuraand K. Sasaki, Phys. Lett.
B729 , 85 (2014).[10] Marc Wagner et al. , 012031.[11] M. Padmanath, C.B. Lang and S. Prelovsek, Phys. Rev. D , 034501 (2015).[12] C. Alexandrou, J. Berlin, J. Finkenrath, T. Leontiou and M. Wagner, Phys. Rev. D ,034502 (2020).[13] C. Hughes, E. Eichten and C.T.H. Davies, Phys. Rev. D , 054505 (2018).[14] A. Francis, R.J. Hudspith, R. Lewis and K. Maltman, Phys. Rev. Lett. , 142001 (2017)[15] A. Francis, R.J. Hudspith, R. Lewis and K. Maltman, Phys. Rev. D , 054505 (2019).[16] P. Junnarkar, N. Mathur and M. Padmanath, Phys. Rev. D , 034507 (2019).[17] L. Leskovec, S. Meinel, M. Pflaumer and M. Wagner, Phys. Rev. D , 014503 (2019).[18] P. Bicudo, and M. Wagner, Phys. Rev. D , 114511 (2013).[19] P. Bicudo, K. Cichy, A. Peters, B. Wagenbach and M. Wagner, Phys. Rev. D , 014507(2015).[20] P. Bicudo, J. Scheunert and M. Wagner, Phys. Rev. D , 034502, (2017).[21] P. Bicudo, M. Cardoso, A. Peters, M. Pflaumer and M. Wagner, Phys. Rev. D , 054510(2017).[22] B.A. Thacker and G.P. Lepage, Phys. Rev. D , 196 (1991).
23] G.P. Lepage, L. Magnea, C. Nakhleh, U. Magnea and K. Hornbostel, Phys. Rev. D , 4052(1992).[24] E. Follana, Q. Mason, C. Davies, K. Hornbostel, G.P. Lepage, J. Shigemitsu, H. Trottier andK. Wong, Phys. Rev. D , 054502 (2007).[25] A.X. El-Khadra, A.S. Kronfeld and P.B. Mackenzie, Phys. Rev. D , 3933 (1997).[26] E.B. Gregory et al. (HPQCD Collaboration), Phys. Rev. D , 014506 (2011).[27] N. Kawamoto and J. Smit, Nucl. Phys. B , 100 (1981).[28] J. Jiang, W. Chen and S. Zhu, Phys. Rev. D , 094022 (2017).[29] R.L. Jaffe, Phys. Rep. , 1 (2005).[30] P. Mohanta and S. Basak, Phys. Rev. D , 094503 (2020).[31] A. Bazavov et al. , Rev. Mod. Phys. , 1349 (2010).[32] K. Orginos and D. Toussaint (MILC), Phys. Rev. D , 014501 (1998).[33] K. Orginos, D. Toussaint, and R. L. Sugar (MILC), Phys. Rev. D , 054503 (1999).[34] S. Godfrey and N. Isgur, Phys. Rev. D , 189 (1985).[35] S. Capstick and N. Isgur, Phys. Rev. D , 2809 (1986).[36] folk.uib.no/nfylk/Hartree/lindex.html [37] D. R. Hartree and W. Hartree, Proc. R. Soc. Lond. A150: 933, (1935).[38] V. Fock, Z. Physik 61, 126148 (1930).[39] Marek Karliner and Jonathan L. Rosner, Phys. Rev. Lett. , 202001 (2017)[40] M. Tanabashi et al. (Particle Data Group), Phys. Rev. D98, 030001 (2018) and 2019 update.[41] G. Yang, J. Ping, and J. Segovia, Phys. Rev. D , 014001 (2020), 014001 (2020)