Diffusion Monte Carlo calculations of fully-heavy (multiquark) bound states

We use a diffusion Monte Carlo method to solve the many-body Schr\"odinger equation describing fully-heavy tetraquark systems. This approach allows to reduce the uncertainty of the numerical calculation at the percent level, accounts for multi-particle correlations in the physical observables, and avoids the usual quark-clustering assumed in other theoretical techniques applied to the same problem. The interaction between particles was modeled by the most general and accepted potential, i.e. a pairwise interaction including Coulomb, linear-confining and hyperfine spin-spin terms. This means that, in principle, our analysis should provide some rigorous statements about the mass location of the all-heavy tetraquark ground states, which is particularly timely due to the very recent observation made by the LHCb collaboration of some enhancements in the invariant mass spectra of $J/\psi$-pairs. Our main results are: (i) the $cc\bar c\bar c$, $cc\bar b\bar b$ ($bb\bar c\bar c$) and $bb\bar b \bar b$ lowest-lying states are located well above their corresponding meson-meson thresholds; (ii) the $J^{PC}=0^{++}$ $cc\bar c\bar c$ ground state with preferred quark-antiquark pair configurations is compatible with the enhancement(s) observed by the LHCb collaboration; (iii) our results for the $cc\bar c\bar b$ and $bb\bar c\bar b$ sectors seem to indicate that the $0^+$ and $1^+$ ground states are almost degenerate with the $2^+$ located around $100\,\text{MeV}$ above them; (iv) smaller mass splittings for the $cb\bar c\bar b$ system are predicted, with absolute mass values in reasonable agreement with other theoretical works; (v) the $1^{++}$ $cb\bar c\bar b$ tetraquark ground state lies at its lowest $S$-wave meson-meson threshold and it is compatible with a molecular configuration.


I. INTRODUCTION
The J/ψ signal was observed simultaneously at Brookhaven [1] and SLAC [2] in 1974; it was a heavy resonance with a surprisingly small decay width. Three years later, an even heavier resonance but equally narrow, the so-called Υ state, was observed at Fermilab [3,4]. The interpretation of the J/ψ and Υ as low-lying bound states of a heavy quark, Q, and its antiquark,Q, with Q either c-or b-quark, explained their narrow decay widths and, in fact, was proved to be crucial to establish Quantum Chromodynamics (QCD) as the strong-interaction sector of the Standard Model of Particle Physics [5,6]. It is also worth mentioning that a related system, the cb bound state (B + c ), has also been found in nature [7]; and that the heaviest of the quarks, i.e. the top-quark, was discovered in 1995 at Fermilab [8], with a mass around 175 GeV and a large decay width that, due to weak interactions, forbids to form narrow tt resonances.
Heavy quarkonia, viz. mesons containing only a heavy valence quark-antiquark pair, opened the possibility to use a non-relativistic (NR) picture of QCD. They can indeed be classified in terms of the quantum numbers of a NR bound state, and the spacings between radial, orbital and spin excitations have a pattern similar to the ones observed in positronium, an e + e − NR bound state well studied in Quantum Electrodynamics (QED) [9,10]. Being baryonic analogues of heavy quarkonia, triply-heavy baryons may provide a complementary window in the un-It is a fact that many precise experimental results are available for conventional heavy quarkonia [13]. In addition, tens of charmonium-and bottomonium-like states, the so-called XYZ states, which cannot fit the quark model picture, have been identified along the last two decades at B-factories (BaBar, Belle, and CLEO), τcharm facilities (CLEO-c and BES) and hadron-hadron colliders (CDF, D0, LHCb, ATLAS, and CMS). So far, there is no consensus about the nature of these exotic states (see Refs. [14][15][16][17] for reviews of the experimental and theoretical status of the subject). Their analysis and new determinations will continue with the upgrade of experiments such as BES III [18], Belle II [19], and HL-and HE-LHC [20]. This will provide a sustained progress in the field as well as the breadth and depth necessary for a vibrant heavy quark research environment.
The ultimate aim of theory is to describe the properties of the XYZ states from QCD's first principles. However, since this task is quite challenging, a more modest goal is to start with the development of QCD moti-vated phenomenological models that specify the colored constituents, how they are clustered and the forces between them. In that line, simultaneously to the experimental measurements, theorists have been proposing for the XYZ states different kinds of color-singlet clusters, made by quarks and gluons, which go beyond conventional mesons (quark-antiquark), baryons (three-quarks) and antibaryons (three-antiquarks); the most famous are glueballs, quark-gluon hybrids and multiquark systems (for a graphic picture of these kinds of hadrons see, for example, Figs. 1, 6, and 7 of Ref. [16]). Related to the last ones, the first and best known exotic state is the X(3872), which was observed in 2003 as an extremely narrow peak in the B + → K + (π + π − J/ψ) channel and at exactly thē D 0 D * 0 threshold [21,22]. It is suspected to be a cncn (n = u or d quark) tetraquark state whose features resemble those of a molecule, but some experimental findings seem to point out the existence in its wavefunction of more compact components such as diquark-antidiquark and quark-antiquark [23][24][25][26][27][28].
Finally, fully-heavy tetraquarks have recently received considerable attention, both experimentally and theoretically. On the experimental side, it is thought that allheavy tetraquark states will be very easy to spot because their masses should be far away from the typical mass regions populated by both conventional heavy mesons and the XYZ states discovered until now. A search for deeply bound bbbb tetraquark states at the LHC was motivated by Eichten et al. in Ref. [29], and it was carried out by the LHCb collaboration [30] determining that no significant excess is found in the µ + µ − Υ(1S) invariant-mass distribution. On the other hand, the LHCb collaboration has recently released in Ref. [31] a study of the J/ψ-pair invariant mass spectrum finding a narrow peak and a broad structure which could originate from hadron states consisting of four charm quarks.
From the theoretical side, concerning the interaction between heavy quarks, chiral symmetry is explicitly broken and thus meson-exchange forces cannot exist in a fully heavy tetraquark system [32,33], which would favor the formation of genuine tetraquark configurations rather than loosely bound hadronic molecules. This is interesting by itself, but also it simplifies the kind of quark-(anti-)quark interactions to be taken into account, justifying the proliferation of theoretical works.
The goal of the present study is to achieve the most general and accurate prediction for the ground states of fully-heavy tetraquarks. In order to comply with the first feature, we are not going to assume any particular clustering between the valence quarks (antiquarks), the interaction between them is the most simple and accepted one: Coulomb + linear-confining + hyperfine spinspin, and it will be implemented non-perturbativelly. 1 The second feature, accurateness, is fulfilled using a diffusion Monte Carlo (DMC) technique for solving the manybody Schrödinger equation wich, in contrast with variational methods, allows to reduce the uncertainty of the numerical calculation at the percent level, since the systematic one associated with the trial wave function is eliminated by the algorithm. Note, too, that it is possible to disentangle the uncertainty inherent of manybody techniques from the theoretical one coming from the model.
The manuscript is arranged as follows. In Sec. II the theoretical framework is presented; we explain first the origin, features and implementation of the computational algorithm and, later, the quark model Hamiltonian and how their parameters are fixed. Section III is mostly devoted to the analysis and discussion of our theoretical results on fully-heavy tetraquarks; note here that we firstly study all-heavy mesons and baryons, comparing our results with those available from variational methods in order to confirm the validation of our approach. Finally, we summarize and give some prospects in Sec. IV.

A. Quark model
The Hamiltonian which describes fully-heavy boundstate systems can be written as where m i is the quark mass, p i is the momentum of the quark, and T CM is the center-of-mass kinetic energy. Since chiral symmetry is explicitly broken in the heavy 1 It is fair to notice that a similar calculation to the one presented herein has been recently released in Ref. [58]; however, important differences must be mentioned: (i) meson-meson and diquarkantidiquark clusters were assumed, (ii) the sextet-antisextet diquark-antidiquark configuration was fully neglected, and (iii) the hyperfine spin-spin interaction was computed perturbatively.
quark sector, the two-body potential, V ( r ij ), can be deduced from the one-gluon exchange and confining interactions; i.e.
It is important to highlight that the potentials above have tensor and spin-orbit contributions, which shall be neglected in this work. This is because our main purpose here is to get a first, unified and non-perturbative reliable description of fully-heavy hadron ground states (from 2to 4-quark systems) within the DMC method. Moreover, this kind of interactions appeared not to be essential for a global description of baryons [59], and beyond [54].
The Coulomb and hyperfine terms are collected in the one-gluon exchange potential, and are given by where α s is the strong coupling constant fixed in a phenomenological way, λ are the SU (3)-color Gell-Mann matrices, and the Pauli spin matrices are denoted by σ. The Dirac delta function of the hyperfine term comes from the Fermi-Breit approximation of the one-gluon exchange interaction. In order to perform non-perturbative calculations, the δ (3) ( r ij ) is usually replaced by a smeared function that, in our case, reads as follows with κ a quark model parameter, and r 0 = A 2mimj mi+mj B a regulator which depends on the reduced mass of the quark-(anti-)quark pair. Confinement is one of the crucial aspects of the strong interaction that is widely accepted and incorporated into any QCD based model. Studies of QCD on a lattice have demonstrated that multi-gluon exchanges produce an attractive linearly rising potential, which is proportional to the distance between infinitely heavy quarks [60]. This phenomenological observation is usually modeled as where b is the confinement strength and ∆ is a global constant fixing the origin of energies. Table I shows the quark model parameters relevant for this work. Note here that we are using the so-called AL1 potential proposed by Silvestre-Brac and Semay in Ref. [61], and applied extensively to the baryon sector in Ref. [59]. It is worth emphasizing that the potential collects nicely the most important phenomenological features of QCD for heavy quarks, and that the parameters were constrained by a simultaneous fit of 36 mesons and 53 baryons with a remarkable agreement with data.

B. Computational algorithm
Quantum Monte Carlo (QMC) methods have been successfully applied to many research areas but quantum chemistry and material science are the ones which have received more attention [62][63][64]. This is because QMC is a natural competitor of other methods where the uncorrelated or Hartree-Fock state does not provide a good description of the many-body ground state. Other applications of QMC algorithms are solid-state physics concerning the dynamics of condensed helium systems [65,66], and studies on the properties of both bosonic and fermionic ultracold quantum gases [67][68][69][70].
Since nuclear Hamiltonians induce strong correlations, QMC methods have appeared to be very valuable in the understanding of nuclei and nucleonic matter. Variational Monte Carlo (VMC) algorithms dealing with nuclear interactions were introduced in the early 1980s [71]. Afterwards, methods based on Green Function Monte Carlo (GFMC) burst into nuclear physics in the late 1980s [72,73], and were applied mostly to spin-isospindependent Hamiltonians. The GFMC technique is very accurate but becomes increasingly more complex when moving toward larger systems, being 12 C the state-ofthe-art studied one [74][75][76]. Diffusion Monte Carlo methods [77] appear to be much more efficient at treating large systems; however, there are unsolved issues when dealing with spin-isospin dependent potentials.
The application of QMC methods to hadron physics has been scarce, basically because most known hadrons consist on 2-and 3-body relativistic bound states. However, the description from different perspectives of composite states with four and more, light or heavy, quarks (antiquarks) has recently received considerable attention since the discovery of the potentially non-relativistic, 4quark, charmonium-like system X(3872). It is in this context that QMC algorithms can contribute to shed some light to the study of tetraquark, pentaquark, hexaquark, etc., non-relativistic systems.
Up to our knowledge, the first QMC study of mesons and baryons was performed by Carlson et al. in Refs. [78,79]. The authors used a VMC method developed previ-ously for nuclear physics problems and their results compared reasonably well with those of the well-known Isgur-Karl's quark model [80][81][82][83]. The Diffusion Monte Carlo method has been applied recently to the fully-beautiful tetraquark system in Ref. [58]; in particular, the authors calculate the ground state energy of the J P C = 0 ++ bbbb system.
The central idea behind the DMC method is to write the Schrödinger equation for n-particles in imaginary time ( = c = 1): where E s is the usual energy shift used in DMC methods, R ≡ ( r 1 , . . . , r n ) stands for the position of n particles and α denotes each possible spin-color channel, with given quantum numbers, for the n-particles system. The function Ψ α (R, t) can be expanded in terms of a complete set of the Hamiltonian's eigenfunctions as where the E i are the eigenvalues of the system's Hamiltonian operator,Ĥ. The ground state wave function, φ 0,α (R), is obtained as the asymptotic solution of Eq. (6) when t → ∞, as long as there is overlap between Ψ α (R, t = 0) and φ 0,α (R), for any α-channel.
A crucial feature of QMC methods is the use of importance sampling techniques [84], in order to reduce the statistical fluctuations to a manageable level. Given a Hamiltonian, H = H 0 + V , of the form with m the mass of each particle which composes the bound-state system, if one works with the function where ψ(R) is the time-independent trial function: where i, j are indices which run over the number of particles. Within this work, we choose φ( r ij ) as the wave function solution of the 2-body Hamiltonian of the system at short interquark distances: φ( r ij ) = e −aij rij , where a ij are determined by the so-called cusp conditions. Equation (6) turns out to be are the so-called local energy and drift force, respectively. The formal solution of Eq. (11) is given by While algorithms based on GFMC implement the whole Green's function, DMC methods rely on reasonable approximations of G α α (R , R, ∆t) for small values of the time step ∆t and iterates repeatedly to obtain the asymptotic solution f α (R, t → ∞). In our case, the Green's function is approximated by (15) where for the first part we follow Ref. [85] and approximate it as which is exact up to order (∆t) 2 . With the expression above, Eq. (14) becomes with Note herein that in the case of having more than one α-channel, we follow the method proposed in Ref. [86] and propagate the quantity such as where we have introduced the weight factor With either one chanel or more, the effective way of applying the DMC method defined by Eq. (22) is the following. Each walker is characterized by the positions, R, and the coefficients of each channel, c α ; all of them fixed by an initial guess. Then, for each walker, (i) Move it, under the drift force F (R), during an interval ∆t/2 with accuracy (∆t) 2 .
(iv) Randomly replicate the walker considering the product of branching ratios and where the last one is derived from Eq. (23), and therefore the coefficients are updated according to The procedure above must be repeated for each walker until the set of them is exhausted. The resulting set of walkers corresponds to the new positions and coefficients {R , c α }. Finally, the whole procedure must be repeated as many times as needed to reach the asymptotic limit t → ∞.

III. RESULTS
A detailed discussion about the particular features of our spectrum will be given in the following subsections. However, a comment is due here on the theoretical uncertainty of our results. There are two kind of theoretical errors: one is inherently connected to the statistical nature of the DMC algorithm and the other one is related with a shortcoming of the quark model approach and lies on the way to fix the model parameters. The statistical error is of the order of 1 MeV and negligible with respect the systematic one related with the quark model. As mentioned above, the set of model parameters are fitted to reproduce a certain number of hadron observables within a determinate range of agreement with experiment. Therefore, it is difficult to assign an error to those parameters and, as a consequence, to the magnitudes calculated when using them. As the range of agreement between theory and experiment is around 10 − 20%, this value can be taken as an estimation of the model uncertainty for fully-heavy tetraquark systems.

A. Fully-heavy mesons
Four fundamental degrees of freedom at the quark level: space, spin, flavor and color are generally accepted in QCD, and any hadron's wave function must be expressed as a product of these four terms: From now on, we shall drop out the trivial flavor wave function because only fully-heavy hadrons are considered in this work and thus, due to the explicit breaking of chiral symmetry, the considered quark model Hamiltonian is blind to the heavy quark flavor.
Concerning the color degree of freedom, any hadron state must be a color singlet one, because no color charges has been observed in nature. The leading Fock state of a fully-heavy meson is constituted by a quark, Q, and antiquark,Q, with Q either c or b-quark. Therefore, the SU (3) color wave function is constructed as follows: where the singlet state is the physically interesting one in this work, and it is given by the well know, symmetric expression, Since a meson is made by distinguishable particles (a quark and an antiquark), there is no restriction in the spin wave functions due to Pauli principle and thus all possibilities can be considered. The S = 0 state is antisymmetric respect the particle exchange 1 ↔ 2, whereas the S = 1 wave functions are all symmetric. One can guess that an excitation of a unit of angular momentum costs an energy around 500 MeV. This effect can be estimated from the experimental M (L = 1) − M (L = 0) mass differences: MeV, but also in the baryon sector as, for instance, N (1535) − N (940) ≈ 570 MeV. Therefore, we shall limit ourselves to analyze states with orbital angular momentum equal to zero. An important consequence is that the space-wave function will always represent an S-wave state, which is completely symmetric.
Once the (trial) wave function of the meson is constructed, we follow the DMC algorithm explained in the section above and obtain, as a proof of concept, the masses of the singlet and triplet S-wave ground states of heavy quarkonia. Our results are shown in Table II, they compare fairly well with the experimental data despite the simplicity of the quark model; this also supports our confidence on the mass prediction for fullyheavy tetraquark bound states. Our numerical technique is contrasted with the same calculation but using a variational approach [54]. As one can see in Table II, there are negligible differences between the two methods; as expected, our values are always below the ones reported in TABLE II. The mass spectra of the heavy quarkonia in units of MeV. The 1st, 2nd and 3rd columns refer to the name and quantum numbers of the considered hadron, 4th column is our result, 5th column is the same calculation but using a variational method [54], and 6th column collects experimental data if exist. n 2S+1 LJ J P C DMC VAR [54] EXP [13] ηc Ref. [54]. This validates our technique against the variational one, which should give an energy's upper limit quite close to the real eigenvalue in the case of meson systems. As we shall see later, the differences between the two numerical approaches will be larger as the number of particles increases. The concept of radial distribution function can be applied to multiquark systems and, in fact, can provide valuable information about the existence of interquark correlations; in particular, 2-body correlations. If the nparticle wave function is defined as ψ( r 1 , . . . , r n ), where spin, flavor and color degrees of freedom have been ignored for simplicity without lost of generalization, the probability of finding particle 1 in position r 1 , particle 2 in position r 2 , . . ., particle n in position r n is: and it is normalize to one, i.e.
Therefore, one can define which expresses the probability of finding 2 particles in positions r 1 and r 2 ; and the radial distribution function as where r indicates now the distance between the two particles considered. Figure 1 shows, for the studied mesons, the radial distribution functions which are pure estimators calculated within our DMC following Ref. [87,88]. Among the features one can observe, the following are of particular interest: (i) the cc states are the most extended objects 0.041 Υ(1S) 0.043 and the system becomes more compact in going from the cc meson to the bb one, being the size of the bb system almost half of the cc-meson's one; (ii) the S = 1 state is slightly more extended than the S = 0 state because the different sign of the spin-spin hyperfine interaction; and (iii) the structural differences between S = 0 and S = 1 states seems to blur as we go to heavier quarks, as expected since the hyperfine mass splitting gets smaller.
Finally, the typical mean-square radii, r 2 , of the studied quarkonium systems are shown in Table III. They compare nicely with the results reported by a well-known non-relativistic QCD effective field theory for quarkonium [89,90]. Moreover, looking at the table, we can confirm that the interparticle distance of the bb system is almost half of the cc one; and the B c 's one is closer to the bottomonium than to the charmonium.

B. Fully-heavy baryons
We turn now our attention to the study of all-heavy ground state baryons and thus no orbital angular momentum excitations must be considered, i.e. only S-wave symmetric states are under scrutiny.
Concerning the SU (3) color wave function, it is con-structed as follows: where the colorless state is fully anti-symmetric and its color wave function is given by the textbook expression With three particles of spin 1/2, one can construct the following spin states: examples of the spin wave functions used herein, without lost of generality. For the baryons the spin wave function must comply Pauli statistics in the case that quarks are the same. Namely, the S sym. = 3/2 S IV. Masses, in MeV, of the ground states of fullyheavy baryons. The 1st, 2nd and 3rd columns refer to the name and quantum numbers of the considered hadron, 4th column is our result, and 5th column is the same calculation but using a variational method [59]. Note that, for comparison, the masses presented here include three-body-force corrections using the value of the constant C reported in Table  1  wave function, which is completely symmetric, must be used for the Ω ccc and Ω bbb baryons; whereas the S sym. = 1/2 MS wave function, which is symmetric respect the two particles that are equal, must be used for the ground states of Ω ccb and Ω cbb , because it is an allowed state with a lower energy than the S = 3/2 case. Table IV shows the calculated masses for the ground states of the Ω QQQ baryons (Q = c or b) in each allowed spin channel, and compares them with those obtained in Ref. [59]. As one can see, there are negligible differences between the two numerical approaches. Unfortunately, there is no experimental data to compare with; therefore, we encourage the design of experimental set-ups at, for instance, the LHC@CERN facility able to detect this kind of particles because the reward could be high and, as mentioned above, triply-heavy baryons are ideally suited to study QCD and, in particular, the heavy quark-(anti-)quark interaction, as it has been the case for heavy quarkonia.
In order to check further the capabilities of the DMC method to describe heavy baryons, we calculate several structural (static) properties to compare them with the results obtained in Ref. [59]. In particular, the following ones: They are presented in Table V and one can see that the values obtained in Ref. [59] are perfectly reproduced. Therefore, as a proof of concept, our numerical technique can be applied to study not only eigenenergies but also structure properties of the hadrons under consideration when the momentum of the prove is less or equal than the typical hadron scale, 1 GeV. Equation (45) gives an idea of the physical size of the baryon, where is the center-of-mass coordinate, with M = m 1 +m 2 +m 3 the total mass of the system. One can see in Table V that the results lie between 0.02 and 0.07 fm 2 . This means that the spatial extension of such baryons goes from 0.15 to 0.26 fm while for the mesons we had triply-heavy baryons are objects only slightly less compact than their heavy quarkonia counterparts. This can be explained by the fact that the color quark-quark interaction is half as intense as the quark-antiquark one but, a priori, one would expect a bigger effect and this could indicate that some kind of diquark correlations be at play.
Equation (46) defines the charge mean-square radii that can be deduced from the value of the baryon's electric form factor at the photon point extracted from lepton-baryon scattering experiments. A result that might be worth to emphasize is the electric charge radius of the positive charged baryon Ω + ccb , which is 0.31 fm, three times smaller than the proton's radius 0.84 fm [13].
If mesonic exchange currents and relativistic effects can be ignored, the magnetic moment operator is given by Eq. (47), which is just the sum of the magnetic moments of each quark, with orbital and spin contributions. Our values for the positive-charged and neutral triply heavy baryons, Ω + ccb and Ω 0 cbb , are respectively 0.475µ N and −0.193µ N . These can be compared with the values of the proton and neutron: 2.79µ N and −1.91µ N , collected in Ref. [13]. Namely, the triply-heavy baryon partner of the nucleon have a magnetic moment 5−10 times smaller. Figure 2 shows the relevant radial distribution functions. Among the observed features, the following are of particular interest. In the cases of Ω ccc and Ω bbb baryons, the quark-quark probability distribution function is broader than its quark-antiquark counterpart. There are two kinds of probability distributions in the and Ω bbb (bottom-right panel ) baryons. The (pink) long-dashed curve plotted in the Ωccc and Ω bbb panels is the same radial distribution function but for the cc and bb mesons, respectively.
cases of Ω ccb and Ω cbb baryons, one referring to the QQpair where the two heavy quarks are equal and the other one when the QQ-pair is made with different species. The same pattern as in heavy quarkonia is observed in triplyheavy baryons where the QQ-pair seems to be more compact as the heavy quark mass is larger, and thus this could facilitate the appearance of strong (cb)-and (bb)diquark correlations inside the Ω ccb and Ω cbb baryons, respectively. Finally, the typical mean-square radii, r 2 ij , of the studied triply-heavy baryons are collected in Table VI. Among other features, it highlights that the cb-and bbpair are closer inside the Ω ccb and Ω cbb baryons, respectively.
For completeness, we have calculated the J P = 3/2 + lowest-lying states of the Ω ccb and Ω cbb baryons. Their masses are 8046 MeV and 11247 MeV, respectively. They compare well with the variational ones reported in Ref. [59]. We expect bigger differences when radial, angular and spin excitations will be compared; this goes beyond the scope of the present manuscript but indicates a possible next step to follow in the future. Finally, the predicted mass splittings ∆m = m(3/2 + ) − m(1/2 + ) are 28 MeV and 32 MeV for the Ω ccb and Ω cbb baryons; they are of the same order of magnitude than those collected in Ref. [91], and references therein.

C. Fully-heavy tetraquarks
Diffusion Monte Carlo methods are designed to compute non-relativistic bound-states of a few-to manyparticles system. This technique has been successfully applied to the field of nuclear physics studying light and middle nuclei, as well as objects of high nuclear density such as neutron stars. However, as explained before, applications of DMC to hadron physics are scarce because most hadrons were understood as relativistic bound states of few (two or three) light quarks (antiquarks). Nowadays, many experimental signals point out the existence of tetra- [31], penta- [92,93] and even hexaquark [94] systems, mostly in heavy quark sectors where non-relativistic dynamics could be thought as a good assumption, and thus Monte Carlo techniques are becoming attractive approaches to apply in hadron physics. After the application of the method to mesons and baryons as a proof-of-concept, our next step is the study of allheavy tetraquarks. Further studies of more complex multiquark systems shall be performed in the future but they go beyond the scope of this manuscript.
Multiquark systems present richer color structures than mesons and baryons. Without assuming any kind of clustering, the color algebra applied to tetraquark states leads to the following irreducible representations where two color singlet states appear. They are usually known as, respectively, the (3 c ⊗3 c ) and (6 c ⊗6 c ) diquarkantidiquark configurations, i.e.
Knowing that3 c diquark and 3 c antidiquark representations are antisymmetric under the transposition of the two particles: where αβγ is the Levi-Civita tensor, with the greek letters going from 1 to 3; and that 6 c diquark and6 c an-tidiquark are symmetric under the same transposition: 2 one can build two orthogonal singlet color tetraquark states Their explicit expressions in terms of the familiar red, green and blue degrees-of-freedom, and without explicitly clustering, can be written down as which is antisymmetric under the exchange of either both quarks or both antiquarks, and which is symmetric under the exchange of either both quarks or both antiquarks. It is also important to notice herein that the two color states defined above can be expressed in another set of color representations called meson-meson, 1 c ⊗ 1 c , and color-hidden, 8 c ⊗ 8 c , states.
If we turn now our attention to the spin degree-offreedom, the QQQQ (Q = c or b) system, made by fermions of spin 1/2, can have total spin S = 0, 1 and 2. There are two linearly independent S = 0 wave functions that can be written as They are, respectively, symmetric (S) and antisymmetric (A) under the exchange of both quarks and both antiquarks. The linearly independent S = 1 wave functions are given by (S z = S is assumed without lost of generality): which correspond to the symmetric − antisymmetric, antisymmetric − symmetric and symmetric − symmetric exchange of quarks and antiquarks. Finally, the S = 2 spin wave function can be written as where S z = S is again assumed without lost of generality. Note that this spin wave function is fully symmetric under the exchange of both quarks and both antiquarks. We shall focus our attention to the ground states of QQQQ, with Q either c-or b-quark in any possible combination. This implies that the space wave function is totally symmetric and, in order to fulfill Fermi-Dirac statistics, the wave-function configurations of Table VII must be taken into account for each tetraquark sector and J P (C) quantum numbers. We perform a coupledchannels calculation, based on the DMC algorithm explained above, for those sectors shown in Table VII which have more than one spin-color configuration. It is worth highlighting that the treatment of spin-dependent potentials and coupled-channels calculations are challenging features for Monte Carlo methods which have avoided their extended application to hadron physics, both are considered herein. Now, let us proceed to describe in detail our theoretical findings for each sector of fully-heavy tetraquarks.
The cccc tetraquark ground states. The masses of the S = 0, 1 and 2 ground states of cccc tetraquarks are shown in Table VIII. We compare first our results with those obtained when solving the same Hamiltonian but using a Rayleigh-Ritz variational method. One can see that both approaches are compatible; however, the differences are larger than the ones shown in the meson and baryon sectors because the variational methods begin to have difficulties as the number of particles grow. Nevertheless, the disagreement between both approaches is between 8 and 20 MeV, which is smaller than the usual uncertainty assigned to any quark model. Note, too, that the variational method only provides an upper limit of the eigenenergy and it depends on the goodness of the trial wave function; DMC results depend less on such details and, in principle, provides an exact estimate of the eigenenergy.
Table VIII provides also a comparison with numerous works that reported results on cccc tetraquarks using a large variety of techniques. It seems that our results 6.35 GeV, 6.44 GeV and 6.47 GeV for the ground states with quantum numbers J P C = 0 ++ , 1 +− and 2 ++ , respectively, are located just in the middle of the ranges covered by all approaches which are [5.97 − 6.80] GeV, [6.05 − 6.90] GeV and [6.09 − 6.96] GeV. These mass ranges are in agreement with the recently observed structures in the invariant mass distribution of J/ψ-pairs [31]. The cccc ground state with quantum numbers J P C = 0 ++ is about 400 − 500 MeV above the η c η c and J/ψJ/ψ thresholds. This suggests that the J P C = 0 ++ state is unstable and can decay into η c η c and J/ψJ/ψ final states through quark rearrangements. The J P C = 1 +− state lies about 400 MeV above the mass threshold of η c J/ψ, while J P C = 2 ++ is about 300 MeV above the mass threshold of J/ψJ/ψ, they can also easily decay into such charmonium pairs through quark rearrangement. Figure 3 shows for the cccc ground states, the relevant radial distribution functions. It is very interesting to observe that quark-antiquark pairs are slightly closer than quark-quark ones in the S = 0 (J P C = 0 ++ ) state whereas the situation gradually changes as spin increases. In any case, the cc and cc radial distribution functions are very similar with a mean value close to ∼ 0.5 fm, indicating that our tetraquark structures are compact objects and not meson-meson molecular states with typical inter-meson distances of 1 fm.
In order to quantify our observations, we have computed the mass mean-square radii of the J P C = 0 ++ , 1 +− and 2 ++ cccc ground states. They are, respectively, 0.087 fm 2 , 0.092 fm 2 and 0.095 fm 2 ; i.e. they are comparable with those obtained for mesons and baryons. Moreover, Table IX shows the typical mean-square radii for each (anti-)quark-quark pair indicating that no interquark distance is very different from the rest and thus all pairs play a similar role with distances of the order of ∼ 0.5 fm.
The bbbb tetraquark ground states. Let us now turn our attention to the J P C = 0 ++ , 1 +− and 2 ++ TABLE VII. Spin-color configurations for fully-heavy tetraquark systems.

System
J P (C) Spin-Color configurations cccc, bbbb ccbb (bbcc)   Table X. We compare first our results with those obtained when solving the same Hamiltonian but using a Rayleigh-Ritz variational method. One can see that even larger differences, between 30 and 50 MeV, are found in the bbbb sector, with our results always below the ones reported in Ref. [54], as expected. From the dynamical point of view, we would expect better agreement in tetraquark sectors where heavier quark masses are involved and thus the issue might be related with the oscillating parameters of the Gaussian basis used in Ref. [54] for the bbbb sector. Another possibility is that, as it will be shown later, interquark distances are smaller in these tetraquarks and thus multi-body correlations could play a more important role. Table X provides also a comparison with numerous works that reported results on bbbb tetraquarks using a large variety of theoretical techniques. Again, our results 19.20 GeV, 19.28 GeV and 19.29 GeV for the ground states with quantum numbers J P C = 0 ++ , 1 +− and 2 ++ , respectively, are located just in the middle of the ranges covered by all approaches which are [18.46 − 20.16] GeV, [18.32 − 20.21] GeV and [18.32 − 20.24] GeV. Figure 4 shows the relevant radial distribution functions for the bbbb ground states. As in the case of fullycharm tetraquarks, the quark-antiquark pairs are closer than quark-quark ones in the S = 0 (J P C = 0 ++ ) state. The situation, however, gradually changes as spin increases and thus the diquark (antidiquark) distance seems to be slightly smaller for J P C = 1 +− and 2 ++ bbbb tetraquarks. The theoretical fact that the J P C = 0 ++ bbbb ground state could prefer to be in   quark-antiquark pairs, together with a predicted mass which is 300 − 400 MeV above the η b η b and Υ(1S)Υ(1S) thresholds, could provide founded reasons to keep looking for (reasonably wide) bumps in the invariant mass of bottomonium-pairs at the LHCb experiment.
We report in Table XI the typical mean-square radii of the studied bbbb tetraquarks. Such table highlights again the change of configuration picture from quarkantiquark pairs to diquark-antidiquark ones, when going from S = 0 to S = 2 ground states of bbbb tetraquarks. Finally, the mass mean-square radii for the J P C = 0 ++ , 1 +− and 2 ++ ground states are, respectively, 0.026 fm 2 , 0.028 fm 2 and 0.029 fm 2 . This indicates that these states are very compact and far away from the picture of mesonmeson molecules.
The ccbb (bbcc) tetraquark ground states. A similar theoretical calculation must be performed for computing the ground states of the J P = 0 + , 1 + and 2 + ccbb (bbcc) tetraquarks. We obtain the masses: with the variational calculation [54]. As one can see, the differences between the two numerical methods are around 20 MeV, and our values are always below those obtained by the variational method. Note, too, that our results are in reasonable agreement with the scarce ones reported by other theoretical methods [33]. Moreover, our predictions are above their lowest open-flavor decay channels for about 240 − 280 MeV and thus they should appear as resonances with relatively large widths in the invariant mass of B c B c , B c B * c , or B * c B * c final states. An interesting feature distinguishes this sector with respect the two already discussed, i.e. the cccc and bbbb sectors. The expected pattern of interquark distances for cc, cb andbb (from larger to shorter) is conserved when changing the total spin S = 0, 1 and 2. An example of the output is drawn in Fig. 5 for the J P = 0 + ccbb ground state. We report, too, the mass mean-square radii for the J P = 0 + , 1 + and 2 + ground states, which Radial distribution functions for the J P = 0 + ground state of ccbb tetraquark. The case bbcc is obviously equal. No further information is obtained for the J P = 1 + and 2 + cases. are 0.046 fm 2 , 0.047 fm 2 , 0.048 fm 2 , respectively. These values lie between the ones reported above for the cccc and bbbb tetraquarks The cccb and bbcb tetraquark ground states. The cccb, bbcb (and cbcb) systems have not been studied before with neither the Hamiltonian used herein nor the variational method. Therefore, we cannot analyze in these cases the goodness of our numerical framework, diffusion Monte Carlo, with respect the variational one. We, however, shall compare our results with those obtained by other theoretical approaches, and numerical techniques, but such works are scarce because the complexity of the coupled-channels calculation needed to study the cccb, bbcb and cbcb tetraquarks.
The cccb and bbcb systems share some common features in terms of heavy quark symmetry and thus they will be discussed together. Table XII shows our spectrum of ground states with quantum numbers J P = 0 + , 1 + and 2 + . As one can see, in both tetraquark sectors, the 0 + and 1 + states are almost degenerate, with the 2 + ground state lying around 100 MeV above them. The same picture is drawn by Ref. [36] while smaller mass splittings are predicted in Ref. [33]. It is interesting to observe, too, that our absolute figures are in disagreement with the other two cases reported in Table XII, with results of Ref. [36] much higher than ours and those of Ref. [33]; the latter being of the same order of magnitude than ours. Another feature needed to be mentioned is that the lowest strong-decay threshold is η c B c (η b B c ) for the cccb (bbcb) tetraquark sector and it is around 400 MeV (400 MeV) below the lowest-lying bound state; therefore, bound states of the cccb and bbcb systems with narrow widths are not favored. Figure 6 shows the radial distribution functions for the J P = 0 + , 1 + and 2 + ground states of the cccb (upper panels) and bbcb (lower panels) systems. The interested reader can observe how the different quark-(anti-)quark rearrangements are taking over as the total spin is getting higher. For instance, in both cccb and bbcb tetraquark sectors, the (purple) dot-dashed curve which represents a quark-antiquark pair is the one that least extends in space for S = 0 and 1 indicating that quarkantiquark-pairs configuration could be important in such cases whereas the S = 2 prefers the diquark-antidiquarkpairs.
The typical mean-square radii, r 2 ij , of the studied cccb (upper rows) and bbcb (lower rows) tetraquarks are collected in Table XIII. For both tetraquark sectors, one can observe that the closest quark-(anti-)quark pair is the cb for S = 0 and S = 1 whereas is thecb for S = 2; indicating the change between quark-antiquark-pairs configuration and the one related with diquark-antidiquark pairs when the total spin grows. Finally, the mass meansquare radii for the J P = 0 + , 1 + and 2 + ground states are, respectively, 0.059 fm 2 , 0.064 fm 2 and 0.066 fm 2 for cccb system, and 0.035 fm 2 , 0.036 fm 2 and 0.037 fm 2 for bbcb system.
The cbcb tetraquark ground states. The cbcb system has no constraints from the Pauli principle. Therefore, looking at Table VII, there are four spin-color configurations for the J P C = 0 ++ ground state, another four in the case of the J P C = 1 +− , only two for the J P C = 1 ++ channel, and two more for the J P C = 2 ++ ground state. The predicted ground-state masses are listed in Table XIV and compared with the available results reported by other theoretical approaches [33,35,36]. Our results, with masses at around 12.5 GeV, are in reasonable agreement with those of Ref. [35]; however, both approaches predict figures which are slightly lower than the ones reported in Ref. [33]. The results published in Ref. [36] are systematically higher than others because the confinement term is ignored completely. As in the cases of cccb and bbcb tetraquarks, the J = 0 and the lowest J = 1  states are almost degenerate with even the J P C = 1 +− cbcb ground state located below the J P C = 0 ++ one. At this stage of our work, it could be necessary to remind that our calculation is parameter-free once the model is fitted to the meson and baryon sector.
The lowest S-wave meson-meson thresholds in the cbcb tetraquark sector are the η c η b (12429 MeV) for the J P C = 0 ++ channel, η c Υ(1S) (12467 MeV) for the J P C = 1 +− channel, and J/ψΥ(1S) (12563 MeV) for the J P C = 1 ++ and 2 ++ channels. We are pre-  dicting tetraquark ground state masses which lie 100 MeV above their lowest open-flavor S-wave mesonmeson threshold, i.e. they could appear as resonance candidates in the invariant mass of hteir corresponding meson-meson channel. It is worth emphasizing that the J P C = 1 ++ ground state is located at exactly, within theoretical uncertainty, its lowest S-wave meson-meson threshold. Figure 7 shows the relevant radial distribution functions for the J P C = 0 ++ , 1 +− , 1 ++ and 2 ++ ground states of the cbcb system. One can see that the radial distribution functions follow the same pattern for all ground states, i.e. the interquark distance of the bb-pair is the shortest one, it is followed by the one of the cc-pair, and finally the distances between the c-quark and either thē b-or b-quark are the same and the largest. However, the left bottom panel of Fig. 7 shows that the J P C = 1 ++ cbcb ground state is particularly different with respect the others, with a more extended radial distribution for the cb and cb pairs. In fact, such distribution is indicating that the interquark distance between the c-and b-quarks is of the order of 1 fm, whereas the cc-and bb-pairs are clustered within a distance of 0.5 fm or less. Therefore, it seems that the J P C = 1 ++ cbcb tetraquark ground state prefers to be in a meson-meson configuration.
In order to quantify our statements above, Table XIII shows the typical mean-square radii, r 2 ij , of the studied cbcb tetraquark ground states. In the case of the J P C = 1 ++ ground state, a mean-square radii larger than 1 fm 2 is obtained for the cb and cb pairs, whereas the others are ∼ 0.1 fm 2 . This may indicate, in contrast with the other cases studied, that cc and bb-pairs tend to form separated clusters within a distance of 1 fm. Furthermore, we report herein the mass mean-square radii for the J P C = 0 ++ , 1 +− , 1 ++ and 2 ++ ground states which are, respectively, 0.053 fm 2 , 0.055 fm 2 , 0.3 fm 2 and 0.076 fm 2 . One can see again that the J P C = 1 ++ cbcb ground state is twice more extended than the others.
Our results suggest that the J P C = 1 ++ ground state is remarkably different with respect to the other cases studied because, on one hand, its mass lies exactly at its lowest S-wave meson-meson threshold and, on the other hand, clusters of QQ-pairs (with Q either c or b) appear with a separation between them of the order of 1 fm. We believe that both features are intimately related and further studies shall be performed, which go beyond the scope of this manuscript.

IV. SUMMARY AND OUTLOOK
Fully-heavy tetraquarks have recently received considerable attention from experiment. The most significant example is the observation made by the LHCb collaboration of some enhancements in the J/ψ-pair invariant mass spectrum whose origin could be linked to hadron states consisting of four charm quarks. Moreover, one should expect that the searching of doubly hiddenbottom and -charm tetraquark states will probably become one of the most attractive experimental goals with the future running of BES III, Belle II, and LHC.
From the theoretical side, several approaches have been proposed to calculate the spectrum of tetraquark systems made up only by heavy quarks. Their main goal was to established theoretically the existence of fully-heavy tetraquarks with narrow widths, i.e. stable. Mixed results have been obtained with some theoretical studies claiming that some lowest-lying all-heavy tetraquarks could be located slightly lower than their respective thresholds of quarkonium pairs and others asserting that fully-heavy tetraquarks are located much higher in mass than their lowest possible meson-meson strong decay channel.
In order to contribute on a better understanding of the multiquark dynamics, we have used a diffusion Monte Carlo method to solve the many-body Schrödinger equation that describes the fully-heavy tetraquark systems. This approach allows to reduce the uncertainty of the numerical calculation, accounts for multi-particle correlations in the physical observables, and avoids the usual quark-clustering assumed in other theoretical techniques applied to the same problem. Moreover, Quantum Monte Carlo computations shall enable to scale the same boundstate problem to other multiquark systems such as pentaquarks, hexaquarks, etc. The used quark model Hamiltonian has a pairwise interaction which is the most general and accepted one: Coulomb + linear-confining + hyperfine spin-spin; therefore, our analysis should provide some rigorous statements about the mass location of the all-heavy tetraquark ground states. Note, too, that such conclusions are parameter-free because the model parameters were constrained by a simultaneous fit of 36 mesons and 53 baryons, with a range of agreement between theory and experiment around 10 − 20%, which can be taken as an estimation of the model uncertainty for fully-heavy tetraquarks.
The cccc, ccbb (bbcc) and bbbb lowest-lying states are located just in the middle of the mass ranges predicted by other theoretical approaches. All states appear above their corresponding meson-meson thresholds and thus the existence of stable cccc, ccbb (bbcc) and bbbb systems with very narrow widths is disfavored; nevertheless, this does not forbid to have resonances in these tetraquark sectors which can be experimentally observed in the near future. Interestingly too is the observation that there is a transition between quark-antiquark pairs and diquarkantidiquark ones when going from S = 0 to S = 2 in the QQQQ with all Q either c-or b-quark, but not in the ccbb (bbcc) sector. However, it is important to clarify that these states are compact ones, with sizes in the order of a typical hadron. At last, the J P C = 0 ++ cccc ground state is predicted to have a mass compatible with the enhancements observed by the LHCb collaboration.
Theoretical studies of the cccb, bbcb and cbcb systems are scarce because the complexity of the needed coupledchannels calculation. For the cccb and bbcb sectors, our results seem to indicate that the 0 + and 1 + ground states are almost degenerate, with the 2 + lowest-lying state located around 100 MeV above them. For the cbcb system, we predict small mass splittings between the studied bound-states and absolute mass values located in between those predicted by other theoretical works. Moreover, we found clear evidence that the J P C = 1 ++ cbcb ground state has a meson-meson molecular configuration which deserves to be investigated further.
Finally, the diffusion Monte Carlo method is, in principle, applicable to a wide range of related problems and open questions. For instance, some natural extensions of the work presented herein could be the analysis of excited states and the exploration of other multiquark systems such as pentaquarks, hexaquarks, etc. On the other hand, answering some complex questions related with few-and many-body hadron physics appears scientifically interesting such as what would be the binding energy per quark in a many-body bound-state system?, is there any limit in the number of quarks and antiquarks that a hadron could host?, or could other constituents play a role in the stability of exotic hadrons?