Diffusion Monte Carlo calculations of fully-heavy (multiquark) bound states
DDiffusion Monte Carlo calculations of fully-heavy (multiquark) bound states
M.C. Gordillo, ∗ F. De Soto, † and J. Segovia ‡ Departamento de Sistemas F´ısicos, Qu´ımicos y Naturales,Universidad Pablo de Olavide, E-41013 Sevilla, Spain (Dated: September 28, 2020)We use a diffusion Monte Carlo method to solve the many-body Schr¨odinger equation describingfully-heavy tetraquark systems. This approach allows to reduce the uncertainty of the numericalcalculation 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 sameproblem. 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. Thismeans that, in principle, our analysis should provide some rigorous statements about the masslocation of the all-heavy tetraquark ground states, which is particularly timely due to the veryrecent observation made by the LHCb collaboration of some enhancements in the invariant massspectra of
J/ψ -pairs. Our main results are: (i) the cc ¯ c ¯ c , cc ¯ b ¯ b ( bb ¯ c ¯ c ) and bb ¯ b ¯ b lowest-lying states arelocated well above their corresponding meson-meson thresholds; (ii) the J PC = 0 ++ cc ¯ c ¯ c groundstate with preferred quark-antiquark pair configurations is compatible with the enhancement(s)observed by the LHCb collaboration; (iii) our results for the cc ¯ c ¯ b and bb ¯ c ¯ b sectors seem to indicatethat the 0 + and 1 + ground states are almost degenerate with the 2 + located around 100 MeV abovethem; (iv) smaller mass splittings for the cb ¯ c ¯ b system are predicted, with absolute mass values inreasonable agreement with other theoretical works; (v) the 1 ++ cb ¯ c ¯ b tetraquark ground state lies atits lowest S -wave meson-meson threshold and it is compatible with a molecular configuration. I. INTRODUCTION
The
J/ψ signal was observed simultaneously atBrookhaven [1] and SLAC [2] in 1974; it was a heavyresonance with a surprisingly small decay width. Threeyears 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 boundstates of a heavy quark, Q , and its antiquark, ¯ Q , with Q either c - or b -quark, explained their narrow decay widthsand, in fact, was proved to be crucial to establish Quan-tum Chromodynamics (QCD) as the strong-interactionsector of the Standard Model of Particle Physics [5, 6].It is also worth mentioning that a related system, the c ¯ b bound state ( B + c ), has also been found in nature [7]; andthat the heaviest of the quarks, i.e. the top-quark, wasdiscovered in 1995 at Fermilab [8], with a mass around175 GeV and a large decay width that, due to weak in-teractions, forbids to form narrow t ¯ t resonances.Heavy quarkonia, viz. mesons containing only a heavyvalence quark-antiquark pair, opened the possibility touse a non-relativistic (NR) picture of QCD. They can in-deed be classified in terms of the quantum numbers of aNR bound state, and the spacings between radial, orbitaland spin excitations have a pattern similar to the onesobserved in positronium, an e + e − NR bound state wellstudied in Quantum Electrodynamics (QED) [9, 10]. Be-ing baryonic analogues of heavy quarkonia, triply-heavybaryons may provide a complementary window in the un- ∗ [email protected] † [email protected] ‡ [email protected] derstanding of the non-relativistic regime of QCD and thestrong interaction between heavy quarks, without takinginto account the usual light-quark complications. Onecan imagine that we could continue with the game ofincorporating valence heavy quarks and/or antiquarks,and thus form tetraquark, pentaquark, hexaquark, etc.,bound-state systems which would help to understand,at least, the generalization of the (heavy) quark-quarkstrong interaction to multi-body states. Note herethat composite hadrons with four and more, light orheavy, quarks were already conjectured in 1964 by Gell-Mann [11] and Zweig [12] as a side product of explainingthe observed spectrum of mesons and baryons from thequark model picture.It is a fact that many precise experimental results areavailable for conventional heavy quarkonia [13]. In addi-tion, tens of charmonium- and bottomonium-like states,the so-called XYZ states, which cannot fit the quarkmodel picture, have been identified along the last twodecades at B-factories (BaBar, Belle, and CLEO), τ -charm facilities (CLEO-c and BES) and hadron-hadroncolliders (CDF, D0, LHCb, ATLAS, and CMS). So far,there is no consensus about the nature of these exoticstates (see Refs. [14–17] for reviews of the experimentaland theoretical status of the subject). Their analysis andnew determinations will continue with the upgrade of ex-periments such as BES III [18], Belle II [19], and HL- andHE-LHC [20]. This will provide a sustained progress inthe field as well as the breadth and depth necessary fora vibrant heavy quark research environment.The ultimate aim of theory is to describe the proper-ties of the XYZ states from QCD’s first principles. How-ever, since this task is quite challenging, a more mod-est goal is to start with the development of QCD moti- a r X i v : . [ h e p - ph ] S e p vated phenomenological models that specify the coloredconstituents, how they are clustered and the forces be-tween them. In that line, simultaneously to the experi-mental measurements, theorists have been proposing forthe XYZ states different kinds of color-singlet clusters,made by quarks and gluons, which go beyond conven-tional mesons (quark-antiquark), baryons (three-quarks)and antibaryons (three-antiquarks); the most famous areglueballs, quark-gluon hybrids and multiquark systems(for a graphic picture of these kinds of hadrons see, forexample, Figs. 1, 6, and 7 of Ref. [16]). Related to the lastones, the first and best known exotic state is the X (3872),which was observed in 2003 as an extremely narrow peakin the B + → K + ( π + π − J/ψ ) channel and at exactly the¯ D D ∗ threshold [21, 22]. It is suspected to be a cn ¯ c ¯ n ( n = u or d quark) tetraquark state whose features resem-ble those of a molecule, but some experimental findingsseem to point out the existence in its wavefunction ofmore compact components such as diquark-antidiquarkand quark-antiquark [23–28].Finally, fully-heavy tetraquarks have recently receivedconsiderable attention, both experimentally and theoret-ically. On the experimental side, it is thought that all-heavy tetraquark states will be very easy to spot becausetheir masses should be far away from the typical mass re-gions populated by both conventional heavy mesons andthe XYZ states discovered until now. A search for deeplybound bb ¯ b ¯ b tetraquark states at the LHC was motivatedby Eichten et al. in Ref. [29], and it was carried out bythe LHCb collaboration [30] determining that no signif-icant excess is found in the µ + µ − Υ(1 S ) invariant-massdistribution. On the other hand, the LHCb collaborationhas recently released in Ref. [31] a study of the J/ψ -pairinvariant mass spectrum finding a narrow peak and abroad structure which could originate from hadron statesconsisting of four charm quarks.From the theoretical side, concerning the interactionbetween heavy quarks, chiral symmetry is explicitly bro-ken and thus meson-exchange forces cannot exist in afully heavy tetraquark system [32, 33], which would fa-vor the formation of genuine tetraquark configurationsrather than loosely bound hadronic molecules. This isinteresting by itself, but also it simplifies the kind ofquark–(anti-)quark interactions to be taken into account,justifying the proliferation of theoretical works.In the literature, we find fully-heavy tetraquark com-putations based on phenomenological mass formulae [34–36], QCD sum rules [32, 37–39], QCD motivated bagmodels [40], NR effective field theories [41, 42], potentialmodels [33, 43–55], non-perturbative functional meth-ods [56], and even some exploratory lattice-QCD calcu-lations [57]. Some works predict the existence of stable QQ ¯ Q ¯ Q ( Q = c or b ) bound states with masses slightlylower than the respective thresholds of quarkonium pairs(see, for instance, Refs. [32, 34, 35, 37, 38, 41, 42, 50].In contrast, there are other studies that predict no sta-ble cc ¯ c ¯ c and bb ¯ b ¯ b tetraquark bound states because theirmasses are larger than two-quarkonium thresholds (see, e.g. , Refs. [36, 43, 45, 47, 57]). To some extent, a bet-ter understanding of the mass locations of fully-heavytetraquark states would be desirable, if not crucial, forour comprehension of their underlying dynamics andtheir experimental hunting.The goal of the present study is to achieve the most general and accurate prediction for the ground states offully-heavy tetraquarks. In order to comply with thefirst feature, we are not going to assume any particu-lar clustering between the valence quarks (antiquarks),the interaction between them is the most simple and ac-cepted one: Coulomb + linear-confining + hyperfine spin-spin, and it will be implemented non-perturbativelly. The second feature, accurateness, is fulfilled using a diffu-sion Monte Carlo (DMC) technique for solving the many-body Schr¨odinger equation wich, in contrast with varia-tional methods, allows to reduce the uncertainty of thenumerical calculation at the percent level, since the sys-tematic one associated with the trial wave function iseliminated by the algorithm. Note, too, that it is pos-sible to disentangle the uncertainty inherent of many-body techniques from the theoretical one coming fromthe model.The manuscript is arranged as follows. In Sec. II thetheoretical framework is presented; we explain first theorigin, features and implementation of the computationalalgorithm and, later, the quark model Hamiltonian andhow their parameters are fixed. Section III is mostly de-voted to the analysis and discussion of our theoretical re-sults on fully-heavy tetraquarks; note here that we firstlystudy all-heavy mesons and baryons, comparing our re-sults with those available from variational methods inorder to confirm the validation of our approach. Finally,we summarize and give some prospects in Sec. IV.
II. THEORETICAL FRAMEWORKA. Quark model
The Hamiltonian which describes fully-heavy bound-state systems can be written as H = n-part. (cid:88) i =1 (cid:18) m i + (cid:126)p i m i (cid:19) − T CM + n-part. (cid:88) j>i =1 V ( (cid:126)r ij ) , (1)where m i is the quark mass, (cid:126)p i is the momentum ofthe quark, and T CM is the center-of-mass kinetic energy.Since chiral symmetry is explicitly broken in the heavy It is fair to notice that a similar calculation to the one presentedherein has been recently released in Ref. [58]; however, importantdifferences must be mentioned: (i) meson-meson and diquark-antidiquark clusters were assumed, (ii) the sextet-antisextetdiquark-antidiquark configuration was fully neglected, and (iii)the hyperfine spin-spin interaction was computed perturbatively. quark sector, the two-body potential, V ( (cid:126)r ij ), can be de-duced from the one-gluon exchange and confining inter-actions; i.e. V ( (cid:126)r ij ) = V OGE ( (cid:126)r ij ) + V CON ( (cid:126)r ij ) . (2)It is important to highlight that the potentials abovehave tensor and spin-orbit contributions, which shall beneglected in this work. This is because our main purposehere is to get a first, unified and non-perturbative reliabledescription of fully-heavy hadron ground states (from 2-to 4-quark systems) within the DMC method. Moreover,this kind of interactions appeared not to be essential fora global description of baryons [59], and beyond [54].The Coulomb and hyperfine terms are collected in theone-gluon exchange potential, and are given by V OGE ( (cid:126)r ij ) = 14 α s ( (cid:126)λ i · (cid:126)λ j ) (cid:34) r ij − π m i m j δ (3) ( (cid:126)r ij )( (cid:126)σ i · (cid:126)σ j ) (cid:35) , (3)where α s is the strong coupling constant fixed in a phe-nomenological way, (cid:126)λ are the SU (3)-color Gell-Mann ma-trices, and the Pauli spin matrices are denoted by (cid:126)σ . TheDirac delta function of the hyperfine term comes from theFermi-Breit approximation of the one-gluon exchange in-teraction. In order to perform non-perturbative calcula-tions, the δ (3) ( (cid:126)r ij ) is usually replaced by a smeared func-tion that, in our case, reads as follows δ (3) ( (cid:126)r ij ) → κ e − r ij /r π / r , (4)with κ a quark model parameter, and r = A (cid:16) m i m j m i + m j (cid:17) B a regulator which depends on the reduced mass of thequark–(anti-)quark pair.Confinement is one of the crucial aspects of the stronginteraction that is widely accepted and incorporated intoany QCD based model. Studies of QCD on a lattice havedemonstrated that multi-gluon exchanges produce an at-tractive linearly rising potential, which is proportional tothe distance between infinitely heavy quarks [60]. Thisphenomenological observation is usually modeled as V CON ( (cid:126)r ij ) = ( b r ij + ∆)( (cid:126)λ i · (cid:126)λ j ) , (5)where b is the confinement strength and ∆ is a globalconstant fixing the origin of energies.Table I shows the quark model parameters relevantfor this work. Note here that we are using the so-calledAL1 potential proposed by Silvestre-Brac and Semay inRef. [61], and applied extensively to the baryon sectorin Ref. [59]. It is worth emphasizing that the potentialcollects nicely the most important phenomenological fea-tures of QCD for heavy quarks, and that the parameterswere constrained by a simultaneous fit of 36 mesons and53 baryons with a remarkable agreement with data. TABLE I. Quark model parameters used herein and takenfrom AL1 potential in Refs. [59, 61].Quark masses m c (GeV) 1.836 m b (GeV) 5.227OGE α s κ A (GeV) B − B b (GeV ) 0.1653∆ (GeV) -0.8321 B. Computational algorithm
Quantum Monte Carlo (QMC) methods have beensuccessfully applied to many research areas but quan-tum chemistry and material science are the ones whichhave received more attention [62–64]. This is becauseQMC is a natural competitor of other methods wherethe uncorrelated or Hartree-Fock state does not pro-vide a good description of the many-body ground state.Other applications of QMC algorithms are solid-statephysics concerning the dynamics of condensed heliumsystems [65, 66], and studies on the properties of bothbosonic and fermionic ultracold quantum gases [67–70].Since nuclear Hamiltonians induce strong correlations,QMC methods have appeared to be very valuable in theunderstanding of nuclei and nucleonic matter. Varia-tional Monte Carlo (VMC) algorithms dealing with nu-clear interactions were introduced in the early 1980s [71].Afterwards, methods based on Green Function MonteCarlo (GFMC) burst into nuclear physics in the late1980s [72, 73], and were applied mostly to spin-isospin-dependent Hamiltonians. The GFMC technique is veryaccurate but becomes increasingly more complex whenmoving toward larger systems, being C the state-of-the-art studied one [74–76]. Diffusion Monte Carlo meth-ods [77] appear to be much more efficient at treating largesystems; however, there are unsolved issues when dealingwith spin-isospin dependent potentials.The application of QMC methods to hadron physicshas been scarce, basically because most known hadronsconsist on 2- and 3-body relativistic bound states. How-ever, the description from different perspectives of com-posite states with four and more, light or heavy, quarks(antiquarks) has recently received considerable attentionsince the discovery of the potentially non-relativistic, 4-quark, charmonium-like system X (3872). It is in thiscontext that QMC algorithms can contribute to shedsome light to the study of tetraquark, pentaquark, hex-aquark, etc., non-relativistic systems.Up to our knowledge, the first QMC study of mesonsand 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 com-pared reasonably well with those of the well-known Isgur-Karl’s quark model [80–83]. The Diffusion Monte Carlomethod has been applied recently to the fully-beautifultetraquark system in Ref. [58]; in particular, the authorscalculate the ground state energy of the J P C = 0 ++ bb ¯ b ¯ b system.The central idea behind the DMC method is to writethe Schr¨odinger equation for n -particles in imaginarytime ( (cid:126) = c = 1): − ∂ Ψ α (cid:48) ( R , t ) ∂t = ( H α (cid:48) α − E s )Ψ α ( R , t ) , (6)where E s is the usual energy shift used in DMC methods, R ≡ ( (cid:126)r , . . . , (cid:126)r n ) stands for the position of n particles and α denotes each possible spin-color channel, with givenquantum numbers, for the n -particles system. The func-tion Ψ α ( R , t ) can be expanded in terms of a complete setof the Hamiltonian’s eigenfunctions asΨ α ( R , t ) = (cid:88) i c i,α e − ( E i − E s ) t φ i,α ( R ) , (7)where the E i are the eigenvalues of the system’s Hamil-tonian operator, ˆ H . The ground state wave function, φ ,α ( R ), is obtained as the asymptotic solution of Eq. (6)when t → ∞ , as long as there is overlap betweenΨ α ( R , t = 0) and φ ,α ( R ), for any α -channel.A crucial feature of QMC methods is the use of im-portance sampling techniques [84], in order to reduce thestatistical fluctuations to a manageable level. Given aHamiltonian, H = H + V , of the form H α (cid:48) α = − ∇ R m δ α (cid:48) α + V α (cid:48) α ( R ) , (8)with m the mass of each particle which composes thebound-state system, if one works with the function f α ( R , t ) ≡ ψ ( R ) Ψ α ( R , t ) , (9)where ψ ( R ) is the time-independent trial function: ψ ( R ) = (cid:89) i A detailed discussion about the particular features ofour spectrum will be given in the following subsections.However, a comment is due here on the theoretical uncer-tainty of our results. There are two kind of theoretical er-rors: one is inherently connected to the statistical natureof the DMC algorithm and the other one is related with ashortcoming of the quark model approach and lies on theway to fix the model parameters. The statistical error isof the order of 1 MeV and negligible with respect the sys-tematic one related with the quark model. As mentionedabove, the set of model parameters are fitted to reproducea certain number of hadron observables within a deter-minate 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 us-ing them. As the range of agreement between theory andexperiment is around 10 − A. Fully-heavy mesons Four fundamental degrees of freedom at the quarklevel: space, spin, flavor and color are generally acceptedin QCD, and any hadron’s wave function must be ex-pressed as a product of these four terms: | ψ hadron (cid:105) = | φ r (cid:105)| χ s (cid:105) | χ f (cid:105) | χ c (cid:105) . (27)From now on, we shall drop out the trivial flavor wavefunction because only fully-heavy hadrons are consideredin this work and thus, due to the explicit breaking of chi-ral symmetry, the considered quark model Hamiltonianis blind to the heavy quark flavor.Concerning the color degree of freedom, any hadronstate must be a color singlet one, because no color chargeshas been observed in nature. The leading Fock state ofa fully-heavy meson is constituted by a quark, Q , andantiquark, ¯ Q , with Q either c or b -quark. Therefore, the SU (3) color wave function is constructed as follows: ⊗ = ⊕ c ⊗ ¯ c = c ⊕ c , (28)where the singlet state is the physically interesting onein this work, and it is given by the well know, symmetricexpression, | χ c (cid:105) meson = 1 √ (cid:16) | r ¯ r (cid:105) + | g ¯ g (cid:105) + | b ¯ b (cid:105) (cid:17) . (29)Since a meson is made by distinguishable particles (aquark and an antiquark), there is no restriction in thespin wave functions due to Pauli principle and thus all | χ S =0 ,S z =0 (cid:105) = 1 √ (cid:16) | ↑↓(cid:105) − | ↓↑(cid:105) (cid:17) , (30) | χ S =1 ,S z =+1 (cid:105) = | ↑↑(cid:105) , (31) | χ S =1 ,S z =0 (cid:105) = 1 √ (cid:16) | ↑↓(cid:105) + | ↓↑(cid:105) (cid:17) , (32) | χ S =1 ,S z = − (cid:105) = | ↓↓(cid:105) , (33)possibilities can be considered. The S = 0 state is anti-symmetric respect the particle exchange 1 ↔ 2, whereasthe S = 1 wave functions are all symmetric.One can guess that an excitation of a unit of angu-lar momentum costs an energy around 500 MeV. Thiseffect can be estimated from the experimental M ( L =1) − M ( L = 0) mass differences: f (1285) − ω (782) ≈ 499 MeV, a (1260) − ρ (770) ≈ 460 MeV, χ c (1 P ) − J/ψ ≈ χ b (1 P ) − Υ(1 S ) ≈ 432 MeV, but also in thebaryon sector as, for instance, N (1535) − N (940) ≈ 570 MeV. Therefore, we shall limit ourselves to ana-lyze states with orbital angular momentum equal to zero.An important consequence is that the space-wave func-tion will always represent an S -wave state, which is com-pletely symmetric.Once the (trial) wave function of the meson is con-structed, we follow the DMC algorithm explained in thesection above and obtain, as a proof of concept, themasses of the singlet and triplet S -wave ground statesof heavy quarkonia. Our results are shown in Table II,they compare fairly well with the experimental data de-spite the simplicity of the quark model; this also sup-ports our confidence on the mass prediction for fully-heavy tetraquark bound states. Our numerical techniqueis contrasted with the same calculation but using a vari-ational approach [54]. As one can see in Table II, thereare negligible differences between the two methods; as ex-pected, our values are always below the ones reported in TABLE II. The mass spectra of the heavy quarkonia in unitsof MeV. The refer to the nameand quantum numbers of the considered hadron, is our result, is the same calculation but using avariational method [54], and collects experimentaldata if exist. n S +1 L J J PC DMC VAR [54] EXP [13] η c S − + . . ± . J/ψ S −− . . ± . B c S − + . . ± . B ∗ c S −− η b S − + . . ± . S ) 1 S −− . . ± . Ref. [54]. This validates our technique against the vari-ational one, which should give an energy’s upper limitquite close to the real eigenvalue in the case of mesonsystems. As we shall see later, the differences betweenthe two numerical approaches will be larger as the num-ber of particles increases.The concept of radial distribution function can be ap-plied to multiquark systems and, in fact, can providevaluable information about the existence of interquarkcorrelations; in particular, 2-body correlations. If the n -particle wave function is defined as ψ ( (cid:126)r , . . . , (cid:126)r n ), wherespin, flavor and color degrees of freedom have been ig-nored for simplicity without lost of generalization, theprobability of finding particle 1 in position (cid:126)r , particle 2in position (cid:126)r , . . . , particle n in position (cid:126)r n is: P ( (cid:126)r , . . . , (cid:126)r n ) = ψ ∗ ( (cid:126)r , . . . , (cid:126)r n ) ψ ( (cid:126)r , . . . , (cid:126)r n ) , (34)and it is normalize to one, i.e. (cid:90) d(cid:126)r · · · d(cid:126)r n P ( (cid:126)r , . . . , (cid:126)r n ) . (35)Therefore, one can define ρ (2) ( (cid:126)r , (cid:126)r ) = (cid:90) d(cid:126)r · · · d(cid:126)r n P ( (cid:126)r , . . . , (cid:126)r n ) , (36)which expresses the probability of finding 2 particles inpositions (cid:126)r and (cid:126)r ; and the radial distribution functionas ρ ( r ) = 4 πr (cid:90) d (cid:126)R ρ (2) ( (cid:126)R + (cid:126)r, (cid:126)R ) (37)where r indicates now the distance between the two par-ticles considered.Figure 1 shows, for the studied mesons, the radial dis-tribution functions which are pure estimators calculatedwithin our DMC following Ref. [87, 88]. Among the fea-tures one can observe, the following are of particular in-terest: (i) the c ¯ c states are the most extended objects r ρ (r) [f m ] r [fm]c – c (S=0)c – c (S=1) r ρ (r) [f m ] r [fm]c – b (S=0)c – b (S=1) r ρ (r) [f m ] r [fm]b – b (S=0)b – b (S=1) FIG. 1. Radial distribution functions, r ρ (2) ( r ) with r = | (cid:126)r − (cid:126)r | the relative coordinate between the 2 particles, for thestudied charmonium ( left panel ), B c ( middle panel ), and bottomonium ( right panel ) states.TABLE III. Mean-square radii, (cid:104) r (cid:105) , of the studied quarko-nium systems, in units of fm .Meson (cid:104) r (cid:105) η c . J/ψ . B c . B ∗ c . η b . S ) 0 . and the system becomes more compact in going from the c ¯ c meson to the b ¯ b one, being the size of the b ¯ b systemalmost half of the c ¯ c -meson’s one; (ii) the S = 1 stateis slightly more extended than the S = 0 state becausethe 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, asexpected since the hyperfine mass splitting gets smaller.Finally, the typical mean-square radii, (cid:104) r (cid:105) , of thestudied quarkonium systems are shown in Table III. Theycompare nicely with the results reported by a well-knownnon-relativistic QCD effective field theory for quarko-nium [89, 90]. Moreover, looking at the table, we canconfirm that the interparticle distance of the b ¯ b systemis almost half of the c ¯ c one; and the B c ’s one is closer tothe bottomonium than to the charmonium. B. Fully-heavy baryons We turn now our attention to the study of all-heavyground state baryons and thus no orbital angular mo-mentum excitations must be considered, i.e. only S -wavesymmetric states are under scrutiny.Concerning the SU (3) color wave function, it is con- structed as follows: ⊗ ⊗ = ⊕ ⊕ ⊕ c ⊗ c ⊗ c = c ⊕ c ⊕ c ⊕ c , (38)where the colorless state is fully anti-symmetric and itscolor wave function is given by the textbook expression | χ c (cid:105) baryon = 1 √ (cid:16) | rgb (cid:105) + | gbr (cid:105) + | brg (cid:105)− | rbg (cid:105) − | grb (cid:105) − | bgr (cid:105) (cid:17) . (39)With three particles of spin 1 / 2, one can construct thefollowing spin states: ⊗ ⊗ = ⊕ ⊕ s ⊗ s ⊗ s = / S ⊕ / MS ⊕ / MA , (40)with | χ S =3 / ,S z =+3 / (cid:105) S = | ↑↑↑(cid:105) , (41) | χ S =1 / ,S z =+1 / (cid:105) MS = 1 √ (cid:16) | ↑↓↑(cid:105) + | ↓↑↑(cid:105)− | ↑↑↓(cid:105) (cid:17) , (42) | χ S =1 / ,S z =+1 / (cid:105) MA = 1 √ (cid:16) | ↑↓↑(cid:105) − | ↓↑↑(cid:105) (cid:17) , (43)examples of the spin wave functions used herein, withoutlost of generality. For the baryonsΩ ++ ccc , Ω + ccb , Ω cbb , Ω − bbb , (44)the spin wave function must comply Pauli statistics in thecase that quarks are the same. Namely, the S sym. = 3 / S TABLE IV. Masses, in MeV, of the ground states of fully-heavy baryons. The refer to thename and quantum numbers of the considered hadron, is our result, and is the same calculationbut using a variational method [59]. Note that, for compari-son, the masses presented here include three-body-force cor-rections using the value of the constant C reported in Table1 of Ref. [59].Baryon S +1 L J J P DMC Ref. [59]Ω ++ ccc S / / + + ccb S / / + cbb S / / + − bbb S / / + ), chargemean-square radii (in e fm ), and magnetic moments (in nu-clear magnetons) for all the baryons listed in Table IV.Observable Approach Ω ++ ccc Ω + ccb Ω cbb Ω − bbb (cid:104) R m (cid:105) DMC 0 . 069 0 . 040 0 . 028 0 . . 069 0 . 040 0 . 028 0 . (cid:104) R c (cid:105) DMC 0 . 138 0 . 096 0 . − . . 138 0 . 097 0 . − . (cid:104) µ (cid:105) DMC 1 . 023 0 . − . − . . 023 0 . − . − . wave function, which is completely symmetric, must beused for the Ω ccc and Ω bbb baryons; whereas the S sym. =1 / MS wave function, which is symmetric respect the twoparticles that are equal, must be used for the groundstates of Ω ccb and Ω cbb , because it is an allowed statewith a lower energy than the S = 3 / QQQ baryons ( Q = c or b ) in each allowedspin channel, and compares them with those obtained inRef. [59]. As one can see, there are negligible differencesbetween the two numerical approaches. Unfortunately,there is no experimental data to compare with; there-fore, we encourage the design of experimental set-upsat, for instance, the LHC@CERN facility able to detectthis kind of particles because the reward could be highand, as mentioned above, triply-heavy baryons are ide-ally suited to study QCD and, in particular, the heavyquark–(anti-)quark interaction, as it has been the casefor heavy quarkonia.In order to check further the capabilities of the DMCmethod to describe heavy baryons, we calculate severalstructural (static) properties to compare them with theresults obtained in Ref. [59]. In particular, the following ones: (cid:104) R m (cid:105) ≡ (cid:104) ψ hadron | (cid:88) i =1 m i M ( (cid:126)r i − (cid:126)R ) | ψ hadron (cid:105) , (45) (cid:104) R c (cid:105) ≡ (cid:104) ψ hadron | (cid:88) i =1 e i ( (cid:126)r i − (cid:126)R ) | ψ hadron (cid:105) , (46) (cid:104) µ (cid:105) ≡ (cid:104) ψ hadron | (cid:88) i =1 e i m i ( (cid:96) iz + 2 s iz ) | ψ hadron (cid:105) . (47)They are presented in Table V and one can see thatthe values obtained in Ref. [59] are perfectly reproduced.Therefore, as a proof of concept, our numerical techniquecan be applied to study not only eigenenergies but alsostructure properties of the hadrons under considerationwhen the momentum of the prove is less or equal thanthe typical hadron scale, 1 GeV.Equation (45) gives an idea of the physical size of thebaryon, where (cid:126)R = m (cid:126)r + m (cid:126)r + m (cid:126)r M , (48)is the center-of-mass coordinate, with M = m + m + m the total mass of the system. One can see in Ta-ble V that the results lie between 0 . 02 and 0 . 07 fm .This means that the spatial extension of such baryonsgoes from 0 . 15 to 0 . 26 fm while for the mesons we had √(cid:104) R m (cid:105) Q ¯ Q = √(cid:104) r (cid:105) / ∈ [0 . , . 2] fm, i.e. triply-heavybaryons are objects only slightly less compact than theirheavy quarkonia counterparts. This can be explained bythe fact that the color quark-quark interaction is halfas intense as the quark-antiquark one but, a priori , onewould expect a bigger effect and this could indicate thatsome kind of diquark correlations be at play.Equation (46) defines the charge mean-square radiithat can be deduced from the value of the baryon’selectric form factor at the photon point extracted fromlepton-baryon scattering experiments. A result thatmight be worth to emphasize is the electric charge ra-dius 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 canbe ignored, the magnetic moment operator is given byEq. (47), which is just the sum of the magnetic momentsof each quark, with orbital and spin contributions. Ourvalues for the positive-charged and neutral triply heavybaryons, Ω + ccb and Ω cbb , are respectively 0 . µ N and − . µ N . These can be compared with the values ofthe proton and neutron: 2 . µ N and − . µ N , collectedin Ref. [13]. Namely, the triply-heavy baryon partner ofthe nucleon have a magnetic moment 5 − 10 times smaller.Figure 2 shows the relevant radial distribution func-tions. Among the observed features, the following areof particular interest. In the cases of Ω ccc and Ω bbb baryons, the quark-quark probability distribution func-tion is broader than its quark-antiquark counterpart.There are two kinds of probability distributions in the Ω ccc r ρ (r) [f m ] r [fm] ccc – c (S=0) Ω ccb r ρ (r) [f m ] r [fm] cccb Ω cbb r ρ (r) [f m ] r [fm] cbbb Ω bbb r ρ (r) [f m ] r [fm] bbb – b (S=0) FIG. 2. Radial distribution functions for the studied Ω ccc ( upper-left panel ), Ω ccb ( upper-right panel ), Ω cbb ( bottom-left panel )and Ω bbb ( bottom-right panel ) baryons. The (pink) long-dashed curve plotted in the Ω ccc and Ω bbb panels is the same radialdistribution function but for the c ¯ c and b ¯ b mesons, respectively. cases of Ω ccb and Ω cbb baryons, one referring to the QQ -pair where the two heavy quarks are equal and the otherone when the QQ -pair is made with different species. Thesame pattern as in heavy quarkonia is observed in triply-heavy baryons where the QQ -pair seems to be more com-pact as the heavy quark mass is larger, and thus thiscould facilitate the appearance of strong ( cb )- and ( bb )-diquark correlations inside the Ω ccb and Ω cbb baryons,respectively.Finally, the typical mean-square radii, (cid:104) r ij (cid:105) , of thestudied triply-heavy baryons are collected in Table VI.Among other features, it highlights that the cb - and bb -pair are closer inside the Ω ccb and Ω cbb baryons, respec-tively.For completeness, we have calculated the J P = 3 / + lowest-lying states of the Ω ccb and Ω cbb baryons. Theirmasses are 8046 MeV and 11247 MeV, respectively. They TABLE VI. Mean-square radii, (cid:104) r ij (cid:105) , of the studied triply-heavy baryons, in units of fm .Baryon (cid:104) r cc (cid:105) (cid:104) r cb (cid:105) (cid:104) r bb (cid:105) Ω ++ ccc + ccb cbb - 0.117 0.073Ω − bbb - - 0.062 compare well with the variational ones reported inRef. [59]. We expect bigger differences when radial, an-gular and spin excitations will be compared; this goesbeyond the scope of the present manuscript but indicatesa possible next step to follow in the future. Finally, the0predicted mass splittings ∆ m = m (3 / + ) − m (1 / + ) are28 MeV and 32 MeV for the Ω ccb and Ω cbb baryons; theyare of the same order of magnitude than those collectedin Ref. [91], and references therein. C. Fully-heavy tetraquarks Diffusion Monte Carlo methods are designed to com-pute non-relativistic bound-states of a few- to many-particles system. This technique has been successfullyapplied to the field of nuclear physics studying light andmiddle nuclei, as well as objects of high nuclear den-sity such as neutron stars. However, as explained be-fore, applications of DMC to hadron physics are scarcebecause most hadrons were understood as relativisticbound states of few (two or three) light quarks (anti-quarks). Nowadays, many experimental signals point outthe existence of tetra- [31], penta- [92, 93] and even hexa-quark [94] systems, mostly in heavy quark sectors wherenon-relativistic dynamics could be thought as a good as-sumption, and thus Monte Carlo techniques are becom-ing attractive approaches to apply in hadron physics. Af-ter the application of the method to mesons and baryonsas a proof-of-concept, our next step is the study of all-heavy tetraquarks. Further studies of more complex mul-tiquark systems shall be performed in the future but theygo beyond the scope of this manuscript.Multiquark systems present richer color structuresthan mesons and baryons. Without assuming any kind ofclustering, the color algebra applied to tetraquark statesleads to the following irreducible representations c ⊗ c ⊗ ¯ c ⊗ ¯ c = (cid:16) × c (cid:17) ⊕ (cid:16) × c (cid:17) ⊕ (cid:16) × c (cid:17) ⊕ c , (49)where two color singlet states appear. They are usuallyknown as, respectively, the (¯3 c ⊗ c ) and (6 c ⊗ ¯6 c ) diquark-antidiquark configurations, i.e. c ⊗ c ⊗ ¯ c ⊗ ¯ c = (cid:16) ¯ c ⊕ c (cid:17) ⊗ (cid:16) c ⊕ ¯ c (cid:17) = (¯ c ⊗ c ) ⊕ ( c ⊗ ¯ c ) ⊕ . . . (50)Knowing that ¯ c diquark and c antidiquark represen-tations are antisymmetric under the transposition of thetwo particles: | ¯3 α (cid:105) = 1 √ (cid:15) αβγ | Q β (1) (cid:105)| Q γ (2) (cid:105) , (51a) | ,α (cid:105) = 1 √ (cid:15) αβγ | ¯ Q β (1) (cid:105)| ¯ Q γ (2) (cid:105) , (51b)where (cid:15) αβγ is the Levi-Civita tensor, with the greek let-ters going from 1 to 3; and that c diquark and ¯ c an- tidiquark are symmetric under the same transposition: | α (cid:105) = d αβγ | Q β (1) (cid:105)| Q γ (2) (cid:105) , (52a) | ¯6 ,α (cid:105) = d αβγ | ¯ Q β (1) (cid:105)| ¯ Q γ (2) (cid:105) ; (52b)one can build two orthogonal singlet color tetraquarkstates | ¯3 (cid:105) = 1 √ (cid:15) αβγ (cid:15) αλσ | Q β (1) (cid:105)| Q γ (2) (cid:105)| ¯ Q λ (3) (cid:105)| ¯ Q σ (4) (cid:105) , (53) | ¯6 (cid:105) = 1 √ d αβγ d αλσ | Q β (1) (cid:105)| Q γ (2) (cid:105)| ¯ Q λ (3) (cid:105)| ¯ Q σ (4) (cid:105) . (54)Their explicit expressions in terms of the familiar red , green and blue degrees-of-freedom, and without explicitlyclustering, can be written down as | ¯3 (cid:105) = 1 √ (cid:16) + | rg ¯ r ¯ g (cid:105) + | gr ¯ g ¯ r (cid:105) − | rg ¯ g ¯ r (cid:105)− | gr ¯ r ¯ g (cid:105) + | rb ¯ r ¯ b (cid:105) + | br ¯ b ¯ r (cid:105)− | rb ¯ b ¯ r (cid:105) − | br ¯ r ¯ b (cid:105) + | gb ¯ g ¯ b (cid:105) + | bg ¯ b ¯ g (cid:105) − | gb ¯ b ¯ g (cid:105) − | bg ¯ g ¯ b (cid:105) (cid:17) , (55)which is antisymmetric under the exchange of either bothquarks or both antiquarks, and | ¯6 (cid:105) = 1 √ (cid:104) + | rr ¯ r ¯ r (cid:105) + | gg ¯ g ¯ g (cid:105) + | bb ¯ b ¯ b (cid:105) + 12 (cid:16) + | rg ¯ r ¯ g (cid:105) + | gr ¯ g ¯ r (cid:105) + | rg ¯ g ¯ r (cid:105) + | gr ¯ r ¯ g (cid:105) + | rb ¯ r ¯ b (cid:105) + | br ¯ b ¯ r (cid:105) + | rb ¯ b ¯ r (cid:105) + | br ¯ r ¯ b (cid:105) + | gb ¯ g ¯ b (cid:105) + | bg ¯ b ¯ g (cid:105) + | gb ¯ b ¯ g (cid:105) + | bg ¯ g ¯ b (cid:105) (cid:17)(cid:105) , (56)which is symmetric under the exchange of either bothquarks or both antiquarks.It is also important to notice herein that the two colorstates defined above can be expressed in another set ofcolor representations | ¯3 (cid:105) = + (cid:114) | (cid:105) − (cid:114) | (cid:105) = − (cid:114) | (cid:105) + (cid:114) | (cid:105) , (57) | ¯6 (cid:105) = + (cid:114) | (cid:105) + (cid:114) | (cid:105) = + (cid:114) | (cid:105) + (cid:114) | (cid:105) , (58) The non-vanishing d αβγ = d αβγ constants are d = d = d = 1 and d = d = d = d = d = d = 1 / √ c ⊗ c , and color-hidden, c ⊗ c ,states.If we turn now our attention to the spin degree-of-freedom, the QQ ¯ Q ¯ Q ( Q = c or b ) system, made byfermions of spin 1 / 2, can have total spin S = 0, 1 and 2.There are two linearly independent S = 0 wave functionsthat can be written as | χ S =0 ,S z =0 (cid:105) SS = 1 √ (cid:16) + 2 | ↓↓↑↑(cid:105) + 2 | ↑↑↓↓(cid:105)− | ↓↑↑↓(cid:105) − | ↑↓↓↑(cid:105) + | ↓↑↓↑(cid:105) − | ↑↓↑↓(cid:105) (cid:17) , (59a) | χ S =0 ,S z =0 (cid:105) AA = 12 (cid:16) − | ↓↑↓↑(cid:105) − | ↓↑↑↓(cid:105)− | ↑↓↓↑(cid:105) + | ↑↓↑↓(cid:105) (cid:17) . (59b)They are, respectively, symmetric ( S ) and antisymmetric( A ) under the exchange of both quarks and both anti-quarks. The linearly independent S = 1 wave functionsare given by ( S z = S is assumed without lost of general-ity): | χ S =1 ,S z =+1 (cid:105) SA = 1 √ (cid:16) | ↑↑↑↓(cid:105) − | ↑↑↓↑(cid:105) (cid:17) , (60a) | χ S =1 ,S z =+1 (cid:105) AS = 1 √ (cid:16) | ↑↓↑↑(cid:105) − | ↓↑↑↑(cid:105) (cid:17) , (60b) | χ S =1 ,S z =+1 (cid:105) SS = 12 (cid:16) + | ↑↑↑↓(cid:105) + | ↑↑↓↑(cid:105)− | ↑↓↑↑(cid:105) − | ↓↑↑↑(cid:105) (cid:17) , (60c)which correspond to the symmetric − antisymmetric,antisymmetric − symmetric and symmetric − symmetricexchange of quarks and antiquarks. Finally, the S = 2spin wave function can be written as | χ S =2 ,S z =+2 (cid:105) SS = | ↑↑↑↑(cid:105) , (61)where S z = S is again assumed without lost of general-ity. Note that this spin wave function is fully symmetricunder the exchange of both quarks and both antiquarks.We shall focus our attention to the ground states of QQ ¯ Q ¯ Q , with Q either c - or b -quark in any possible com-bination. This implies that the space wave function istotally symmetric and, in order to fulfill Fermi-Diracstatistics, the wave-function configurations of Table VIImust be taken into account for each tetraquark sectorand J P ( C ) quantum numbers. We perform a coupled-channels calculation, based on the DMC algorithm ex-plained above, for those sectors shown in Table VII whichhave more than one spin-color configuration. It is worthhighlighting that the treatment of spin-dependent poten-tials and coupled-channels calculations are challengingfeatures for Monte Carlo methods which have avoidedtheir extended application to hadron physics, both areconsidered herein.Now, let us proceed to describe in detail our theoreticalfindings for each sector of fully-heavy tetraquarks. The cc¯c¯c tetraquark ground states. The massesof the S = 0, 1 and 2 ground states of cc ¯ c ¯ c tetraquarksare shown in Table VIII. We compare first our resultswith those obtained when solving the same Hamiltonianbut using a Rayleigh-Ritz variational method. One cansee that both approaches are compatible; however, thedifferences are larger than the ones shown in the mesonand baryon sectors because the variational methods be-gin to have difficulties as the number of particles grow.Nevertheless, the disagreement between both approachesis between 8 and 20 MeV, which is smaller than the usualuncertainty assigned to any quark model. Note, too, thatthe variational method only provides an upper limit ofthe eigenenergy and it depends on the goodness of thetrial wave function; DMC results depend less on such de-tails and, in principle, provides an exact estimate of theeigenenergy.Table VIII provides also a comparison with numerousworks that reported results on cc ¯ c ¯ c tetraquarks using alarge variety of techniques. It seems that our results6 . 35 GeV, 6 . 44 GeV and 6 . 47 GeV for the ground stateswith quantum numbers J P C = 0 ++ , 1 + − and 2 ++ , re-spectively, are located just in the middle of the rangescovered by all approaches which are [5 . − . 80] GeV,[6 . − . 90] GeV and [6 . − . 96] GeV. These massranges are in agreement with the recently observed struc-tures in the invariant mass distribution of J/ψ -pairs [31].The cc ¯ c ¯ c 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 ++ stateis unstable and can decay into η c η c and J/ψJ/ψ finalstates 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 themass threshold of J/ψJ/ψ , they can also easily decay intosuch charmonium pairs through quark rearrangement.Figure 3 shows for the cc ¯ c ¯ c ground states, the rele-vant radial distribution functions. It is very interestingto observe that quark-antiquark pairs are slightly closerthan quark-quark ones in the S = 0 ( J P C = 0 ++ ) statewhereas the situation gradually changes as spin increases.In any case, the cc and c ¯ c radial distribution functionsare very similar with a mean value close to ∼ . (cid:38) J P C = 0 ++ ,1 + − and 2 ++ cc ¯ c ¯ c ground states. They are, respectively,0 . 087 fm , 0 . 092 fm and 0 . 095 fm ; i.e. they are com-parable with those obtained for mesons and baryons.Moreover, Table IX shows the typical mean-square radiifor each (anti-)quark–quark pair indicating that no in-terquark distance is very different from the rest and thusall pairs play a similar role with distances of the order of ∼ . The bb¯b¯b tetraquark ground states. Let us nowturn 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 cc ¯ c ¯ c, bb ¯ b ¯ bcc ¯ b ¯ b ( bb ¯ c ¯ c ) 0 +(+) | χ S =0 (cid:105) SS ⊗ | ¯3 c c (cid:105) AA | χ S =0 (cid:105) AA ⊗ | c ¯6 c (cid:105) SS +( − ) | χ S =1 (cid:105) SS ⊗ | ¯3 c c (cid:105) AA +(+) | χ S =2 (cid:105) SS ⊗ | ¯3 c c (cid:105) AA cc ¯ c ¯ b , bb ¯ c ¯ b + | χ S =0 (cid:105) SS ⊗ | ¯3 c c (cid:105) AA | χ S =0 (cid:105) AA ⊗ | c ¯6 c (cid:105) SS + | χ S =1 (cid:105) SS ⊗ | ¯3 c c (cid:105) AA | χ S =1 (cid:105) SA ⊗ | ¯3 c c (cid:105) AA | χ S =1 (cid:105) AS ⊗ | c ¯6 c (cid:105) SS + | χ S =2 (cid:105) SS ⊗ | ¯3 c c (cid:105) AA cb ¯ c ¯ b +(+) | χ S =0 (cid:105) SS ⊗ | ¯3 c c (cid:105) AA | χ S =0 (cid:105) SS ⊗ | c ¯6 c (cid:105) SS | χ S =0 (cid:105) AA ⊗ | ¯3 c c (cid:105) AA | χ S =0 (cid:105) AA ⊗ | c ¯6 c (cid:105) SS +( − ) | χ S =1 (cid:105) SS ⊗ | ¯3 c c (cid:105) AA | χ S =1 (cid:105) SS ⊗ | c ¯6 c (cid:105) SS √ (cid:0) | χ S =1 (cid:105) SA − | χ S =1 (cid:105) AS (cid:1) ⊗ | ¯3 c c (cid:105) AA √ (cid:0) | χ S =1 (cid:105) SA − | χ S =1 (cid:105) AS (cid:1) ⊗ | c ¯6 c (cid:105) SS +(+) 1 √ (cid:0) | χ S =1 (cid:105) SA + | χ S =1 (cid:105) AS (cid:1) ⊗ | ¯3 c c (cid:105) AA √ (cid:0) | χ S =1 (cid:105) SA + | χ S =1 (cid:105) AS (cid:1) ⊗ | c ¯6 c (cid:105) SS +(+) | χ S =2 (cid:105) SS ⊗ | ¯3 c c (cid:105) AA | χ S =2 (cid:105) SS ⊗ | c ¯6 c (cid:105) SS TABLE VIII. Predicted masses, in MeV, for the cc ¯ c ¯ c system computed with the Diffusion Monte Carlo technique and comparedwith those obtained with the variational approach [54]. For completeness, we also compare our results with those of otherframeworks. J PC DMC VAR [54] [33] [36] [45] [32] [43] [95] [34] [46] [37, 38] [50] [35] [41]0 ++ < + − ++ (cid:104) r ij (cid:105) , of the studied fully-charmed tetraquarks, in units of fm . J PC (cid:104) r cc (cid:105) (cid:104) r c ¯ c (cid:105) (cid:104) r c ¯ c (cid:105) ++ + − ++ ground states of bb ¯ b ¯ b tetraquarks shown in Table X. Wecompare first our results with those obtained when solv-ing the same Hamiltonian but using a Rayleigh-Ritz vari-ational method. One can see that even larger differ-ences, between 30 and 50 MeV, are found in the bb ¯ b ¯ b sec-tor, with our results always below the ones reported inRef. [54], as expected. From the dynamical point of view,we would expect better agreement in tetraquark sectorswhere heavier quark masses are involved and thus theissue might be related with the oscillating parameters of the Gaussian basis used in Ref. [54] for the bb ¯ b ¯ b sector.Another possibility is that, as it will be shown later, in-terquark distances are smaller in these tetraquarks andthus multi-body correlations could play a more importantrole. Table X provides also a comparison with numerousworks that reported results on bb ¯ b ¯ b tetraquarks using alarge variety of theoretical techniques. Again, our re-sults 19 . 20 GeV, 19 . 28 GeV and 19 . 29 GeV for the groundstates with quantum numbers J P C = 0 ++ , 1 + − and 2 ++ ,respectively, are located just in the middle of the rangescovered by all approaches which are [18 . − . 16] GeV,[18 . − . 21] GeV and [18 . − . 24] GeV.Figure 4 shows the relevant radial distribution func-tions for the bb ¯ b ¯ b ground states. As in the case of fully-charm tetraquarks, the quark-antiquark pairs are closerthan quark-quark ones in the S = 0 ( J P C = 0 ++ )state. The situation, however, gradually changes asspin increases and thus the diquark (antidiquark) dis-tance seems to be slightly smaller for J P C = 1 + − and2 ++ bb ¯ b ¯ b tetraquarks. The theoretical fact that the J P C = 0 ++ bb ¯ b ¯ b ground state could prefer to be in3 S=0 cc – c – c r ρ (r) [f m ] r [fm] ccc – c S=1 cc – c – c r ρ (r) [f m ] r [fm] ccc – c S=2 cc – c – c r ρ (r) [f m ] r [fm] ccc – c FIG. 3. Relevant radial distribution functions for the studied cc ¯ c ¯ c tetraquark ground states.TABLE X. Predicted masses, in MeV, for the bb ¯ b ¯ b system computed with the Diffusion Monte Carlo technique and comparedwith those obtained with the variational approach [54]. For completeness, we also compare our results with those of otherframeworks. J PC DMC VAR [54] [33] [36] [37, 38] [34] [35] [41] [58] [32] [57] [41]0 ++ < + − ++ (cid:104) r ij (cid:105) , of the studied fully-bottom tetraquarks, in units of fm . J PC (cid:104) r bb (cid:105) (cid:104) r b ¯ b (cid:105) (cid:104) r b ¯ b (cid:105) ++ + − ++ quark-antiquark pairs, together with a predicted masswhich is 300 − 400 MeV above the η b η b and Υ(1 S )Υ(1 S )thresholds, could provide founded reasons to keep look-ing for (reasonably wide) bumps in the invariant mass ofbottomonium-pairs at the LHCb experiment.We report in Table XI the typical mean-square radiiof the studied bb ¯ b ¯ b tetraquarks. Such table highlightsagain the change of configuration picture from quark-antiquark pairs to diquark-antidiquark ones, when goingfrom S = 0 to S = 2 ground states of bb ¯ b ¯ b tetraquarks.Finally, the mass mean-square radii for the J P C = 0 ++ ,1 + − and 2 ++ ground states are, respectively, 0 . 026 fm ,0 . 028 fm and 0 . 029 fm . This indicates that these statesare very compact and far away from the picture of meson-meson molecules. The cc¯b¯b (bb¯c¯c) tetraquark ground states. Asimilar theoretical calculation must be performed forcomputing the ground states of the J P = 0 + , 1 + and 2 + cc ¯ b ¯ b ( bb ¯ c ¯ c ) tetraquarks. We obtain the masses: m ( J P = 0 + ) = 12865 MeV , (62a) m ( J P = 1 + ) = 12908 MeV , (62b) m ( J P = 2 + ) = 12926 MeV , (62c)which compare as follows m ( J P = 0 + ) = 12886 MeV , (63a) m ( J P = 1 + ) = 12924 MeV , (63b) m ( J P = 2 + ) = 12940 MeV , (63c)with the variational calculation [54]. As one can see,the differences between the two numerical methods arearound 20 MeV, and our values are always below thoseobtained by the variational method. Note, too, that ourresults are in reasonable agreement with the scarce onesreported by other theoretical methods [33]. Moreover,our predictions are above their lowest open-flavor decaychannels for about 240 − 280 MeV and thus they shouldappear as resonances with relatively large widths in theinvariant mass of B c B c , B c B ∗ c , or B ∗ c B ∗ c final states.An interesting feature distinguishes this sector withrespect the two already discussed, i.e. the cc ¯ c ¯ c and bb ¯ b ¯ b sectors. The expected pattern of interquark distancesfor cc , c ¯ b and ¯ b ¯ b (from larger to shorter) is conservedwhen changing the total spin S = 0, 1 and 2. An ex-ample of the output is drawn in Fig. 5 for the J P = 0 + cc ¯ b ¯ b ground state. We report, too, the mass mean-squareradii for the J P = 0 + , 1 + and 2 + ground states, which4 S=0 bb–b–b r ρ (r) [f m ] r [fm] bbb–b S=1 bb–b–b r ρ (r) [f m ] r [fm] bbb–b S=2 bb–b–b r ρ (r) [f m ] r [fm] bbb–b FIG. 4. Relevant radial distribution functions for the studied bb ¯ b ¯ b tetraquark ground states. S=0 cc–b–b r ρ (r) [f m ] r [fm] ccc–b – b – b FIG. 5. Radial distribution functions for the J P = 0 + ground state of cc ¯ b ¯ b tetraquark. The case bb ¯ c ¯ c is obviouslyequal. No further information is obtained for the J P = 1 + and 2 + cases. are 0 . 046 fm , 0 . 047 fm , 0 . 048 fm , respectively. Thesevalues lie between the ones reported above for the cc ¯ c ¯ c and bb ¯ b ¯ b tetraquarks The cc¯c¯b and bb¯c¯b tetraquark ground states. The cc ¯ c ¯ b , bb ¯ c ¯ b (and cb ¯ c ¯ b ) systems have not been stud-ied before with neither the Hamiltonian used herein northe variational method. Therefore, we cannot analyzein these cases the goodness of our numerical framework,diffusion Monte Carlo, with respect the variational one.We, however, shall compare our results with those ob-tained by other theoretical approaches, and numericaltechniques, but such works are scarce because the com-plexity of the coupled-channels calculation needed tostudy the cc ¯ c ¯ b , bb ¯ c ¯ b and cb ¯ c ¯ b tetraquarks.The cc ¯ c ¯ b and bb ¯ c ¯ b systems share some common fea-tures in terms of heavy quark symmetry and thus they TABLE XII. Predicted masses, in MeV, for the cc ¯ c ¯ b and bb ¯ c ¯ b systems computed with the diffusion Monte Carlo techniqueand compared with those obtained by Refs. [33, 36]. J P DMC [33] [36] DMC [33] [36]0 + + + will be discussed together. Table XII shows our spectrumof ground states with quantum numbers J P = 0 + , 1 + and2 + . As one can see, in both tetraquark sectors, the 0 + and 1 + states are almost degenerate, with the 2 + groundstate lying around 100 MeV above them. The same pic-ture is drawn by Ref. [36] while smaller mass splittingsare predicted in Ref. [33]. It is interesting to observe,too, that our absolute figures are in disagreement withthe other two cases reported in Table XII, with resultsof Ref. [36] much higher than ours and those of Ref. [33];the latter being of the same order of magnitude thanours. Another feature needed to be mentioned is thatthe lowest strong-decay threshold is η c B c ( η b B c ) for the cc ¯ c ¯ b ( bb ¯ c ¯ b ) tetraquark sector and it is around 400 MeV(400 MeV) below the lowest-lying bound state; therefore,bound states of the cc ¯ c ¯ b and bb ¯ c ¯ b systems with narrowwidths are not favored.Figure 6 shows the radial distribution functions forthe J P = 0 + , 1 + and 2 + ground states of the cc ¯ c ¯ b (up-per panels) and bb ¯ c ¯ b (lower panels) systems. The inter-ested reader can observe how the different quark–(anti-)quark rearrangements are taking over as the total spinis getting higher. For instance, in both cc ¯ c ¯ b and bb ¯ c ¯ b tetraquark sectors, the (purple) dot-dashed curve whichrepresents a quark-antiquark pair is the one that leastextends in space for S = 0 and 1 indicating that quark-antiquark-pairs configuration could be important in suchcases whereas the S = 2 prefers the diquark-antidiquark-pairs.The typical mean-square radii, (cid:104) r ij (cid:105) , of the studied cc ¯ c ¯ b S=0 cc – c–b r ρ (r) [f m ] r [fm] ccc – cc – b – c – b S=1 cc – c–b r ρ (r) [f m ] r [fm] ccc – cc – b – c – b S=2 cc – c–b r ρ (r) [f m ] r [fm] ccc – cc – b – c – b S=0 bb – c–b r ρ (r) [f m ] r [fm] bbb – cb–b – c – b S=1 bb – c–b r ρ (r) [f m ] r [fm] bbb – cb–b – c – b S=2 bb – c–b r ρ (r) [f m ] r [fm] bbb – cb–b – c – b FIG. 6. Relevant radial distribution functions for the J P = 0 + , 1 + and 2 + ground states of the cc ¯ c ¯ b (upper panels) and bb ¯ c ¯ b (lower panels) systems. (upper rows) and bb ¯ c ¯ b (lower rows) tetraquarks are col-lected in Table XIII. For both tetraquark sectors, one canobserve that the closest quark–(anti-)quark pair is the c ¯ b for S = 0 and S = 1 whereas is the ¯ c ¯ b for S = 2; in-dicating the change between quark-antiquark-pairs con-figuration and the one related with diquark-antidiquarkpairs when the total spin grows. Finally, the mass mean-square radii for the J P = 0 + , 1 + and 2 + ground statesare, respectively, 0 . 059 fm , 0 . 064 fm and 0 . 066 fm for cc ¯ c ¯ b system, and 0 . 035 fm , 0 . 036 fm and 0 . 037 fm for bb ¯ c ¯ b system. The cb¯c¯b tetraquark ground states. The cb ¯ c ¯ b sys-tem has no constraints from the Pauli principle. There-fore, looking at Table VII, there are four spin-color con-figurations for the J P C = 0 ++ ground state, another fourin 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 ++ groundstate. The predicted ground-state masses are listed in Ta-ble XIV and compared with the available results reportedby other theoretical approaches [33, 35, 36]. Our results,with masses at around 12 . cc ¯ c ¯ b and bb ¯ c ¯ b tetraquarks, the J = 0 and the lowest J = 1 TABLE XIII. Mean-square radii, (cid:104) r ij (cid:105) , of the studied cc ¯ c ¯ b (upper rows) and bb ¯ c ¯ b (lower rows) tetraquarks, in units offm . J P (cid:104) r cc (cid:105) (cid:104) r c ¯ c (cid:105) (cid:104) r c ¯ b (cid:105) (cid:104) r c ¯ b (cid:105) + + + J P (cid:104) r bb (cid:105) (cid:104) r b ¯ c (cid:105) (cid:104) r b ¯ b (cid:105) (cid:104) r c ¯ b (cid:105) + + + states are almost degenerate with even the J P C = 1 + − cb ¯ c ¯ b ground state located below the J P C = 0 ++ one. Atthis stage of our work, it could be necessary to remindthat our calculation is parameter-free once the model isfitted to the meson and baryon sector.The lowest S -wave meson-meson thresholds in the cb ¯ c ¯ b tetraquark sector are the η c η b (12429 MeV) forthe J P C = 0 ++ channel, η c Υ(1 S ) (12467 MeV) forthe J P C = 1 + − channel, and J/ψ Υ(1 S ) (12563 MeV)for the J P C = 1 ++ and 2 ++ channels. We are pre-6 TABLE XIV. Predicted masses, in MeV, for the cb ¯ c ¯ b systemcomputed with the Diffusion Monte Carlo technique and com-pared with those obtained by other theoretical frameworks. J PC DMC [33] [35] [36]0 ++ + − ++ ++ (cid:104) r ij (cid:105) , of the studied cb ¯ c ¯ b tetraquarks, in units of fm . J PC (cid:104) r cb (cid:105) (cid:104) r c ¯ c (cid:105) (cid:104) r c ¯ b (cid:105) (cid:104) r b ¯ b (cid:105) ++ + − ++ ++ dicting tetraquark ground state masses which lie (cid:46) 100 MeV above their lowest open-flavor S -wave meson-meson threshold, i.e. they could appear as resonancecandidates in the invariant mass of hteir correspondingmeson-meson channel. It is worth emphasizing that the J P C = 1 ++ ground state is located at exactly, withintheoretical uncertainty, its lowest S -wave meson-mesonthreshold.Figure 7 shows the relevant radial distribution func-tions for the J P C = 0 ++ , 1 + − , 1 ++ and 2 ++ groundstates of the cb ¯ c ¯ b system. One can see that the radial dis-tribution functions follow the same pattern for all groundstates, i.e. the interquark distance of the b ¯ b -pair is theshortest one, it is followed by the one of the c ¯ c -pair, andfinally the distances between the c -quark and either the¯ b - or b -quark are the same and the largest. However, theleft bottom panel of Fig. 7 shows that the J P C = 1 ++ cb ¯ c ¯ b ground state is particularly different with respect theothers, with a more extended radial distribution for the cb and c ¯ b pairs. In fact, such distribution is indicatingthat the interquark distance between the c - and b -quarksis of the order of 1 fm, whereas the c ¯ c - and b ¯ b -pairs areclustered within a distance of 0 . J P C = 1 ++ cb ¯ c ¯ b tetraquark ground stateprefers to be in a meson-meson configuration.In order to quantify our statements above, Table XIIIshows the typical mean-square radii, (cid:104) r ij (cid:105) , of the studied cb ¯ c ¯ b tetraquark ground states. In the case of the J P C =1 ++ ground state, a mean-square radii larger than 1 fm is obtained for the cb and c ¯ b pairs, whereas the others are ∼ . . This may indicate, in contrast with the othercases studied, that c ¯ c and b ¯ b -pairs tend to form separatedclusters within a distance of 1 fm. Furthermore, we reportherein the mass mean-square radii for the J P C = 0 ++ ,1 + − , 1 ++ and 2 ++ ground states which are, respectively, 0 . 053 fm , 0 . 055 fm , 0 . and 0 . 076 fm . One can seeagain that the J P C = 1 ++ cb ¯ c ¯ b ground state is twicemore extended than the others.Our results suggest that the J P C = 1 ++ ground stateis remarkably different with respect to the other casesstudied because, on one hand, its mass lies exactly at itslowest S -wave meson-meson threshold and, on the otherhand, clusters of Q ¯ Q -pairs (with Q either c or b ) appearwith a separation between them of the order of 1 fm. Webelieve that both features are intimately related and fur-ther studies shall be performed, which go beyond thescope of this manuscript. IV. SUMMARY AND OUTLOOK Fully-heavy tetraquarks have recently received con-siderable attention from experiment. The most signif-icant example is the observation made by the LHCbcollaboration of some enhancements in the J/ψ -pair in-variant mass spectrum whose origin could be linked tohadron states consisting of four charm quarks. Moreover,one should expect that the searching of doubly hidden-bottom and -charm tetraquark states will probably be-come one of the most attractive experimental goals withthe future running of BES III, Belle II, and LHC.From the theoretical side, several approaches havebeen proposed to calculate the spectrum of tetraquarksystems made up only by heavy quarks. Their maingoal was to established theoretically the existence offully-heavy tetraquarks with narrow widths, i.e. sta-ble. Mixed results have been obtained with some theo-retical studies claiming that some lowest-lying all-heavytetraquarks could be located slightly lower than their re-spective thresholds of quarkonium pairs and others as-serting that fully-heavy tetraquarks are located muchhigher in mass than their lowest possible meson-mesonstrong decay channel.In order to contribute on a better understanding of themultiquark dynamics, we have used a diffusion MonteCarlo method to solve the many-body Schr¨odinger equa-tion that describes the fully-heavy tetraquark systems.This approach allows to reduce the uncertainty of thenumerical calculation, accounts for multi-particle corre-lations in the physical observables, and avoids the usualquark-clustering assumed in other theoretical techniquesapplied to the same problem. Moreover, Quantum MonteCarlo computations shall enable to scale the same bound-state problem to other multiquark systems such as pen-taquarks, hexaquarks, etc.The used quark model Hamiltonian has a pairwiseinteraction which is the most general and acceptedone: Coulomb + linear-confining + hyperfine spin-spin;therefore, our analysis should provide some rigorousstatements about the mass location of the all-heavytetraquark ground states. Note, too, that such conclu-sions are parameter-free because the model parameterswere constrained by a simultaneous fit of 36 mesons and7 J PC =0 ++ cb – c–b r ρ (r) [f m ] r [fm] cbc–cc – b – b – b J PC =1 +- cb – c–b r ρ (r) [f m ] r [fm] cbc–cc – b – b – b J PC =1 ++ cb – c–b r ρ (r) [f m ] r [fm] cbc–cc – b – b – b J PC =2 ++ cc – c–b r ρ (r) [f m ] r [fm] cbc–cc – b – b – b FIG. 7. Relevant radial distribution functions for the J PC = 0 ++ , 1 + − , 1 ++ and 2 ++ ground states of the cb ¯ c ¯ b system. 53 baryons, with a range of agreement between theoryand experiment around 10 − cc ¯ c ¯ c , cc ¯ b ¯ b ( bb ¯ c ¯ c ) and bb ¯ b ¯ b lowest-lying states arelocated just in the middle of the mass ranges predictedby other theoretical approaches. All states appear abovetheir corresponding meson-meson thresholds and thusthe existence of stable cc ¯ c ¯ c , cc ¯ b ¯ b ( bb ¯ c ¯ c ) and bb ¯ b ¯ b systemswith very narrow widths is disfavored; nevertheless, thisdoes not forbid to have resonances in these tetraquarksectors which can be experimentally observed in the nearfuture. Interestingly too is the observation that there isa transition between quark-antiquark pairs and diquark-antidiquark ones when going from S = 0 to S = 2 inthe QQ ¯ Q ¯ Q with all Q either c - or b -quark, but not in the cc ¯ b ¯ b ( bb ¯ c ¯ c ) sector. However, it is important to clar-ify that these states are compact ones, with sizes in theorder of a typical hadron. At last, the J P C = 0 ++ cc ¯ c ¯ c ground state is predicted to have a mass compatible withthe enhancements observed by the LHCb collaboration.Theoretical studies of the cc ¯ c ¯ b , bb ¯ c ¯ b and cb ¯ c ¯ b systemsare scarce because the complexity of the needed coupled-channels calculation. For the cc ¯ c ¯ b and bb ¯ c ¯ b sectors, ourresults seem to indicate that the 0 + and 1 + ground statesare almost degenerate, with the 2 + lowest-lying state lo-cated around 100 MeV above them. For the cb ¯ c ¯ b sys-tem, we predict small mass splittings between the stud-ied bound-states and absolute mass values located in be-tween those predicted by other theoretical works. More-over, we found clear evidence that the J P C = 1 ++ cb ¯ c ¯ b ground state has a meson-meson molecular configuration8which deserves to be investigated further.Finally, the diffusion Monte Carlo method is, in prin-ciple, applicable to a wide range of related problems andopen questions. For instance, some natural extensionsof the work presented herein could be the analysis ofexcited states and the exploration of other multiquarksystems such as pentaquarks, hexaquarks, etc. On theother hand, answering some complex questions relatedwith few- and many-body hadron physics appears sci-entifically interesting such as what would be the bindingenergy per quark in a many-body bound-state system?, isthere any limit in the number of quarks and antiquarksthat a hadron could host?, or could other constituents play a role in the stability of exotic hadrons? ACKNOWLEDGMENTS We thank A. Lovato, J.M. Morgado, C.D. Roberts, J.Rodr´ıguez-Quintero for constructive discussions and con-tinuous support. This work has been partially funded bythe Ministerio Espa˜nol de Ciencia e Innovaci´on undergrant No. PID2019-107844GB-C22 and FIS2017-84114-C2-2-P; the Junta de Andaluc´ıa under contract No. Op-erativo FEDER Andaluc´ıa 2014-2020 UHU-1264517; butalso PAIDI FQM-205 and -370. The authors acknowl-edges, too, the use of the computer facilities of C3UPOat the Universidad Pablo de Olavide, de Sevilla. [1] J. Aubert et al. (E598), Phys. Rev. Lett. , 1404 (1974).[2] J. Augustin et al. (SLAC-SP-017), Phys. Rev. Lett. ,1406 (1974).[3] S. Herb, D. Hom, L. Lederman, J. Sens, H. Snyder, et al. ,Phys. Rev. Lett. , 252 (1977).[4] W. R. Innes, J. Appel, B. Brown, C. Brown, K. Ueno, et al. , Phys. Rev. Lett. , 1240 (1977).[5] T. Appelquist and H. Politzer, Phys. Rev. Lett. , 43(1975).[6] A. De Rujula and S. Glashow, Phys. Rev. Lett. , 46(1975).[7] F. Abe et al. (CDF), Phys. Rev. Lett. , 2432 (1998),arXiv:hep-ex/9805034.[8] F. Abe et al. (CDF), Phys. Rev. Lett. , 225 (1994),arXiv:hep-ex/9405005.[9] A. Ore and J. Powell, Phys. Rev. , 1696 (1949).[10] M. Deutsch, Phys. Rev. , 455 (1951).[11] M. Gell-Mann, Phys. Lett. , 214 (1964).[12] G. Zweig, CERN-TH-412, NP-14146 (1964).[13] M. Tanabashi et al. (Particle Data Group), Phys. Rev.D , 030001 (2018).[14] N. Brambilla et al. , Eur. Phys. J. C , 1534 (2011),arXiv:1010.5827 [hep-ph].[15] S. L. Olsen, Front. Phys. (Beijing) , 121 (2015),arXiv:1411.7738 [hep-ex].[16] S. L. Olsen, T. Skwarnicki, and D. Zieminska, Rev. Mod.Phys. , 015003 (2018), arXiv:1708.04012 [hep-ph].[17] N. Brambilla, S. Eidelman, C. Hanhart, A. Nefediev,C.-P. Shen, C. E. Thomas, A. Vairo, and C.-Z. Yuan,(2019), arXiv:1907.07583 [hep-ex].[18] D. Asner et al. , Int. J. Mod. Phys. A , S1 (2009),arXiv:0809.1869 [hep-ex].[19] A. Bevan et al. (BaBar, Belle), Eur. Phys. J. C , 3026(2014), arXiv:1406.6311 [hep-ex].[20] A. Cerri et al. , “Report from Working Group 4: Op-portunities in Flavour Physics at the HL-LHC and HE-LHC,” in Report on the Physics at the HL-LHC,and Per-spectives for the HE-LHC , Vol. 7, edited by A. Dainese,M. Mangano, A. B. Meyer, A. Nisati, G. Salam, andM. A. Vesterinen (2019) pp. 867–1158, arXiv:1812.07638[hep-ph].[21] S.-K. Choi et al. (Belle), Phys.Rev.Lett. , 262001(2003), arXiv:hep-ex/0309032. [22] B. Aubert et al. (Babar), Phys.Rev. D71 , 071103 (2005),arXiv:hep-ex/0406022.[23] D. Acosta et al. (CDF), Phys.Rev.Lett. , 072001(2004), arXiv:hep-ex/0312021 [hep-ex].[24] V. Abazov et al. , Phys.Rev.Lett. , 162002 (2004),arXiv:hep-ex/0405004 [hep-ex].[25] R. Aaij et al. (LHCb), Eur.Phys.J. C72 , 1972 (2012),arXiv:1112.5310 [hep-ex].[26] S. Chatrchyan et al. (CMS), JHEP , 154 (2013),arXiv:1302.3968 [hep-ex].[27] C. Collaboration (CMS), “Evidence for χ c (3872) inPbPb collisions and studies of its prompt production at √ s NN = 5 . 02 TeV,” (2019).[28] R. Aaij et al. (LHCb), (2020), arXiv:2005.13419 [hep-ex].[29] E. Eichten and Z. Liu, (2017), arXiv:1709.09605 [hep-ph].[30] R. Aaij et al. (LHCb), JHEP , 086 (2018),arXiv:1806.09707 [hep-ex].[31] R. Aaij et al. (LHCb), (2020), arXiv:2006.16957 [hep-ex].[32] W. Chen, H.-X. Chen, X. Liu, T. Steele, and S.-L. Zhu,Phys. Lett. B , 247 (2017), arXiv:1605.01647 [hep-ph].[33] M.-S. Liu, Q.-F. L, X.-H. Zhong, and Q. Zhao, Phys.Rev. D , 016006 (2019), arXiv:1901.02564 [hep-ph].[34] M. Karliner, S. Nussinov, and J. L. Rosner, Phys. Rev.D , 034011 (2017), arXiv:1611.00348 [hep-ph].[35] A. Berezhnoy, A. Luchinsky, and A. Novoselov, Phys.Rev. D , 034004 (2012), arXiv:1111.1867 [hep-ph].[36] J. Wu, Y.-R. Liu, K. Chen, X. Liu, and S.-L. Zhu, Phys.Rev. D , 094015 (2018), arXiv:1605.01134 [hep-ph].[37] Z.-G. Wang, Eur. Phys. J. C , 432 (2017),arXiv:1701.04285 [hep-ph].[38] Z.-G. Wang and Z.-Y. Di, Acta Phys. Polon. B , 1335(2019), arXiv:1807.08520 [hep-ph].[39] L. Reinders, H. Rubinstein, and S. Yazaki, Phys. Rept. , 1 (1985).[40] L. Heller and J. Tjon, Phys. Rev. D , 755 (1985).[41] M. N. Anwar, J. Ferretti, F.-K. Guo, E. Santopinto,and B.-S. Zou, Eur. Phys. J. C , 647 (2018),arXiv:1710.02540 [hep-ph].[42] A. Esposito and A. D. Polosa, Eur. Phys. J. C , 782(2018), arXiv:1807.06040 [hep-ph].[43] J. Ader, J. Richard, and P. Taxil, Phys. Rev. D , 2370 (1982).[44] S. Zouzou, B. Silvestre-Brac, C. Gignoux, andJ. Richard, Z. Phys. C , 457 (1986).[45] R. J. Lloyd and J. P. Vary, Phys. Rev. D , 014009(2004), arXiv:hep-ph/0311179.[46] N. Barnea, J. Vijande, and A. Valcarce, Phys. Rev. D , 054004 (2006), arXiv:hep-ph/0604010.[47] J.-M. Richard, A. Valcarce, and J. Vijande, Phys. Rev.C , 035211 (2018), arXiv:1803.06155 [hep-ph].[48] J.-M. Richard, A. Valcarce, and J. Vijande, Phys. Rev.D , 054019 (2017), arXiv:1703.00783 [hep-ph].[49] J. Vijande, A. Valcarce, and N. Barnea, Phys. Rev. D , 074010 (2009), arXiv:0903.2949 [hep-ph].[50] V. Debastiani and F. Navarra, Chin. Phys. C , 013105(2019), arXiv:1706.07553 [hep-ph].[51] X. Chen, Eur. Phys. J. A , 106 (2019),arXiv:1902.00008 [hep-ph].[52] X. Chen, Phys. Rev. D , 094009 (2019),arXiv:1908.08811 [hep-ph].[53] X. Chen, (2020), arXiv:2001.06755 [hep-ph].[54] G.-J. Wang, L. Meng, and S.-L. Zhu, Phys. Rev. D ,096013 (2019), arXiv:1907.05177 [hep-ph].[55] G. Yang, J. Ping, L. He, and Q. Wang, (2020),arXiv:2006.13756 [hep-ph].[56] M. A. Bedolla, J. Ferretti, C. D. Roberts, and E. San-topinto, (2019), arXiv:1911.00960 [hep-ph].[57] C. Hughes, E. Eichten, and C. T. H. Davies, Phys. Rev.D , 054505 (2018), arXiv:1710.03236 [hep-lat].[58] Y. Bai, S. Lu, and J. Osborne, Phys. Lett. B , 134930(2019), arXiv:1612.00012 [hep-ph].[59] B. Silvestre-Brac, Few Body Syst. , 1 (1996).[60] G. S. Bali, H. Neff, T. Duessel, T. Lippert, andK. Schilling (SESAM), Phys. Rev. D , 114513 (2005),arXiv:hep-lat/0505012.[61] C. Semay and B. Silvestre-Brac, Z. Phys. C , 271(1994).[62] B. Hammond, W. Lester, and P. Reynolds, Monte CarloMethods in ab Initio Quantum Chemistry (World Scien-tific, Singapore, 1994).[63] W. Foulkes, L. Mitas, R. Needs, and G. Rajagopal, Rev.Mod. Phys. , 33 (2001).[64] M. Nightingale and C. J. Umrigar, Quantum Monte CarloMethods in Physics and Chemistry (Springer, Vienna,2014).[65] K. E. Schmidt and D. M. Ceperley, “Monte carlo tech-niques for quantum fluids, solids and droplets,” in The Monte Carlo Method in Condensed Matter Physics ,edited by K. Binder (Springer Berlin Heidelberg, Berlin,Heidelberg, 1992) pp. 205–248.[66] D. Ceperley, Rev. Mod. Phys. , 279 (1995).[67] J. Carlson, S.-Y. Chang, V. Pandharipande, andK. Schmidt, Phys. Rev. Lett. , 050401 (2003),arXiv:physics/0303094.[68] S. Giorgini, L. P. Pitaevskii, and S. Stringari, Rev. Mod.Phys. , 1215 (2008), arXiv:0706.3360 [cond-mat.other]. [69] F. De Soto, C. Carbonell-Coronado, and M. C. Gordillo,Phys. Rev. A , 023633 (2014).[70] C. Carbonell-Coronado, F. D. Soto, and M. C. Gordillo,New Journal of Physics , 025015 (2016).[71] J. Lomnitz-Adler, V. Pandharipande, and R. Smith,Nucl. Phys. A , 399 (1981).[72] J. Carlson, Phys. Rev. C , 2026 (1987).[73] J. Carlson, Phys. Rev. C , 1879 (1988).[74] A. Lovato, S. Gandolfi, R. Butler, J. Carlson, E. Lusk,S. C. Pieper, and R. Schiavilla, Phys. Rev. Lett. ,092501 (2013), arXiv:1305.6959 [nucl-th].[75] A. Lovato, S. Gandolfi, J. Carlson, S. C. Pieper, andR. Schiavilla, Phys. Rev. Lett. , 182502 (2014),arXiv:1401.2605 [nucl-th].[76] A. Lovato, S. Gandolfi, J. Carlson, S. C. Pieper,and R. Schiavilla, Phys. Rev. C , 062501 (2015),arXiv:1501.01981 [nucl-th].[77] K. Schmidt and S. Fantoni, Phys. Lett. B , 99 (1999).[78] J. Carlson, J. B. Kogut, and V. Pandharipande, Phys.Rev. D , 233 (1983).[79] J. Carlson, J. Kogut, and V. Pandharipande, Phys. Rev.D , 2807 (1983).[80] N. Isgur and G. Karl, Phys. Rev. D , 4187 (1978).[81] N. Isgur and G. Karl, Phys. Rev. D , 2653 (1979),[Erratum: Phys.Rev.D 23, 817 (1981)].[82] N. Isgur and G. Karl, Phys. Rev. D , 1191 (1979).[83] S. Capstick and N. Isgur, AIP Conf. Proc. , 267(1985).[84] M. H. Kalos, D. Levesque, and L. Verlet, Phys. Rev. A , 2178 (1974).[85] J. Boronat and J. Casulleras, Phys. Rev. B , 8920(1994).[86] J. S´anchez-Baena, J. Boronat, and F. Mazzanti, Phys.Rev. A , 053632 (2018).[87] J. Casulleras and J. Boronat, Physical Review B ,36543661 (1995).[88] E. Krotscheck and J. Navarro, MicroscopicApproaches to Quantum Liquids in Con-fined Geometries , 074024(2013), arXiv:1302.3528 [hep-ph].[90] C. Peset, A. Pineda, and M. Stahlhofen, JHEP , 017(2016), arXiv:1511.08210 [hep-ph].[91] G. Yang, J. Ping, P. G. Ortega, and J. Segovia, Chin.Phys. C , 023102 (2020), arXiv:1904.10166 [hep-ph].[92] R. Aaij et al. (LHCb), Phys. Rev. Lett. , 072001(2015), arXiv:1507.03414 [hep-ex].[93] R. Aaij et al. (LHCb), Phys. Rev. Lett. , 222001(2019), arXiv:1904.03947 [hep-ex].[94] H. Clement, Prog. Part. Nucl. Phys. , 195 (2017),arXiv:1610.05591 [nucl-ex].[95] Y. Iwasaki, Prog. Theor. Phys.54