Masses and Magnetic Moments of Doubly Heavy Tetraquarks via Diffusion Monte Carlo Method

We present mass spectrum and magnetic moments of the $\bar{n}\bar{n}QQ$ states, where $n=u,d,s$ and $Q=c,b$. We solve four-body Schr\"odinger equation with a quark potential model by using diffusion Monte Carlo (DMC) method. The quark potential is based on the Coulomb, confinement and spin-spin interaction terms. We find mass and magnetic moment of the $T_{cc}^+$ state as $M_{T_{cc}^+}=3892 ~\text{MeV}$ and $\mu=0.28 \mu_N$, respectively. We also find the mass and magnetic moment of $T_{bb}^-$ as $M_{T_{bb}^-}=10338 ~\text{MeV}$ and $\mu=-0.32 \mu_N$, respectively. We find some bound state candidates of doubly heavy tetraquark systems with $I(J^P)=0(1)^+$ $nn \bar b \bar b$, $I(J^P)=0(0)^+$ $nn \bar c \bar b$, $I(J^P)=0(1)^+$ $nn \bar c \bar b$, and $I(J^P)=1/2(1)^+$ $ns \bar b \bar b$. We compare our results with other approaches in the literature.

After the breakthrough works of Gell-Mann [1] and Zweig [2] in 1964, hadrons can be classified into two groups in terms of their valence quarks and antiquarks: quark-antiquark (q q) configuration as mesons and three-quark configuration (qqq) as baryons.This classification was called as quark model and reproduced reasonably well many hadron observables such as baryon and meson spectrum, strong decays, magnetic moments of baryons, etc.However, the Quantum Chromodynamics (QCD) allows exotic hadrons beyond the scope of quark model.The existence of multiquark states is consistent with QCD.In fact, in Refs.[1,2], multiquark states beyond conventional ones (mesons and baryons) were conjectured as combinations of quarks.The era of conventional hadrons directed in a different route when the observation of χ c1 (3872) by the Belle Collaboration [3] announced in 2003.This state has a narrow width, decays into J/ψ2π and J/ψ3π and does not fall into the traditional charmonium spectra.Many unconventional states have been reported by Belle, BaBar, CLEO, BESIII, CDF, D0, LHCb, ATLAS, and CMS collaborations since 2003.These unconventional charmonium and bottomonium states are named as XYZ mesons.For a recent naming convention of these states, see Ref. [4].These unconventional states mostly decay into final states which include (cc) and (b b) pairs.At this point, heavy quarkonium physics becomes important for elucidating these states.The current status and prospects of the heavy quarkonium physics can be found in Refs.[5][6][7][8][9][10].
Hadron spectroscopy constitutes basis to test predictions of QCD.Hadron spectroscopy is related to masses and decays of the hadrons.The abundance and clear signals in the charmonium and bottomonium states prompt hadron spectroscopy.Charmonium and bottomonium states are heavy quarkonium states.Heavy quarkonium is a good laboratory for QCD: at high energies (small distances) where a perturbative calculations are possible and at low energies (large distances) nonperturbative effects take place [119].Since charm and bottom quarks are heavy, nonrelativistic approach of QCD can be used for these mesons.
There are various methods to study physical properties of exotic states.QCD motivated phenomenological models can be used to classify and specify the components (coloured constituents) and forces between them in order to study these exotic states.In that line, many important problems related to these states can be addressed by solving Schrödinger equation.Since these exotic states include at least more than three quarks, solving many-body Schrödinger equation is quite challenging.In general four-quark states are often taken as diquark-antidiquark states or meson-meson molecules.Such quark clusterings yield a two-body Schrödinger equation in the context of potential models.In this present work we rely on a different method for solving four-quark Schrödinger equation.In our model, none of the valence quarks will be assumed to be clustered together.
The interaction in the context of potential model for this system is the widely accepted and used one: Coulomb + linear confining + spin-spin term.We solve four-quark Schrödinger equation by using Diffusion Monte Carlo (DMC) method which takes into account full correlations between the constituents of the system and is capable of producing exactly the lowest eigenenergy since the effect of trial wave function is diminished by the method.
The order of the paper is as follows: We describe the potential model, numerical approach in Sec.II.We present and explain the mass spectrum and magnetic moments of the doubly heavy tetraquark states in Sec.III.The conclusion is reserved for Sec.IV.

A. Potential Model
Phenomenological models with a simple relativistic kinetic energy term and a scalar potential term which incorporates the so called linear confinement, and a term related to short distance which incorporates color-Coulomb interaction stemmed from QCD, give good results and descriptions of the observed spectra of both heavy and light quark mesons and baryons.
The Hamiltonian which describes doubly heavy tetraquark states reads as follows [120] with the potential where r 0 (m i , m j ) = A Ä 2mimj mi+mj ä −B , A and B are constant parameters, κ and κ ′ are parameters, r ij is the interquark distance |r i − r j |, σ i are the Pauli matrices and λi are Gell-Mann matrices.κ and κ ′ are related to the hyperfine interactions, just like the replacement of δ(⃗ r) with a smearing function in the Breit-Fermi approximation.There are four potentials referred to the p and r c : This potential was developed under the nonrelativistic quark model (NRQM) and used for exploratory studies.It compose of "Coulomb + linear" or "Coulomb + 2/3-power'" term and a strong but smooth hyperfine term.As mentioned in Ref. [120], potentials with r c ̸ = 0 give results slightly better than those with r c = 0, where r c = 0 means κ(r) and κ(r) ′ are constant parameters.Therefore we choose AL2 potential with the following parameter listed in Table I:

B. Diffusion Monte Carlo Method
Many important aspects of microscopic phenomena can be provided by the Schrödinger equation.Unfortunately, there are very few instances in which the Schrödinger equation can be solved exactly.For the rest which include more realistic cases, one needs to use approximations and numerical descriptions.Diffusion Monte Carlo (DMC) technique may help on this endeavour.
The imaginary-time Schrödinger equation (in ℏ = c = 1) can be written as with the Hamiltonian defined as In the above equations, R ≡ (⃗ r 1 , . . ., ⃗ r n ) represents the positions of n particles, E R is the shifted energy, and α denotes each possible spin-color channel of the system.The wave function in Eq. ( 4) can be expanded as where E i are the eigenvalues and Φ i,α (R) are the eigenfunctions.If E R is taken close to the ground state energy E 0 , i.e.E R ≃ E 0 , the wave function Ψ α (R, t) will approach to Φ 0,α (R) when sufficiently long time elapses, i.e. t → ∞.
An important feature of the quantum monte carlo (QMC) method is to use sampling technique to reduce statistical fluctuation [121].In this approach, we use the function where ψ(R) is the time-independent trial function which in principle should be chosen as close to Φ 0,α (R) as possible.Since ψ(R) is a trial function, it should be chosen according to the sake of the calculations.In this present work, we formulate it as [122] where the values of the adjustable constants a ij are set to minimize fluctuation.
With this new definition of the trial function (Eq.7), the Schrödinger equation turns out to be where is the drift force.The solution to Eq. ( 9) in path integral formalism reads as with Note herein that in the case of having more than one α-channel, we propagate the quantity [123] such as where the weight factor is introduced as The application of DMC method to hadron physics has been scarce.For an introduction of the DMC method see Ref. [124] and for applications to hadron physics see Refs.[122,[125][126][127][128][129][130].

C. Wave function
The total wave function of a tetraquark state can written as by taking into account possible color, spin and flavor configurations.The construction of the total wave function needs elaboration of Pauli principle.The possible spin wave functions for the (q 1 q 2 q3 q4 ) tetraquark states with the (qq Q Q) configuration are with the spin alignments where subscript outside the paranthesis and bracket denotes the spin of the hadron.The only color singlet configurations of q 1 q 2 q3 q4 tetraquark states are 6 c ⊗ 6c or 3c ⊗3 c color configurations with the respective color wave functions For the qq Q Q configuration, we can write Eq. ( 20) as and Concerning the color degree of freedom, the heavy quark is treated as an SU (3) singlet.
Armed with these considerations, we give the wave function bases in the notation of |(q 1 q 2 ) color spin ( Q3 Q4 ) color spin ⟩ total spin : The magnetic moment of a hadron arises from its charge distribution and angular momentum distribution.Not only it presents valuable information about the internal structure of the hadron but also the shape of the hadron.Furthermore, magnetic moment is an excellent probe for quark-gluon dynamics inside the hadron since it is the leading-order response of a bound state to a weak external magnetic moment [132].
The magnetic moment of a doubly-heavy tetraquark system can be written as where μm = Qi 2mi σz i is the magnetic moment operator, Qi is the electric charge operator of the i-th quark, σ z i is the third component (z-component) of Pauli matrix and ϕ IJ is the eigenvector of those states.

III. NUMERICAL RESULTS AND DISCUSSION
The DMC algorithm can be carried out to obtain the desired solution to Eq. ( 15) equation with the formalism described in Ref. [122].We follow the same algorithm.
The main idea can be described as follows: The kinetic and potential energy terms of the Schrödinger equation correspond to the diffusion term (source/sink) in the diffusion equation where u(r, t) is the density of the diffusing matter and D(u(r, t), r) is the diffusion term or coefficient.The Schrödinger equation with imaginary time is a kind of diffusion equation with source/sink mechanism.It can be solved by simulating random walks which are subject to birth/death processes resulting from the source/sink term.In the DMC algorithm, the wave function is represented by random walks of many particles.Each walker will walk randomly in finite small time steps and experience death or birth.After sufficiently long time, the diffusion of the walkers will approach the ground state wave function.
There are two theoretical errors in this model which have different sources.The first one is statistical related to DMC method itself and the second one is related to quark model parameters.We utilize Monte Carlo simulations comprising 50000 steps.We use 1 × 10 4 walkers in the simulations to implement the movement of each quark.Any implementation of DMC method leads to a point where the distribution of the walkers should be uncorrelated.To avoid this situation, we consider in the averages the values for every 500 steps and we drop the first 10 5 values to prevent any effect of chosen initial configuration.The resulting energy is averaged over the last 10000 steps with a standard deviation ≃ 1 MeV.
Before we present our results, we give possible corresponding thresholds in Table II regarding the states that are studied.We also present our results according to the strangeness quantum number.

System
I Thresholds The quark content for these tetraquark states are nncc, nn bb , and nnc b.In Table ??, our results for the mass and magnetic moment values are listed.
It can be seen from the table that, the mass splittings in the nncc system are 221 MeV for I = 1 and 238 MeV for I = 0.In the nn bb system, for I = 1 the mass splitting is 124 MeV and for I = 0 it is 286 MeV.For the mixed heavy quark system, nnc b, the mass splittings are 224 MeV for I = 1 and 276 MeV.One can also compare mass spectra with the corresponding thresholds which is shown in Fig. 1.
We can see from the figure that obtained mass values for I = 1 case are above corresponding thresholds.They should be resonance states.A special attention should be made on the I = 0 states.In the case of I = 0, one of the J P = 1 + states is high above than the DD * threshold.The other is just a few MeV above.The mass and magnetic moment of the T + cc are predicted to be M = 3892 MeV and µ = 0.28µ N , respectively.The obtained mass agree well with the observed T + cc mass, M = 3875 MeV.Ref. [60] calculated I(J P ) = 0(1 + ) T + cc magnetic moment as µ = 0.18µ N in compact configuration and µ = 0.13µ N in deuteron-like configuration.In Ref. [88], the magnetic moments of T + cc state are obtained in both diquark-antidiquark and molecular configurations via light-cone QCD sum rule formalism.They obtained magnetic moment as µ = 0.66 +0.34  −0.23 µ N with diquark-antidiquark configuration and µ = 0.43 +0. 23 −0.22 µ N with molecular configuration.Ref. [131] obtained magnetic moment of T + cc as µ = 0.88µ N by using MIT bag model.Ref. [133] calculated magnetic moment as µ = −0.09µN by using constituent quark model.Ref. [134] obtained a magnetic moment of µ = 0.732µ N by using heavy antiquark-diquark symmetry in the constituent quark model.The calculated results for the magnetic moment of T + cc in general are not consistent with each other.In Figure 2, mass spectrum of nn bb states can be seen.The results of I = 1 case are above corresponding thresholds.They should be resonance states.In the case of I = 0, the state which has a mass of 10624 MeV is 20 MeV above B B * threshold whereas the other is 266 MeV below.The obtained mass of this state in our model is M = 10338 MeV.This should be a deeply bound state.Considering magnetic moment of this state which is named as T − bb , we obtain µ = −0.32µN .Ref. [60] calculated I(J P ) = 0(1 + ) T − bb magnetic moment as µ = 0.64 µ N in compact configuration and as µ = 0.49µ N in deuteron-like configuration.In Ref. [92], magnetic moments of T bb states with J P = 1 + are obtained in light-cone QCD sum rule using molecular picture.The results follow as µ = (1.72 ± 0.67)µ N for B − B * − , µ = (1.38 ± 0.56)µ N for B 0 B * − , µ = (−0.44 ± 0.17)µ N for B − B * 0 , and µ = (−0.77± 0.28)µ N for B 0 B * 0 .Ref. [131]   In Figure 3, mass spectrum of nnc b states can be seen.All the results corresponding to I = 1 are above corresponding thresholds.These states should be resonances.We have bound state candidates in the I = 0 channel.The J P = 0 + state which has a mass of 7108 MeV is 39 MeV below D B threshold.In J P = 1 + states, the state which has a mass of 7084 is 108 MeV below D B * threshold and 204 MeV below D * B threshold.The state having a mass of 7204 MeV is 84 MeV below D * B threshold.We also obtain magnetic moments of I = 0 and I = 1 nnc b states.Focusing on I = 0 states, we obtain µ = 0.74µ N and µ = 0.45µ N for two low-lying J P = 1 + states.Ref. [60] obtained I(J P ) = 0(1 + ) T bc magnetic moment as µ = 0.66µ N in compact configuration and as µ = 0.59µ N in deuteron-like configuration.Our results agree with these results.The quark content for these tetraquark states are nscc and ns bb .We present mass spectrum and magnetic moment results in Table IV.As can be seen from table, in the nscc system the mass splitting is 285 MeV for J P = 0 + states and 145 MeV for J P = 1 + states.In the ns bb system, the mass splitting is 131 MeV for J P = 0 + states and 94 MeV for J P = 1 + states.In Figure 4, mass spectrum of nscc and ns bb tetraquarks can be seen.In the nscc systems, all obtained mass values are above corresponding thresholds.They should be resonance states.In the ns bb systems, the J P = 1 + state with a mass of M = 10684 is below corresponding thresholds B * s B/ B * Bs .This could be a bound state candidate.The J P = 2 + state is just 7 MeV above corresponding B * B * s state.This state seems to be a resonance state but also could be a shallow bound state candidate.Focusing on this state, in Ref. [93], J P = 1 + T QQqs and T QQss states are studied as diquark-antidiquark tetraquarks and magnetic moments are calculated with using heavy axial vector diquark-light scalar antidiquark type current denoted as J 1 , heavy scalar diquark-light axial vector antidiquark type current denoted as J 2 , heavy tensor diquark-light axial vector antidiquark type current denoted as J 3 and heavy vector diquark-light pseudotensor antidiquark type current denoted as J 4 .According to their work, J 1 gives µ = (−0.46± 0.07)µ N , J 2 gives µ = (−0.34± 0.12)µ N , J 3 gives µ = (−2.37 ± 0.50)µ N and J 4 gives µ = (−2.90± 0.36)µ N for T bbūs state.For T bb ds state, J 1 gives µ = (−0.59± 0.09)µ N , J 2 gives µ = (−0.68± 0.24)µ N , J 3 gives µ = (−1.37 ± 0.17)µ N and J 4 gives µ = (−1.45± 0.33)µ N .Our results are compatible with this work.In the same study considering T ccūs state, J 1 gives µ = (0.22 ± 0.07)µ N , J 2 gives µ = (0.18 ± 0.06)µ N , J 3 gives µ = (1.09± 0.36)µ N and J 4 gives µ = (1.04 ± 0.35)µ N .For T cc ds state, J 1 gives µ = (0.15 ± 0.05)µ N , J 2 gives µ = (0.36 ± 0.11)µ N , J The quark content for these tetraquark states are sscc and ss bb .We present mass spectrum and magnetic moment results in Table V.As can be seen from table, in the sscc system the mass splitting are 269 MeV for J P = 0 + states, 104 MeV for ss bb J P = 0 + states.In the ssc b systems, the mass splittings are 180 MeV for J P = 0 + states and 105 MeV for J P = 1 + states.In Figures 5, 6, and 7 mass values of sscc, ss bb , and ssc b states are located.As can be seen from figures all the states are above their corresponding thresholds, therefore they should be resonance states.Considering magnetic moments, Ref. [93] calculated magnetic moments of T ccss and T bbss states in light cone QCD sum rule formalism.The results are µ = (1.52 ± 0.46)µ N for J 3 type current and µ = (1.90 ± 0.56)µ N for J 4 type current.Our magnetic moment results of sscc are compatible with these results.For the T bbss state, they obtained µ = (−2.90± 0.23)µ N for J 3 type current and µ = (−3.20 ± 0.48)µ N for J 4 type current.Our results are not as compatible as in the previous case.For the wave functions in this work, the general ground state can be the mixing of ϕ 2 χ 3 and ϕ 1 χ 6 for the isotriplet tetraquarks with J P = 0 + .In the case of isosinglet tetraquarks with J P = 1 + , the wave function can be mixing of ϕ 2 χ 5 and ϕ 1 χ 4 for nn Q Q tetraquarks and mixing of ϕ 2 χ 2 , ϕ 2 χ 5 , ϕ 1 χ 4 for nn Q Q tetraquarks.In the nnc b tetraquarks, for I(J P ) = 1(0 + ) the wave function can be mixing of ϕ 2 χ 3 and ϕ 1 χ 6 , for I(J P ) = 0(0 + ) the wave function can be mixing of ϕ 1 χ 3 and ϕ 2 χ 6 , for I(J P ) = 1(1 + ) the wave function can be mixing of ϕ 2 χ 2 , ϕ 2 χ 4 , ϕ 1 χ 5 , and for I(J P ) = 0(1 + ) the wave function can be mixing of ϕ 2 χ 2 , ϕ 2 χ 5 , ϕ 1 χ 4 .
We compare our results with some of the works that studied same states.In Table VI we present our results with some other works conducted in different frameworks.As can be seen from table, our results are in agreement with the cited works.We also compare our mass results in Table VII with a recent work conducted in DMC [130].The agreement of the results can be seen.
Finally, we compare our results of magnetic moments with the available results in Table VIII.Our results are in general compatible with the cited work.

IV. CONCLUDING REMARKS
We employ a diffusion Monte Carlo (DMC) method to investigate systematically masses and magnetic moments of doubly heavy tetraquark states in constituent quark model.This method has some advantages compared to other methods; (i) reduces the uncertainty of the numerical procedure, (ii) accounts for multiparticle correlations, and (iii) avoids quark clustering in the relevant states.
Most QCD states arise as bound states or as resonances of multihadron configurations.Studying mass spectrum is an important aspect for the general view point of the related states.To obtain possible mass locations around the corresponding thresholds is not only crucial for understanding the underlying dynamics of the exotic hadrons and the nature of strong interactions of QCD but also useful for experimental researches for their existence.
In the nncc systems, the I(J P ) = 0(1 + ) nncc state with mass 3892 MeV is the candidate of the experimentally observed T + cc state.Its mass is only 16 MeV above DD * threshold.T + cc is argued to be a shallow bound state.Although the level of the uncertainty is less than 1 MeV, one cannot conclude this state is bound or not.Our results suggest I(J P ) = 0(1 + ) nncc could be a shallow bound state.All the other predicted states with I(J P ) = 0(1 + ) and I(J P ) = 1(1 + ) for nncc tetraquarks are well above corresponding thresholds.The magnetic moment presents information about the inner quark structure and can be used to distinguish the preferred configuration from predictions of different theoretical models.Our prediction of the magnetic moment for T + cc state is µ = 0.28µ N .In the nn bb systems, we find a deeply bound state candidate, namely I(J P ) = 0  regard to corresponding threshold.Some of the recent lattice QCD calculations [108,113,135] did not find any bound state indications of I(J P ) = 0(0 + ) and I(J P ) = 0(1 + ) udc b states but cannot rule out their existence.However, Refs.[105,136] support I(J P ) = 0(1 + ) udc b existence.A very recent lattice QCD study also found evidence for J = 0 shallow bound state and found hints for J = 1 broad resonance of udc b [137].
In the nscc and ns bb systems, we find all the nscc states above corresponding thresholds.They could be resonances.In the s threshold.This could be a shallow bound state candidate.In Ref. [113], a clear evidence for the existence of us bb with J P = 1 + is found.Lattice QCD studies present important evidences about the existence of the states.For the I = 0 T bb , I = 0, 1 T bc and T bbs tetraquark states, lattice QCD results support deeply bound state or very narrow resonances (see a recent review on lattice QCD results and references therein [138].Our results are in agreement with those predictions.
In the sscc, ss bb , and ssc b systems, all the states are above their corresponding thresholds.They all could be resonances.We have no bound state candidates in these systems.
The experimental progress on the multiquark states open a new era for our understanding of QCD.We hope that our results will be useful for not only for experimental studies but also other theoretical studies.

FIG. 2 .
FIG. 2. Mass spectrum of nn bb states with I = 0 (left), I = 1 (right) with corresponding J P values.Red color refers to threshold values whereas blue one refers to obtained results.

FIG. 3 .
FIG. 3. Mass spectrum of nnc b states with I = 0 (left), I = 1 (right) with corresponding J P values.Red color refers to threshold values whereas blue one refers to obtained results.

3 FIG. 4 .
FIG. 4. Mass spectrum of nscc (left) and ns bb (right) tetraquarks with corresponding J P values.Red color refers to threshold values whereas one refers to obtained results.

FIG. 5 .
FIG. 5. Mass spectrum of sscc tetraquark states.Red color refers to threshold values whereas blue one refers to obtained results.

FIG. 6 .
FIG.6.Mass spectrum of ss bb tetraquark states.Red color refers to threshold values whereas blue one refers to obtained results.

FIG. 7 .
FIG. 7. Mass spectrum of ssc b tetraquark states.Red color refers to threshold values whereas blue one refers to obtained results.

( 1 +
) nn bb state with mass 10338 MeV.It is just 266 MeV below B B * threshold.All the other nncc states are above their corresponding thresholds and could be resonances.The existence T − bb was conjectured long time ago and recent lattice QCD studies confirm the existence of ud bb tetraquark with I(J P ) = 0(1 + ) [102, 106, 107, 109].We obtain magnetic moment of T − bb as µ = −0.32µN .In the nnc b systems, we observe some bound state candidates.For I = 0, low-lying J P = 0 + nnc b has a mass of 7108 MeV and is 39 MeV below D B threshold.One of the J P = 1 + nnc b has a mass of 7084 MeV and is 108 below D B * threshold and 204 MeV below D * B * threshold.The J P = 1 + nnc b state with mass of 7204 MeV is 12 MeV above D B * threshold and 84 MeV below D * B * threshold.The two low-lying nnc b states with J P = 1 + could be a bound or loosely bound state with

TABLE II .
Corresponding thresholds for the doubly heavy tetraquark systems of this work.
obtained magnetic moment of T − [134].Mass spectrum of nncc states with I = 0 (left), I = 1 (right) and corresponding J P values.Red color refers to threshold values whereas blue one refers to obtained results.µ=0.18µ N whereas Ref.[134]obtained µ = −0.124µN for T − bb .As in the T + cc case, the calculated results for the magnetic moment of T − bb in general are not consistent with each other.

TABLE IV .
Mass spectrum and magnetic moments of doubly heavy nscc and ns bb tetraquarks, where n = u or d.Mass results are in units of MeV.Magnetic moments (in units of nuclear magneton µN , µp = 2.79285µN ) are organized in the order of I3 = 1/2, −1/2 for I = 1.

TABLE V .
Mass spectrum and magnetic moments of doubly heavy nscc and ns bb tetraquarks, where n = u or d.Mass results are in units of MeV.Magnetic moments are presented in units of nuclear magneton µN , µp = 2.79285µN .

TABLE VI .
[29]]rison of calculated mass (in MeV) with different calculations conducted in different models for double heavy tetraquarks.The masses before and after slash stand for color states splitting.Refs.[131],[66],[48]employ the CMI model and Ref.[29]employ relativistic quark model.

TABLE VII .
[130]rison of calculated mass (in MeV) with calculations conducted in DMC."NB" represents no bound solution.We only took the results which are bound state solutions from Ref.[130].

TABLE VIII .
Comparison of magnetic moments (in units of nuclear magneton µN , µp = 2.79285µN ).This could be a bound state candidate.The J P = 2 + state is 7 MeV above B * B * case of ns bb systems, the low-lying J P = 1 + state is 8 MeV below B * Bs threshold and 11 MeV below B * s B threshold.