Beauty-full Tetraquarks

In this article we present a calculation of the $bb\bar{b}\bar{b}$ tetraquark ground-state energy using a diffusion Monte Carlo method to solve the non-relativistic many-body system. The static potential for the four quark system is chosen to align with the flux-tube picture of QCD. Using this approach, we find that the $0^{++}$ state has a mass of $18.69\pm 0.03~\mbox{GeV}$, which is around 100 MeV below twice the $\eta_b$ mass. This bound state can behave as a four-lepton resonance via its decay to $\Upsilon(1S) \Upsilon(1S)^* \rightarrow \ell^+ \ell^- \ell^+ \ell^-$.


I. INTRODUCTION
In the seminal work of Ref. [1], Applequist and Politzer interpreted heavy quark bound states as positronium-like atoms subject to non-relativistic quantum mechanics calculations. The spectroscopy of quarkonia can then be well understood by solving Schrödinger's equation based on the static potential between two heavy quarks [2][3][4] mediated by the asymptotically-free quantum chromodynamic (QCD) interactions. Ever since, two-body heavy quark systems have been used to understand the long distance behavior of QCD.
Multi-quark states were first proposed in 1964 by Gell-Mann as an explanation for the observed spectrum of mesons and baryons [5]. This now universally-accepted picture of mesons (baryons) as two-(three-)quark states has been hugely successful, and suggests that the mass of ordinary matter can be explained by the binding energy between quarks. More recently, enormous progress has been made both theoretically and experimentally in our understanding of four-quark states containing at least one light quark [6,7]. On the other hand, four-quark states containing only heavy quarks have not been directly confirmed by experimental searches. With new data from CMS, ATLAS, and LHCb at the LHC, a multiquark state containing only bottom and/or charm quarks is very likely to be tested. From the theoretical side, such heavy quark states provide a unique environment to examine the non-relativistic QCD effective potential of many-body systems. In this paper, we concentrate on a potential tetraquark state comprised of two b and twō b quarks.
In QED, the equivalent system is the di-positronium molecule (Ps 2 ), first postulated by Wheeler in 1946 [8]. Using the variational method, the binding energy against dissociation into two positronium atoms is calculated to be 0.435 eV [9][10][11], which is around 3.2% of the binding energy of two positronium atoms. However, it wasn't until 2007 that the use of positron traps and accumulators led to the experimental confirmation of the Ps 2 molecular state [12]. For the neutral positronium atoms, the electric dipole-dipole interaction can be used to generate the splitting between the ground and excited states. After the standard quantum mechanical perturbative calculation, the additional binding energy can be interpreted in terms of the R −6 London-Van der Waals force [13].
For QCD, there is no equivalent chromo-electric dipole interaction for a color-neutral meson, but one still has the transition-dipole interaction from color-neutral state to color-octet state. Depending on the relation between the inter-meson distance and intra-meson binding energy, second-order perturbative calculations reveal a similar R −6 London-Van der Waals force [14,15] or an R −7 Casimir-Polder force [16][17][18]. As emphasized in Ref. [14], the Van der Waals force arises at leading order in α s and depends only on the geometric ratio of inter-meson and intra-meson binding energies. In the heavy quark limit where the QCD effective potential is Coulomb-like, one can use the Van der Waals force to calculate the additional inter-meson binding energy and obtain a ratio of O(1%) of the total intra-meson binding energy, similar to the Ps 2 case. For bottom quarks with finite mass, one should also include a long-range linear contribution to the potential [19].
Several methods have been proposed in the literature to calculate the energies of tetraquark systems comprised of purely heavy quarks. One could, for example, rely on the QCD sum rule method as in Ref. [20]. However, the small separation between quarks in the four-b system requires the use of higher-derivative moments to reliably estimate the non-perturbative bound state physics, which in turn calls into question the convergence of the perturbative expansion for such moments. Even for the two-b system, this approach can become less trustable than the similar charmonium calculations [21]. A separate approach is to treat the bb(bb) system as a composite diquark bound state and then calculate the interdiquark binding energy [22,23]. However, this simplified approach turns out to be inadequate as the average distance between the two diquarks is comparable to the separation between their constituent quarks, necessitating a complete four-body calculation.
In this paper, we adopt a phenomenological potential with its parameters determined by fitting to the twobody bb spectrum and verified by lattice simulation. We then numerically solve the Schrödinger equation to obtain the ground state energy and approximate wave function for the four-b tetraquark. After that, we discuss spindependent (SD) corrections and obtain a final estimate for the ground state mass.

II. CORNELL POTENTIAL, FLIP-FLOP, AND BUTTERFLY
For a two-body qq system, the Cornell potential V (r) = −4α s /(3 r) + r/a 2 has been widely used to understand bottomonia and charmonia spectroscopy [3,4]. The −1/r Coulomb term is understood simply as the spin-independent contribution from one-gluon exchange. The r/a 2 linear term is the long-range contribution due to QCD confinement. When one extends this phenomenological potential to many-body system, all fluxtube configurations must be checked to minimize the static potential. For three-body systems like baryons, lattice QCD calculations predict that only two flux-tube configurations matter; which one to choose depends on whether an interior angle of the triangle formed by the three quarks is greater than 120 • [24].
For the four-body system there are in general three relevant configurations, with the flux-tubes reconfiguring in such a way as to minimize the total potential. The first two configurations, which in combination is sometimes referred to as the "flip-flop" configuration, is shown in the left panel of Fig. 1. Depending on the relative distance between (r 13 , r 24 ) and (r 14 , r 23 ), one has the dimeson configuration of either The third, so-called "butterfly" configuration, is shown in the right panel of Fig. 1. This is the diquark-diquark configuration with two bottom quarks, b 1 and b 2 , forming a color-anti-triplet (or less-favored sextet), and the two anti-bottom quarks,b 3 andb 4 , forming a color-triplet. For the butterfly configuration, the middle two connecting points should be chosen to minimize the total fluxtube length. The effective potential for the di-meson system is either (14,23) , which differs only by exchangingb 3 ↔b 4 . The flip-flop potential is defined as the minimum of the two di-meson configurations for any four-particle phase space and is given by (13,24) , V di−meson (14,23) . ( The effective potential for the butterfly configuration is Here, L min is the minimum value of the total flux-tube length for all possible connecting points in the right panel of Fig. 1. The formula for the sextet diquark configuration is obtained by replacing the coefficients (−1/3, −2/3) by (−5/6, 1/3). As the sextet configuration has both attractive and repulsive Coulomb forces, it generically leads to a larger potential energy and thus will not contribute to the ground state energy calculation. The total four-quark potential is defined to be the minimum of the three possible configurations and is given by There is no a priori reason to expect that the two potential parameters, α s and 1/a, should have the same values for four-quark states as in two-quark states. Fortunately, the lattice QCD studies of Ref. [25] find consistency in this approach for a large set of four-quark spatial configurations. Therefore, we will use fitted values of α s and 1/a from the two-body quarkonia spectrum to calculate the solution to the four-body Schrödinger equation and obtain the ground state energy.

III. NUMERICAL CALCULATIONS BASED ON DIFFUSION MONTE CARLO
There are a plethora of ways to calculate the energy levels of many-body non-relativistic systems. For the four-lepton di-positronium molecule, variational methods with a very large amount of trial functions have been used to accurately obtain the ground state energy [9][10][11]. Minimizing the Hamiltonian for so many variational parameters, however, can be extremely computationally expensive (see Ref. [26] for a recent attempt). For our numerical calculation, we will instead adopt the more efficient but perhaps less accurate Diffusion Monte Carlo (DMC) method (see Ref. [27] for an introduction). To test the numerical calculation, we successfully reproduced the binding energy for the di-positronium molecule.
The central idea behind the DMC method is to replace real time by an imaginary time and adjust the guessed ground state energy based on the behavior of the wave function. The time-dependent wave function evolves as n e −iEnt Ψ n ( x), where E n and Ψ n ( x) are the true energy eigenvalues and eigenfunctions of the Hamiltonian, respectively. Scaling the energy eigenvalues E n by a constant guessed value E g makes no physical change to wave function. Thus, after making the substitutions E n → E n − E g and t → −iτ , the evolution of the wave function becomes n e −(En−Eg)τ Ψ n ( x). When E g ≅ E 0 , only the ground state wave function will be stable while the excited states will diffuse away. If E g > E 0 , the wave function will diverge, while for E g < E 0 even the ground state will diffuse away. By adjusting the value of E g based on the behavior of the wave function, we are able to obtain the correct ground state energy.
Practically, the wave function in the DMC algorithm is represented by random walks of many particles in phase space. To observe the behavior of the wave function with respect to E g , a "birth-death" mechanism is implemented such that when E g is too large, the particles will replicate themselves to increase the total number of particles, and vise versa. For each step, the algorithm adjusts E g based on the change in the total number of particles until it converges to the ground state energy E 0 . The stabilized walker distribution then gives us the ground state wave function.
For the system at hand, we also need to find L min for the butterfly configuration. This is similar to the wellknown Fermat and Torricelli problem to link three points with a minimal network, which is a special case of the Steiner tree problem in mathematics. The added connecting points in the middle are known as Steiner points, and the original points are called terminals. In twodimensional space, the Steiner tree problem has an analytical solution, however it is still an NP-complete problem in higher dimensions. To find the positions of the two Steiner points the Steiner configuration for the two vertices in the right panel of Fig. 1, we adopt Smith's algorithm [28]. This iterative algorithm determines an equation for each Steiner point, k, from summing all possible links surrounding it to both terminals and other Steiner points j is the position of point j after the i-th iteration. Usually after only 20 iterations this algorithm gives a solution that matches the true solution to a very high precision.

IV. BINDING ENERGY FOR SPIN-INDEPENDENT POTENTIAL
To calculate the tetraquark mass, we will adopt two sets of benchmark parameters for m b , Both of them can provide a good fit to the bottomonium spectra [3,4], although the second one requires a universal shift of energy levels of around −0.77 GeV, to take into account that the dynamic bottom-quark mass could be different from the bare heavy quark mass.
To check the stability of our numerical calculation, we fix the DMC parameters-initial particle number, time step, and simulation length-and modify the total potential by an overall factor of κ. This is equivalent to scaling the reduced mass and ground state energy while keeping numerical running conditions fixed. In dimensionless units, the binding energy is anticipated to scale as κ 2 , which is clear from Fig. 2. Furthermore, the ratio of additional binding energy, ∆E, over the di-meson binding energy, E 0 , for the flip-flop and flip-flop+butterfly potentials should be independent of κ, which is also approximately true up to small numerical fluctuations. From the lower panel of Fig. 2, it is clear that the flip-flop and flip-flop+butterfly configurations account for an additional 20% and 30% binding energy, respectively.  In Table I, we show the additional binding energy due to the different four-particle effective potential configurations. Within the error of our calculation, we have found that an additional ∼ 50 MeV and ∼ 80 MeV of binding energy can be attributed to the flip-flop or flip-flop+butterfly configurations. The total binding energy for the flip-flop+butterfly configuration is around 330 MeV for BM-I and 300 MeV for BM-II (without accounting for the constant energy shift of −0.77 × 2 GeV). These binding energies are large enough to suggest that the four b state should be treated as a true tetraquark system rather than a weakly-coupled molecular system.
Before we move on to discuss the spin-dependent cor- rections, we show the wave functions for the tetra-quark state. To generate the wave functions in Fig. 3 and simplify the multiple-dimension numerical integration, we have treated the wave-function as approximately flat in r 12 and r 34 and kept the remaining three variables: r 13 , r 14 and a relative angle between them. In the upper panel of Fig. 3, we show the squared wave function times r 2 13 for different static potentials, while in the lower panel we show the value of the squared wave function near the origin. One can see that the squared wave function for the flip-flop+butterfly potential is around 0.59 of the value for dissociated di-meson configuration at the origin. There is only around 5% difference between the flip-flop and flip-flop+butterfly squared wave functions. We also note that because of the symmetriesb 3 ↔b 4 and b 1 ↔ b 2 , the wave functions are identical when plotted in terms of r 14 , r 23 or r 24 .

V. SPIN-DEPENDENT CORRECTIONS
For a two particle system, short-range spin-dependent interactions contain a local δ-function force [29,30] with ∆C 2 as the difference of quadratic Casimir. For the four particle system, we anticipate the ground state spatial wave function to be predominantly S-wave so additional corrections proportional to L · S should be negligible. From Fig. 3, it is apparent that the wave function for the flip-flop configuration closely approximates the wave function for the flip-flop+butterfly configuration. Thus, to implement the SD correction for the 0 ++ tetraquark, we will focus on the flip-flop configuration.
The general ground state wave function for the flip-flop di-meson configuration is Here, R 1 represents the region with V di−meson (13,24) < V di−meson (14,23) and otherwise for R 2 ; (b 1b3 ) 1 represents the color-singlet contraction of b 1 andb 3 . The product of the spatial and color wave functions is symmetric under the interchange of b 1 ↔ b 2 andb 3 ↔b 4 , so the spin-zero wave functions should be The relative minus between the above two terms is necessary to satisfy the Pauli exclusion principle, and provides the lowest ground state energy after hyperfine splitting.
Calculating the matrix element Ψ|H SD |Ψ , the spindependent correction for the 0 ++ ground state is approximately for BM-I with α s (2 m b ) ≈ 0.2. The result for BM-II is similar. The symmetries b 1 ↔ b 2 andb 3 ↔b 4 imply that the contributions of R 1 and R 2 to the matrix element are the same, and thus we double the expression for region R 1 in our calculation. Here, the relative error is taken to be O(α s ) ≈ 20% [31] and should only be used as guidance from the theoretical calculation. Altogether, the energy for the ground state 0 ++ mode is which is below the energy threshold of 2M (η b ) = 18.798 GeV and 2M [Υ(1S)] = 18.920 GeV.
Within the framework of our calculation, one could also calculate heavier states including spin-one and spintwo excitations. The detailed mass spectrum requires one to calculate the spatial excitation energy as well as the full wave-functions in terms of all degrees of freedom. These are conceptually straightforward, but numerically complicated. Furthermore, one could also apply our calculation procedure to other four-heavy quark system like cccc and ccbb. For instance, the ground state, 0 ++ , of cccc is estimated to be below twice of J/Ψ mass and can have the similar decay of 0 ++ → Ψ(1S)Ψ(1S) * → ℓ + ℓ − ℓ + ℓ − , which has a smaller branching ratio and different S/B for the experimental searches.
In summary, based on the static potential of the fluxtube model for four heavy quark interactions, we have used a diffusion Monte Carlo algorithm to numerically solve the many-body non-relativistic Schrödinger equation. We have found a ground state, 0 ++ , with a mass of 18.69 ± 0.03 GeV, which is approximately 100 MeV below twice the mass of η b expected for a disassociated di-meson ground state. Here, the error of 30 MeV is chosen to include potentially next-to-leading order SD corrections and should not be taken too seriously.
We thank Joshua Berger, Geoffrey Bodwin, Estia Eichten, and Zhen Liu for useful discussion. This work is supported by the U. S. Department of Energy under the contract DE-FG-02-95ER40896.