Construction of $bb\bar{u}\bar{d}$ tetraquark states on lattice with NRQCD bottom and HISQ up/down quarks

We construct $bb\bar{u}\bar{d}$ states on lattice using NRQCD action for bottom and HISQ action for the light up/down quarks. The NRQCD-HISQ tetraquark operators are constructed for"bound"$[bb][\bar{u}\bar{d}]$ and"molecular"$[b\bar{u}] [b\bar{d}]$ states. Corresponding to these different operators, two different appropriately tuned light quark masses are needed to obtain the desired spectra. We explain this requirement of different $m_{u/d}$ in the light of relativised quark model involving Hartree-Fock calculation. The mass spectra of double bottom tetraquark states are obtained on MILC $N_f=2+1$ Asqtad lattices at three different lattice spacings. Variational analysis has been carried out to obtain the relative contribution of"bound"and"molecular"states to the energy eigenstates.


I. INTRODUCTION
The multiquark hadronic states other than the mesons and baryons are relatively new entrants particularly in the heavy quark sector. The signature of some of such states containing four or more quarks and/or antiquarks have been observed in experiments [1][2][3][4][5][6].
Such states are characterized by J P C quantum numbers that cannot be arrived at from the quark model. (Such four quark states are popularly referred to as tetraquark, which is used to denote either bound or often both bound and mesonic molecular states. In this paper we use the term tetraquark in the latter sense.) However, heavy hadronic tetraquark states QQll 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 Ql − Ql threshold. Of late, the observations of Z − (4430, 1 + ) of minimal quark content being ccdū [3] have been reported along with the 1 + states like Z b (10610) and Z ′ b (10650), having minimal quark content of four quarks (containing a bb pair) that are a few MeV above the thresholds of B ⋆B (10604.6) and B ⋆B⋆ (10650.2) [1,2]. The proximity of Z b , Z ′ b to the B ⋆B⋆ threshold values perhaps suggest molecular, instead of bound, nature of the states.
Around the same time, lattice QCD has been employed to investigate the bound and/or molecular nature of the heavy tetraquark states, not only to understand the above experimentally observed states but also to identify other possible bound tetraquark states in both 0 + and 1 + channels. In the charm sector, some early lattice studies involve T cc and T cs tetraquark states [9], cccc [10], X(3872) and Y (4140) [11] and more recently D ⋆ s0 (2317) [12]. The bottom sector received intense attention where, instead of B * B or B * B * , relatively simpler BB, BB * systems are studied. The lattice investigations this far involve four bottom bbbb [13] and two bottom tetraquark statesbbl 1 l 2 , where l 1 , l 2 ∈ c, s, u, d, [14][15][16][17]. An important observation of these lattice studies is that the possibility of the existence of bbl 1l2 tetraquark bound states increases with decreasing light quark masses, while they become less bound with decreasing heavy (anti)quark mass.
Besides the usual lattice simulations, the heavy tetraquark systems have also been studied using QCD potential and Born-Oppenheimer approximation [18][19][20][21]. The main idea in this method is to investigate tetraquark states with two heavy (anti)quarks, which wasbb in the study, and two lighter quarks using quantum mechanical Hamiltonian containing screened Coulomb potential. This approach has been used to explain our two different choices of light u/d quark masses for different classes of tetraquark states.
In this work our goal is to construct tetraquark states, having quark content bbl 1l2 in 1 + both below and above B − B * threshold, by a combination of lattice operators and tuning 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 HISQ action [24] for l 1 , l 2 = u/d. Here we also explore through variational/GEVP analysis how the trial states created by our operators contribute to the energy eigenstates.
First, we briefly review the salient features and parameters of both NRQCD and HISQ actions along with the steps involve in combining the relativistic u/d HISQ propagators with the NRQCD b quark propagators in the section II. We have considered two different kind of operators -the local heavy diquark and light antidiquark (often referred as "good diquark" configuration) and molecular meson-meson, we described these constructions in the section III. We collect our spectrum results in section IV that contains subsections on quark mass tuning (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 our results in section V.

II. QUARK ACTIONS
Lattice QCD simulations with quarks require quark mass to be am l ≪ 1, where a is the lattice spacing. In the units of the lattice spacings presently available, the b quark mass is not small i.e. am b ≮ 1. As is generally believed, the typical velocity of a b quark inside a hadron is nonrelativistic v 2 ∼ 0.1 and is much smaller than the bottom mass. This makes NRQCD our action of choice for b quarks on lattice. We have used O(v 6 ) NRQCD action [23], where the Hamiltonian is and O(v 6 ) terms are in δH with coefficients c 1 through c 7 , where ∆ ± and ∆ 2 are discretized symmetric covariant derivative and lattice Laplacian respectively. Both the derivatives are O(a 4 ) improved as are the chromoelectric E and chromomagnetic B fields. The b quark propagator is generated by time evolution of the Hamiltonian H, The tree level value of all the coefficients c 1 , c 2 , c 3 , c 4 , c 5 , c 6 and c 7 is 1. Here n is the factor introduced to ensure numerical stability at small am b , where n > 3/2m 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 momentum directly from the exponential fall-off of their correlation functions. Instead, we calculate the kinetic mass M k of heavy-heavy mesons from its energy-momentum relation, which to O(p 2 ) is [25], where M k is the kinetic mass of the meson which is calculated considering E(p) at different values of lattice momenta p = 2 nπ/L. The b quark mass is tuned from the spin average of kinetic masses of Υ and η b , and matching them with the experimental spin average value, The experimental value to which M bb is tuned to, however, is not 9443 MeV that is obtained from spin averaging Υ (9460 MeV) and η b (9391 MeV) experimental masses, but to an appropriately adjusted value of 9450 MeV [26], which we denote as M mod phys in the equation (6) below. The hadron mass is then obtained from where E latt is the lattice zero momentum energy in MeV, n b is the number of b-quarks in the bottom hadron.
The u/d light quarks comfortably satisfy the criteria am l ≪ 1 and, therefore, we can use a relativistic lattice action. We use HISQ action for the u/d quarks, which is given in [24], Because HISQ action reduces O(α s a 2 ) discretization error found in Asqtad action, it is well suited for u/d (and s) quarks. The parameter ǫ in the coefficient of Naik term can be appropriately 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 propagators do not have any spin structure. The full 4 × 4 spin structure is regained by multiplying the propagators by Kawamoto-Smit multiplicative phase factor [27],

III. TETRAQUARK OPERATORS
In the present paper, we have considered two kinds of tetraquark operators -the local heavy diquark and light antidiquark and molecular meson-meson. The b quark, being nonrelativistic, is expressed in terms of two component field ψ h . We convert it to a four component spinor Q having vanishing two lower components, which help us to combine b field and relativistic four component light quark fields in the usual way. The heavy-light meson operator, that we will make use of in the operator construction, is written as where l(x) stands for the light quark fields,Q = Q † γ 4 and depending on pseudoscalar and vector mesons Γ = γ 5 and γ i respectively.
Because of the vanishing lower components, the states with Q can only be projected to the positive parity states. The local double bottom tetraquark operators that we can construct for bbl 1l2 system are, where l 1 = l 2 and l 1 , l 2 ∈ 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 brackets are color contracted. The diquark-antidiquark 1 + four quark state bbl 1l2 with l 1 = l 2 in (13) can actually be defined in two ways [28], where a, b are color indices. The subscripts Q * andπ in the operator O Q * π are in3 c and 3 c respectively, while Q andπ * in the operator O Qπ * are in 6 c and6 c . But both O Q * π and O Qπ * correspond to the 1 + state. Of these the O Q * π is our desired "bound" tetraquark operator because one-gluon-exchange interaction is attractive for a heavy quark pair in3 c diquark configuration [8] and spin dependent attraction exists for light quark pairs in "good diquark" configuration characterized by color3 c , spin J = 0 and isospin I = 0 or 1/2 [29].
The two terms in O Q * π contribute identically in the final correlator, hence we consider only the first term in the calculation. The generic form of the temporal correlation among the operators at zero momentum is, where X, Y can be any of D, M 1 , M 2 in equations (11), (12) and (13). For example, the explicit forms of the zero momentum correlators, including cross-correlator, when X and Y are M 1 = B * B and D = Q * π , are Above in the equations (17) and (18) has to be rotated to the MILC basis before implementing the equations (16), (17) and (18).
The unitary matrix needed to do this transformation is given by [30],

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 these lattices can be found in [31]. It uses Symanzik-improved Lüscher-Weisz action for the gluons and Asqtad action [32,33] for the sea quarks. The lattices we choose have a fixed ratio of am l /am s = 1/5 with lattice spacings 0.15 fm, 0.12 fm and 0.09 fm and they correspond to the same physical volume. We have not determined the lattice spacings independently but use those given in [31]. In the Table I we listed the ensembles used in this work.

A. Quark mass tuning
For bbūd mass calculation, we need nonperturbative tuning of both m b and m u/d . With the help of equation (5), the tuning of m b has been carried out by calculating the spin average Υ and η b kinetic masses and comparing the same with the spin average and suitably adjusted experimental Υ and η b masses as discussed in section II. The tuned bare am b quark masses for lattices used in this work are given in Table II. But the tuning of am u/d is rather tricky and the way to achieve this numerically has been discussed in the study of bottom baryon spectra in [30]. In the present paper, we try to understand the light quark tuning in more details with the help of relativised quark model [34,35] , and therefore, we use experimental Λ b mass 5620 MeV to tune the bare am u/d . For the operators O M 1 /M 2 , the diquark part is formed between heavy quark and light In such case we tend to use B-meson mass 5279 MeV to tune the am u/d . The result of u/d quark mass tuning that are made use of in this work is given in the Table II. The relativised quark model [34,35] helps us to numerically calculate the masses of B meson and Λ b baryon using the light (anti)quark mass as parameter. The molecular tetraquark state 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 solve the problem by considering the radial part of the Schrödinger equation numerically using suitably modified Herman-Skillman code [36].
Here U(r) = rψ(r) and the potential V (r) is given by The B meson mass M B is, therefore, determined from the energy eigenvalue E, where r 12 is the relative distance between two light quarks "orbiting" the two heavy quarks and 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.6 [34].
To solve the Hamiltonian (23), we consider the trial wave function, which is spacesymmetric and spin-antisymmetric, in terms of Slater determinant 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 be written as where, we have used In contrast to the helium atom, the presence of linear r-terms in the Hamiltonian leads to additional exchange-energy terms in the calculation. With these linear r-terms in, the

Hartree-Fock equation becomes
We solve for E in equation (26) iteratively and, eventually, the Λ b mass is calculated from The PDG value of Λ b (5620) is obtained by setting the m u/d to 0.157 GeV. In Table III, we compare the nonperturbatively tuned m u/d on our lattices with those obtained by solving the equations (20) and (26). The bare lattice light quark masses cannot be directly compared to the parameter m u/d in these equations mainly because of the use of 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   [14][15][16][17] are found to be below this threshold as is obvious from the Fig. 3.

Mass in MeV
To this end, in Fig. 4 we plot the effective masses of these two states obtained at different lattice spacings together with the lattice thresholds for easy comparison. The colored bands represent fitted am eff values. The superscripts Λ b and B denote the light quark tuning.  In the Table IV, we present our results of the tetraquark states corresponding to the operators given in the expressions (11,12,13). We use two-exponential uncorrelated fit to    (7) 10655 (8) (2) 5284(3) 5292(8) gives us an estimate of the binding energy, In the last column, we calculate our lattice average of ∆M X in MeV and compare with some of the previous lattice results. To our knowledge, the binding energies of the M 1 , M 2 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 excited states, 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 we can use their result as a reference. As is generally understood, these operators are expected to have overlap with the desired ground and excited states of the tetraquark system of our interest. The variational analysis can be performed to determine the eigenvalues and the eigenvectors from the mixed states formed by lattice operators. This is typically achieved by constructing a correlation matrix involving the lattice operators O X , where X, Y can either be all or any two combinations of D, M 1 , M 2 in the expressions (11 -13). The terms n|O † X |Ω are the coefficients of expansion of the trial states O † X |Ω , where |Ω is the vacuum state, and written in terms of the energy eigenstates |n as, Presently, we are interested in expressing the energy eigenstates in terms of the trial states to understand the contribution of each to the former. If we confine ourselves to the first few energy eigenstates, we can write The v X m are equivalent to the eigenvector components obtained by solving a GEVP w.r.t a suitably chosen reference time t 0 [12], The eigenvalues λ m (t) are directly related to the energy of the m-th state, ground and the first few excited states, of our system through the relation The component of eigenvectors v m (t, t 0 ) gives information about the relative overlap of the three local operators to the m-th eigenstate. The eigenvectors v m 's are normalized to 1.
To determine the parameter t 0 , we solve the GEVP and found that the ground and excited state energies are almost independent for t 0 = 3, 5, 7, 9 as demonstrated in the Fig. 5. In this plot, we showed our results for B-tuned am u/d for all the operators O X but the results are similar with Λ b tuning and, hence, not shown. We chose t 0 = 5 for our calculations.  The GEVP analysis has been carried out in two steps because of differences in the tuning of am u/d for the "molecular" states M 1 , M 2 and "good" diquark state D. In the Table VI In the Fig. 6 However, the first excitation |1 , for which our data is rather noisy to reliably estimate ∆M, the |M 1 and |M 2 states appear to have comparable contribution and are broadly distributed over different time slices and vary significantly over configurations. This is evident from the histogram plot in the Fig. 7. This may have a bearing with the fact that above the threshold, the Z b tetraquark can couple to multiple decay channels resulting in a broad spectrum. that define the energy eigenstate Including  |0 , E 0 = 1.784 (12)  despite the excited state energy E 1 is below the threshold. Our data for |2 is too noisy to extract much information. Based on this Λ b tuned 3 × 3 GEVP analysis, the binding energy for the ground state obtained is −186 (15) MeV.

V. SUMMARY
In this work we have attempted to study two possible bbūd tetraquark states -one which is bound and where most other lattice results are centered and, the other which is reported in PDG as Z b and Z ′ b . The experimentally observed states are believed to contain a bb pair but ours is bb pair which is considered as theoretically simpler. However, the possible molecular nature of Z b and Z ′ b suggests that our molecular states should have similar masses. For the bottom quarks we have used NRQCD action while HISQ action for the u/d quarks. This NRQCD-HISQ combination has been employed earlier in [26] for bottom meson and recently in [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 1 , M 2 with the expectation that they will contribute to the states above B − B * threshold. There are not many lattice results on the states above threshold possibly because of the complication that they can couple to multiple decay channels besides B * B and B * B * . Our motivation here is to obtain a tentative estimate of the M 1 , M 2 states above threshold and their relative overlap with the bound state below the threshold.
An important component of the present investigation is the tuning of the light u/d quark mass. 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 [bl] mesonic wave functions we find it is necessary to tune am l to match bl meson observed mass. Similarly, for light-light diquark [l 1 l 2 ], where l 1 and l 2 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] and presently with double bottom tetraquark we attempted the same. In order to understand and 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, and the Λ 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 a parameter which is tuned to reproduce the experimental B and Λ b masses. We find that the meson and baryon systems are solved for two different light quark masses which justifies our need for two different tunings. However, the actual numbers from these two sets of light quark masses, one from solving the Schrödinger equation and the other lattice tuned, cannot be 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 1 and M 2 . Naive calculation of O D O † D 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, it is natural to construct correlation matrices to solve the generalized eigenvalue problem in order to obtain the first few lowest lying energies. Besides, the components of the eigenvectors provide 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 when expressed in terms of trial states as shown in the equation (31). The GEVP analysis reveals that tetraquark molecular state just above the threshold by only 17 (14)