A Parity content of HL interpolating operators

# Tetraquark bound states in the heavy-light heavy-light system

## Abstract

A calculation of the interaction potential of two heavy-light mesons in lattice QCD is used to study the existence of tetraquark bound states. The interaction potential of the tetraquark system is calculated on the lattice with 2+1 flavours of dynamical fermions with lattice interpolating fields constructed using colorwave propagators. These propagators provide a new method for constructing all-to-all spatially smeared the interpolating fields, a technique which allows for a better overlap with the ground state wavefunction as well as reduced statistical noise. Potentials are extracted for 24 distinct channels, and are fit with a phenomenological non-relativistic quark model potential, from which a determination of the existence of bound states is made via numerical solution of the two body radial Schrödinger equation.

1

## I Introduction

The calculation of hadronic forces from first principles allows insight into how interactions of the fundamental quark and gluonic degrees of freedom manifest themselves at the hadronic level. Lattice QCD is an excellent tool for calculating hadronic observables in the low energy regime. Although lattice calculations in euclidean space are not well suited for the study of real-time scattering processes, two methods can be used to extract interaction information from the lattice. One method, developed by Lüscher (1), relates the elastic scattering phase shift of a two particle system in a finite periodic box with the energy levels of the system. An alternate method, used in the present work, extracts the interaction energy as a function of hadron separation. This method is only applicable for systems of hadrons containing more than one heavy quarks which can be treated in the static approximation providing a definite spatial position for the hadrons.

Phenomenologically, two heavy-light meson systems (which we will denote as HLHL) have become interesting in the study of tetraquark bound states (2) (3) (4). It has long been known that the binding of a (with ) system increases with the mass ratio of the heavy to light quark flavours (5), thus and systems are excellent candidates in the search for exotic four quark bound states. In Ref. (4) a distinction was made between two types of tetraquark bound states: molecular, in which the four quarks exhibit a single physical two-meson (singlet-singlet) component, and the more exotic compact bound states. The latter would involve a complicated color space structure in which quark pairs form color vectors which then combine to form a colorless four quark state (4). In spite of this complicated color structure, compact bound states can be interpreted as a mixture of various two meson (color singlet) components (6). The expected features that would characterize a molecular bound state would be a small binding energy and a bound state RMS radius greater than that of the sum of the two particle sizes, i.e.:

 ΔR≡RMS4qRMSM1+RMSM2>1

A compact state, on the other hand, would be more tightly bound and have a smaller RMS radius than the molecular state. In Ref. (7) doubly heavy four quark states were modeled as hadronic molecules interacting via a meson exchange potential. Several of the doubly bottom bound states were found to be deeply bound and spatially compact, making them excellent candidates for tetraquark bound states. It is with these ideas in mind that we may begin to search for the signature of compact bound states on the lattice.

A recent lattice calculation of the HLHL interaction energy (8) has in fact hinted at the possibility of a bound tetraquark state in one channel that exhibits a particularly wide and deep potential well when compared with other channels, although no exhaustive determination of a bound state was undertaken. Our work goes beyond this presenting a quantitative determination of a bound state energy in the HLHL system from a lattice calculation.

An inherent difficulty in making comparisons between theoretical models and lattice calculations performed in the static limit stems from the omission of the heavy quark spin in the static limit. As , the integer valued () angular momemtum eigenstates of a single heavy light meson map onto a single static limit eigenstate with . The energies of the non-static angular momentum eigenstates also converge to a single energy corresponding to the eigenstate. Although the two spaces map onto each other, there is not a simple one to one correspondence between static limit eigenstates and their non-static counterparts, and care must be taken in making identifications between the two spaces. Previous lattice studies of HLHL interaction energy ((9), (10) for example) performed in the quenched approximation and included uncontrolled systematic errors because of this. Recently dynamical quarks have been used to calculate the HLHL interaction energy using a complete set of quantum numbers which exploits the full set of symmetries of the HLHL system (11).

With our choice of quantum numbers (presented in section II) we are able to draw a connection between the quantum numbers and the qualitative behavior of the states. Additionally, by way of symmetry arguments, we are able to relate our static-limit states to non-static angular momentum eigenstates.

## Ii Background

### ii.1 Heavy-Light states

The quark model view of a heavy-light meson is of a heavy anti-quark coupled to a light quark . The relevant quantum numbers to describe such a state are total angular momentum and its projection along some axis (here arbitrarily chosen to be ) , and the parity as well as the relevant flavor quantum numbers. For our interests, we choose and . Therefore all states then have bottomness , and are otherwise classified by total isospin and the third component of isospin . Throughout this work, we make the assumption that we fit our correlation functions with a sufficiently large such that contributions from excited states have died out and we extract only the ground state energy. Furthermore, we assume that states with non-zero orbital angular momentum are at sufficiently high energies as to have a negligible contribution to the ground state energies which we extract. We are then free to speak of the spin and total angular momentum interchangeably.

In heavy quark effective theory, spin dependent contributions enter into the heavy quark action at order , and in the static limit () the heavy quark acts as a static color source. This means that the spin of the HL meson comes only from the light degrees of freedom. Because of this, the physical HL meson states with become degenerate in the static limit, with both represented by a single state. The relevant angular momentum classification is then . With the above flavor assignments, the lowest energy excitations of the B spectrum with (coupling to the static B) are and , and for (coupling to the static ), the ground state (neglecting excited states).

### ii.2 Heavy-Light Heavy-Light states

When constructing states with a pair of HL mesons, care must be taken in determining a relevant set of quantum numbers that fully exploit the symmetries of the problem. The flavor quantum numbers for a Heavy-Light Heavy-Light (HLHL) system are straightforward, and for a (with ) there are two isospin combinations, an isospin triplet with and an singlet. For a HLHL pair separated by a vector the rotational symmetry is broken to rotations around the separation axis. Total angular momentum is therefore no longer a conserved quantity, though its projection along the axis of separation (arbitrarily taken to be ) is still conserved. The system will also be symmetric or antisymmetric under parity as well as reflections through a plane containing the separation axis, which we shall call . This last transformation can be accomplished by a parity transformation followed by a rotation of about an axis perpendicular to the reflection plane. States with are not invariant under this transformation (being mapped onto each other), but their average is an eigenstate of . Lastly we choose to classify HLHL states by intrinsic parity , defined to be the product of the intrinsic parities of the two light quarks, and (full) parity , defined as the product of the intrinsic parity transformation and coordinate inversion of the two particle spatial wavefunction. We will use both parity quantum numbers in our classification of states.

## Iii Methodology

### iii.1 HL and HLHL interpolating fields

A general interpolating operator coupling to a single heavy-light state is given by:

 OHL(→x)=¯Q(→x)Γq(→x) (1)

with chosen to achieve the desired angular momentum and parity quantum numbers. For pseudoscalar HL states, (with ), corresponding to a particle in the static limit with , which we will refer to simply as . meson states with correspond to a state with , which we shall refer to as . We make the choice for and for . As it will be useful in the analysis of HLHL states, it should be noted that for these choices of , correlation functions constructed from interpolating fields will consist of only upper (positive parity) components in the Dirac basis of the light quarks while those constructed from will consist of only lower (negative parity) components. This is explicitly shown in Appendix A. The states are classified by the additional flavor quantum numbers for .

For HLHL states, we want to create states with definite and displacement at the source and sink. To do this, we want to couple only our light quarks in spinor space to specify the quantum numbers of the state while allowing the heavy quarks to act only as color sources. Our general HLHL operator is then given by:

 Missing or unrecognized delimiter for \right (2)

where the light quark wavefunctions are combined in such a way as to achieve the set of quantum numbers of the system. The explicit construction of these wavefunctions is described in Appendix B. For simplicity we restrict ourselves to identical source and sink interpolating fields neglecting any cross correlators between states. Isospin is a good quantum number on the flavor lattices with which we work, and we choose our interpolating fields to be isospin eigenstates with and . At large spatial separations, we expect the energy of the four quark state to asymptotically approach the energy of it’s dominant two meson component2. States with will tend towards the energy of a combination at large spatial separations. There are two possible combinations of the light quark parities that yield : . In light of the above discussion of parity content of single HL states, we project our interpolating fields to contain only negative or positive parity spinor components and retain these as distinct interpolating fields. The expectation is that interpolating fields constructed from lower spinor components will exhibit a significantly higher ground state energy in relation to those constructed from upper components. The reason for this is that the interpolating field are constructed by the product of two meson interpolating fields, thus should exhibit an asymptotic energy (as ) near twice that of the single energy. Similarly the interpolating field is constructed from the product of two meson interpolating fields tending asymptotically as towards a ground state energy of twice that of a single meson. We differentiate all interpolating fields by their dominant asymptotic content in the tabulation of interpolating fields in Table 1.

## Iv Details of the lattice calculation

We work with colorwave propagators (described below) calculated on anisotropic () lattices generated by the Hadron Spectrum Collaboration (12) with a pion mass of roughly 380 MeV. The fermion action used was the clover Wilson action with stout link smearing, not smeared in the temporal direction. The gauge action was Symanzik tree level tadpole-improved without a rectangle in the temporal direction, preserving temporal ultra-locality. The spatial and temporal lattice spacings for these lattices are fm and fm. The pion mass on this ensemble is 0.0681(4) in temporal lattice units. The Chroma Software package for Lattice QCD (13) was used to generate both colorwave and heavy propagators. The calculation of the HL and HLHL energies was performed using 305 gauge field configurations with eight sources spaced evenly in the temporal direction. Ground state energies were extracted using single exponential correlated fits, with an appropriate determined from the quality of the fit.

### iv.1 Colorwave Formalism

#### Two quark states

Consider a general operator for a two quark mesonic state:

 O(→x)=¯q1(→x)Γq2(→x) (3)

where we assume for simplicity that the two quarks have different flavors. We seek to calculate the correlation function with localized interpolating fields: (averaged over spatial source and sink locations to increase statistics)

 C(t,t0) = ∑x,y⟨O(y)O†(x)⟩ (4) = ∑x∑ytr(S1(x,t0|y,t)ΓS2(y,t|x,t0)Γ)

Following the methodology presented in (14), we now consider any complete set of orthonormal states which satisfy:

 Missing or unrecognized delimiter for \right (5)

By inserting the completeness relation of eq. 5 twice into the two point function of eq. 4:

 C(t,t0)= ∑x,x′∑y,y′⟨S1(x,t0|y′,t)δ(y−y′)ΓS2(y,t|x′,0)δ(x−x′)Γ⟩ = ∑x,x′∑y,y′⟨S1(x,t0|y′,t)∑iϕ∗i(y)ϕi(y′)ΓS2(y,t|x′,t0)∑jϕ∗j(x)ϕj(x′)Γ⟩ = ∑i,jSj,i1(t0,t)ΓSi,j2(t,t0)Γ (6)

where we have defined:

 Si,j(t,t0)≡∑x,yϕ∗i(y)S(y,t;x,t0)ϕj(x) (7)

A convenient choice for the is a plane wave basis: . The delta functions here operate on color and spin. With this choice of basis, we define to be colorwave propagators. The use of these propagators allows us to implement spatial smearing at the source and sink of our correlation functions. In the limit where all momenta are summed over in equation 7, all to all point-point propagators are recovered. However, introducing a maximum momentum cutoff we are able to introduce and control the effective amount of spatial smearing3. The effect of restricting the plane wave basis to (summing over a momentum space volume) is illustrated in Fig. 1 where effective masses for single HL meson correlation functions4 are presented. It’s evident that the noise of the signal decreases by increasing the momentum space cutoff (as this increases the statistics contributing to the correlation function).

Each effective mass plateau appears to begin at roughly the same point independent of , and thus a common fit range of was chosen for all values of . In Fig. 1 we can see that as increases the overlap with excited states drops resulting lower values for the effective mass at earlier times. This indicates that a small radial smearing of the quarks field results interpolating fields that have better overlap with the ground state of the system. Such behavior is likely due to the fact that the a non-relativistic HL meson in the static limit is a highly localized object whose wavefunction is confined to a small spatial region.

In light of this behavior and in order to reduce computational cost associated with increasing the momentum cutoff, a value of was chosen for calculations of the HLHL system.

#### HLHL States

We begin with a correlation function for two heavy-light mesons separated by as described above:

 CHLHL(t,→r)= ∑x⟨OHLHL(→x,→r,t)O†HLHL(→x,→r,t0)⟩ (8) = ∑x⟨¯Q(→x,t)¯Q(→x+→r,t)q(→x,t)q(→x+→r,t)¯q(→x+→r,t0)¯q(→x,t0)Q(→x+→r,t0)Q(→x,t0)⟩

Each heavy quark source can only be contracted with the sink at the same spatial location, and upon contraction we work only with the Wilson line portion of the heavy quark propagator, as we want the quantum numbers of the system to be determined entirely by the light degrees of freedom. There are two possible light quark contractions, one where the light quarks contract with source and sink at the same spatial location (direct), and one where the light quarks contract at the other spatial location (crossed). Performing these contractions, we have (omitting the overall color trace):

 CHLHL(t,→r)= ∑xγ5W†(→x;t,t0)γ5γ5W†(→x+→r;t,t0)γ5 ×trd[S(→x+→r,t;→x+→r,t0)S(→x,t;→x,t0)−S(→x+→r,t;→x,t0)S(→x,t;→x+→r,t0)] (9)

Here, denotes the trace over Dirac space spinor indices and is the Wilson line

 W(→x;t,t0)=t∏t′=t0U†4(→x,t′) (10)

We now introduce our partially fourier transformed light quark propagators as:

 S(x′1,t;x1,t0)=∑p′1,p1eip′1x′1S(p′1,t;p1,t0)e−ip1x1 (11)

where sums over momenta have been restricted to as described in the previous section.

Using this, the above correlator can be rewritten as:

 CHLHL(t,→r)=∑p1p′1p2p′2∑x γ5W†(→x;t,t0)γ5γ5W†(→x+→r;t,t0)γ5×ei(p′1−p1+p′2−p2)xei(p′2−p2)r ×[S(p′2,t;p2,t0)S(p′1,t;p1,t0)−S(p′2,t;p1,t0)S(p′1,t;p2,t0)] (12)

Defining

 D(→r,t,t0,ω)≡∑xγ5W†(→x;t,t0)γ5γ5W†(→x+→r;t0,t)γ5ei(ω)x (13)

with , our the final form of our HLHL correlation function becomes:

 CHLHL(t,→r)=∑p1p′1p2p′2 D(→r,t,t0,ω)×ei(p′2−p2)r ×[S(p′2,t;p2,t0)S(p′1,t;p1,t0)−S(p′2,t;p1,t0)S(p′1,t;p2,t0)] (14)

With this method, we calculate the costly first using a parallel code (parallelization over space time) and then perform the far less expensive contractions with the colorwave propagators for our complete operator basis on a scalar workstation class machine.

## V HLHL results

For we have 24 unique HLHL corresponding to the operators enumerated in Table 1. Each potential curve is calculated by taking the jackknife difference between the energy of the HLHL state for various and the energy of the expected two meson asymptotic state:

 V(→r)=EHLHL(→r)−EB(1)−EB(1) (15)

The statistical uncertainty for each point is determined from jackknife statistical analysis. The systematic uncertainties are determined by adjusting the chosen fit range by one time slice in each direction and averaging the observed deviations in the energy. The systematic uncertainty for both and are determined independently and then added in quadrature to determine the systematic uncertainty on .

We find three different asymptotic values for the various states as illustrated in Fig. [2]. The lowest lying asymptotic value corresponds to states with a positive intrinsic parity with all spin components in the correlation function projected to the upper spin components, while the highest asymptotic value corresponds to states with positive intrinsic parity and all spins projected to the lower components. This asymptotic behavior is in line with our expectation that the spin projection of our positive intrinsic parity operators helps to increase the coupling to the lower energy state or the higher energy state. The energy difference between the highest and lowest asymptotic values is roughly twice the energy difference between the single HL and states, indicating that they are both tending asymptotically towards their expected two meson asymptotic energies at long distances. The slight overshoot of the highest asymptotic state beyond it’s expected value of twice the energy for fm may be indicative of contamination from mixing of the HL with a state. All states exhibit an asymptotic tendency towards the sum of the single HL and energies as expected.

As the states with the lowest asymptotic energy values trend most cleanly towards their expected asymptotic value (indicating the least contamination from excited states), we will focus mainly on these states which we present in Fig. [3].

Several aspects of these potential curves should be noted: First, we find that the product of exchange parity and intrinsic parity , which is the symmetry of the two meson spatial wavefunction under spatial inversion, directly corresponds to the attractiveness or repulsiveness of the state. This is in agreement with (8). Second, the exhibits a significantly deeper and wider potential well when compared with the two other attractive channels. This qualitative difference was acknowledged in (8), and the quantum numbers of this channel are consistent with a bound state predicted in a phenomenological model in (4).

## Vi Bound States

As the HLHL system has been predicted to be an excellent candidate for bound tetraquark states, we seek a quantitative method for extracting such a bound state (if one exists) from our lattice calculation. Our method is as follows: We fit our lattice potential to a phenomenological quark model potential as described in (15). We make the choice to focus on the channel, as previous work has hinted at the possibility of a bound state here. As a control, we also perform the fit for the channel as well. In our fit, we neglect the points as the finite value of the potential at is a lattice artifact stemming from the ultraviolet cut off introduced by the lattice discretization, leaving us with 7 data points for each potential curve, and two free parameters from the fit model. The model with the extracted fit parameters is then taken to be the interaction potential between two B mesons in the continuum limit. The two body (one-dimensional) Schrodinger equation is then solved numerically with this interaction potential to determine the existence of any negative energy (bound) states. It should be noted here that the solutions to the Schrodinger equation will converge to their continuum values as the continuum limit of the lattice calculation is taken. As we have only a single lattice spacing available to work with this continuum extrapolation is not an option, and it should be understood that the results presented in this section are at finite lattice spacing.

### vi.1 Potential Model

We have limited our displacements fm, therefore long range effective interactions due to meson exchange do not provide a good description of the HLHL system. In reference (15), a quark model picture of a two meson interaction was used to derive an interaction potential for the HLHL system, which included color coulomb, spin-spin, linear confinement interactions. Details of the derivation of the potential model can be found in the aforementioned reference, and we will only highlight several modifications we make when fitting this potential model to our numerical results. The quark model HLHL potential has the form:

 Missing or unrecognized delimiter for \right (16)

with:

 Vcc(αs,β,r) =−4αs9r[1+(2π)1/2βr−4Erf(βr2)]e−β2r2/2 (17) Vss(αs,β,¯m,r) =227(2π)1/2αsβ3¯m2e−β2r2/2 (18) Vlc(b,β,r) =b3β[βre−β2r2/2+2(2π)1/2e−β2r2/2 −(βr+2βr)Erf(βr2)e−β2r2/2−2π1/2e−3β2r2/4] (19)

Here, is the strong coupling constant, is the spatial width of the quark model single HL meson wavefunction, is the mass of the light quark in the scheme, and is the QCD string tension. The coefficients and , which contain the spin information of the HLHL state, are defined as matrix elements between initial (unprimed) and final (primed) two meson states and will be discussed further below. It should be noted that the above potential model acquires an overall minus sign if the isospin wavefunction of the two meson state is antisymmetric. Additionally, the potential is a function of and not , as any tensor interaction terms are neglected in this model.

### vi.2 Fit Model

When applying the above model to our lattice data, we must make several modifications to the above quark model potential. Due to the use of periodic boundary conditions in the calculation, interactions with image “charges” lying past the boundary must be accounted for. We must also consider the possibility that there will be long range meson exchange interactions that were neglected in our choice of potential model. To account for these long range interactions, we extend the original model by adding a simple Yukawa like term for one pion exchange:

 VYuk(r)=VBBDS(r)+ge−mπrr (20)

Here we take to be the mass of the pion on the gauge field configurations used in the calculation ( GeV). The parameter is discussed below.

In principle, interactions with each of the infinitely many image charges contribute to the potential and must be included. In practice however, we may restrict ourselves to contributions where the image of the first meson is ( fm) away from the second and vice versa. This approximation is valid as the contribution of these truncated images (at separations of ) to the potential (with the choice of parameters outlined below) is . With the inclusion of the image charges our potential model then becomes:

 VYukIm=VYuk(r)+2∑ri<3L/2VYuk(ri) (21)

The addition of these image charges modify the potential at long distance as illustrated in Fig. 4

The final modification made to the potential model is a modification of the spin dependent coefficients and . The original presentation of this phenomenological potential model in Ref. (15) sought to provide a comparison with the lattice calculations of the time, which had an incomplete classification of the HLHL states in terms of the total isospin and spin of the system, while also maintaining a connection with the physical B meson states. Because of this, classification of the various potentials was made in terms of the physical and (first angular excitation of the meson) with respect to the quantum numbers and .

The difference in angular momentum spaces of the non-static and static limit prevents a direct interpretation of the lattice data from the present work in terms of physical and states, and our classification of states makes it difficult to reconcile the previous classification with ours. We therefore choose to recalculate the spin dependent coefficients of the potential model relevant for the static limit BB system we study on the lattice, the results of which are presented in Table 2 (For details of the calculation, see Appendix C). The previous determination of these coefficients for the HLHL system included spin degrees of freedom for the heavy quarks in the two meson states allowing for better classification of the potential in terms of non-static limit states. We choose to neglect the spin degrees of freedom of the heavy quarks in our determination, effectively fully implementing the static limit for the the potential model. Thus the spin degrees of freedom of our two meson kets are just those of the spin of the light degrees of freedom of our HLHL state. The evaluation of these coefficients however requires knowledge of the total angular momentum of the two meson state, a point that has been neglected until now. As we seek to fit the and states, we need to determine if these particular states are in a symmetric angular momentum triplet, or an antisymmetric angular momentum singlet. In order to make this identification, we must rely on the overall symmetries of the state in question. We know that the parity of a given state is the product of the intrinsic parity and the symmetry of the spatial wavefunction. From this relationship, and with knowledge of the symmetry of the isospin spatial wavefunction, we can infer the symmetry of the angular momentum wavefunction:

 SymJ=(−)(SymI)(Pi)(P), (22)

where and the symmetries of the angular momentum and isospin wavefunctions. The overall negative sign appears from exchanging fermions in the parity operation. Using this we are able to identify the channel with as a state, and the channel with as a state. The spin dependent coefficients can then be recalculated for our states and are shown in Table [2].

### vi.3 Fitting Procedure and Bound State Determination

In fitting the potential model of eq. 20 to our lattice data, we use two free fit parameters: and and take the remaining parameters , and to be , , and respectively as in Ref. (15). A fit is performed for each of 305 jackknife ensembles, allowing for an accurate way to estimate the error on the extracted fit parameters, shown in Table [2]. As we are ultimately interested in the energy levels allowed by the potential model, and not the model parameters themselves, we will only briefly comment on the fit parameters. It is immediately obvious that is not well determined for the channel. It’s also interesting that the fit parameter is significantly smaller for the channel, indicating a much narrower spatial distribution of the two meson wavefunction.

Once the fit parameters have been extracted they are then inserted into the two body radial Schrodinger equation to determine if any bound states exist. As we are restricting ourselves to states, the two body Schrodinger equation can be written as:

 [−ℏ22¯md2dr2+VYuk(r)]u(r)=Eu(r) (23)

where is the reduced mass of a two B meson system (with the single meson mass taken from the Particle Data Group (16)), and is the potential model presented in the preceding section excluding the image terms.

Eq. 23 is then solved numerically as an eigenvalue problem with a spatial discretization of 0.01 fm and a spatial cutoff of 10 fm (corresponding to a sphere with fm), and the boundary condition that . This spatial volume provides ample space for the potential to decay to zero. The eigenvalue spectrum is then analyzed for each of the two states discussed above. While the channel exhibits a near continuum of positive eigenvalues (discrete only because of the numerical solution method), the channel does admit a single bound state with energy MeV (with the uncertainty determined by carrying through the jackknife analysis from the fit parameters and solving eq. 23 for each of the 305 sets). Aside from the binding energy, we can also calculate the RMS radius for the two meson wavefunction from the wavefunctions above:

 rRMS≡⟨r2⟩1/2=[∑ir2i|u(ri)|2∑i|u(ri)|2]1/2 (24)

For the bound state wavefunction , we find an RMS radius of 0.383(6) fm, the error again estimated by jackknife analysis.

Although no previous calculation of the binding energy in this particular static-limit channel exists (lattice or otherwise), Ref. (4) does quote binding energies and RMS radii for a doubly bottom channel which is consistent (in the static limit) with the quantum numbers of our static limit channel. This reference uses two different potential models to calculate binding energies: the constituent quark cluster model CQC and the the Bhaduri-Cohler-Nogami or BCN model. The BCN model includes the same interactions as those used in Ref. (15) to derive the potential used to fit our lattice results (namely, color coulomb, linear confinement and spin-spin). Furthermore, the BCN parameters corresponding to string tension , strong coulpling , and constituent quark mass used in (4) are very similar to those used in our potential model (compare our to ). These binding energies should provide a relevent point of comparison for our results provided our lattice discretization errors have minor effects on the extracted potential model fit parameters. In comparison, we find our values for the binding energy and RMS radius to be consistent with the values quoted in (4) from the BCN model , providing a good cross check that our lattice calculation has identified a bound state in the static limit channel. The fact that the bound state identified in that work has an RMS radius that is smaller than the sum of the individual mesonic RMS radii is indicative of the compact nature of that bound state. Additionally, as illustrated in Ref. (11) (see eqns. 4), the static limit HLHL tetraquark state can be written as a linear combination of products of two single meson wavefunctions in different spin states. This is consistent with the idea that although the compact tetraquark state may have a complicated color space structure composed of color vectors, this state can always be decomposed into a linear combination of products of two single meson wavefunctions.

## Vii Conclusions

We have computed using lattice QCD the interaction potential between two b-meson states in the limit of static b quarks. With this lattice potential parametrized with a functional form motivated by the quark model description of the two b-meson interaction, we have determined the bound state energies in the heavy-light-heavy-light (HLHL) tetraquark system. To perform this study we introduced colorwave propagators for calculating meson correlation functions and extended the formalism to the HLHL system in order to provide a novel way for an efficient calculation of HLHL correlation functions for several channels. The effect of limiting the colorwave plane wave basis on the ground state overlap of single HL correlation functions was explored, and a choice for the momentum cutoff was made to optimize the quality of the signal versus the computational cost. For a single HL meson, results indicate that a more localized interpolating field has a better overlap on the ground state, suggesting the compact nature of the HL meson.

HLHL potentials were calculated for 24 distinct channels, exhibiting three distinct asymptotic values as corresponding to the different ways and mesons can be combined. The tendency of the HLHL energy to overshoot the expected asymptotic value of and may be due to contamination from excited states and the possibility of mixing with a state. It was determined that the attractiveness or repulsiveness of the HLHL potential corresponds directly to the symmetry of the two meson spatial wavefunction under spatial inversion, in agreement with Ref. (8). The asymptotic behavior of the various HLHL states was shown to be dependent on the intrinsic parity of the state. While the states have only one asymptotic value (corresponding to a single two meson component), the channels have two asymptotic values corresponding to both and two meson components. By examining the construction of single HL correlation functions, it was determined that we could increase overlap with the and two meson wavefunctions by projecting the correlation functions to include only positive or negative parity components of the Dirac basis quark spinors.

The existence of bound states was then explored for the channel as it exhibited a wider and deeper potential when compared with the other attractive potentials. Analysis was also carried out for the for the purposes of comparison. A modified version of the potential model described in Ref. (15) was used to fit the lattice data, and two fit parameters (the gaussian width of the HL meson wavefunction) and (the Yukawa interaction constant) were extracted from, the fit. Inserting the potential with the extracted fit parameters into the two body Schrödinger equation, we then solved numerically for the eigenvalues of the hamiltonian, searching for any negative energy eigenstates. A single negative energy bound state was found in the channel, with an energy of MeV and RMS radius fm. These results were found to be consistent with results presented in Ref. (4) for the state (which maps onto our channel in the static limit). The errors quoted on these results are statistical only. One needs to account for several systematic errors such as corrections ( the b quark mass), lattice spacing effects as well as dependence on the light quark mass. 5

###### Acknowledgements.
We would like to thank W. Detmold, S. Meinel, and E. Mastropas for useful discussions. This work was supported in part by DOE grants DE-AC05-06OR23177 (JSA) and DE-FG02-04ER41302 and DE-FG02-04ER41302 as well as the JSA Jefferson Lab Graduate Fellowship Program.

## Appendix A Parity content of HL interpolating operators

Here we show that correlation functions for our () states are composed entirely of upper(lower) components of the Dirac basis components of the light quark flavors. We begin with a general HL correlation function with arbitrary source and sink operators (neglecting color indices and working in the Dirac basis):

 CHL(t)i,j =∑→x⟨OBi(→x,t)O†Bj(→x,0)⟩ =∑→x⟨¯Q(→x,t)Γiq(→x,t)¯q(→x,0)ΓjQ(→x,0)⟩ =∑→xtr(γ5(SH(→x,t;0))†γ5ΓiSL(→x,t;0)Γj) (25)

is a heavy quark propagator given by:

 SH(→x,t;t0)=(1+γ42)W(→x,t;t0)=P+W(→x,t;t0) (26)

where is a Wilson line from to . Substituting this, we have:

 CHL(t) =∑→xtr(γ5(P+W†(→x,t;0))γ5ΓiSL(→x,t;0)Γj) =∑→xtrc(W†(→x,t;0)trd(ΓiP−ΓjSL(→x,t;0))) (27)

where we have used . For , we project to only the lower components of the Dirac basis light quark propagator, while for we project only to the upper components of the Dirac propagator.

## Appendix B Construction of light quark wavefunctions

To determine two quark wavefunctions in spin and flavor space yielding the quantum numbers , we begin with states of definite :

 Missing or unrecognized delimiter for \bigg (28)

where are the projections of spin and isospin along the z-axis and the intrinsic parities of the light quarks, and the , are the Clebsch-Gordon for angular momentum and isospin. From these operators, we average over states and determine from the quantum numbers and and the spatial symmetry of the operator. It should be noted here that there are two combinations of that contribute to the HLHL states, and we make the decision to keep these as distinct operators.

Linear combinations of the above operators are then taken to produce states of definite exchange parity , the necessary combinations determined by summing over sets of the above operators that map onto each other under with the appropriate weight

 [q1q2]∣∣∣(I,Iz,|Jz|,P⊥,P,Pi)=∑p1,p2WPp1,p2[q1(p1)q2(p2)]∣∣∣(I,Iz,|Jz|,P⊥,Pi) (29)

## Appendix C Determination of spin coefficients for potential model

Here we present our derivation of the spin coefficients and presented in Table 2. In Ref. (17), an interaction potential for two meson states is calculated by including spin-spin, color coulomb and linear confinment interactions in a two quark interaction hamiltonian. By considering these interactions between each of the quark quark pairs in a 4 quark (2 meson) scattering state, transfer matrix elements are calculated and then Fourier transformed to give a corresponding position space potential. In Ref. (15), this method was applied to the HLHL system. When calculating the spin dependent portion of the potential, all but one of the interaction diagrams (referred to as “Transfer 2”) can be neglected because the spin of the heavy quarks is neglected in the static approximation. This diagram includes an insertion of the interaction hamiltonian between the two light quarks, as illustrated in in Fig. 6. The spin dependent contribution of this diagram to the potential can be factorized such that all the dependence enters through two coefficients, which are defined as matrix elements between the initail and final two meson states:

 CI =⟨CD|I|AB⟩ (30) CS⋅S =⟨CiDj|Si⋅Sj|AiBj⟩ (31)

Where here is understood to be the identity operator in spinor space. Upon inspection of the diagram, it’s clear that the matrix element of will not always trivially be unity due to the quark interchange between the initial and final two meson state.

With respect to Fig. 6, these matrix elements as outlined in (17) are defined explicitly as:

 CI =⟨CD|I|AB⟩ =χλCsc,s¯cχλDsd,s¯d[δs¯a,s¯cδs¯b,s¯d]χλAsa,s¯aχλBsb,s¯b (32) CS⋅S =⟨CiDj|Si⋅Sj|AiBj⟩ =χλCsc,s¯cχλDsd,s¯d[δs¯a,s¯cδs¯b,s¯d]χλAsa,s¯aχλBsb,s¯b (33)

For our purposes, we wish to entirely neglect the spin of the heavy quarks in the above matrix elements. Because of this, the Clebsch-Gordan coefficients etc. (relating the spin of the two quark state to the meson state) are all unity. The states between which we wish to calculate these matrix elements are two particle angular momentum eigenstates , of which we are only interested in and . To account the light quark exchange in Fig. 6, we note the following relations:

 |sa,sb⟩∣∣∣J=0,Jz=0=1√2(|↑↓⟩−|↓↑⟩)=−1√2(|↓↑⟩−|↑↓⟩)=−|sc,sd⟩∣∣∣J=0,Jz=0 (34)

and

 |sa,sb⟩∣∣∣J=1,Jz=0=1√2(|↑↓⟩+|↓↑⟩)=1√2(|↓↑⟩+|↑↓⟩)=|sc,sd⟩∣∣∣J=1,Jz=0 (35)

From the above relations, it is easy to calculate the matrix elements of interest for our problem (for the states and ):

 ⟨1,0|c,dI|1,0⟩a,b =⟨1,0|a,bI|1,0⟩a,b=1 (36) ⟨0,