Lattice QCD static potentials of the meson-meson and tetraquark systems computed with both quenched and full QCD
P. Bicudo M. Cardoso O. Oliveira P. J. Silva







bicudo@tecnico.ulisboa.pt



mjdcc@cftp.ist.utl.pt


CFTP, Instituto Superior Técnico, Universidade de Lisboa


orlando@fis.uc.pt



psilva@teor.fis.uc.pt 


CFisUC, Department of Physics, University of Coimbra, P-3004 516 Coimbra, Portugal

12.38.Gc, 11.15.Ha




We revisit the static potential for the Q Q Q̅Q̅ system using SU(3) lattice simulations, studying both the colour singlets groundstate and first excited state. We consider geometries where the two static quarks and the two anti-quarks are at the corners of rectangles of different sizes. 
We analyse the transition between a tetraquark system and a two meson system with a  two by two correlator matrix.
We compare the potentials computed with quenched QCD and with dynamical quarks. 
We also compare our simulations with the results of previous studies and
analyze quantitatively fits of our results with anzatse inspired in the string flip-flop model and in its possible colour excitations. 


Lattice QCD static potentials of the meson-meson and tetraquark systems

 computed with both quenched and full QCD
    P. J. Silva
    December 30, 2023
==================================================================================================================








§ INTRODUCTION




Our current understanding of strong interaction phenomenology, being the hadron spectrum or the form factors associated to
transitions between hadrons, relies on the description of the quark and gluon interaction within Quantum Chromodynamics. 
Despite the efforts of several decades, the non-perturbative nature of QCD still ensconce several properties of its fundamental particles. 
Indeed, we still do not understand the confinement mechanism,  which prevents the observation of free quarks and gluons in nature, 
and still do not have a satisfactory answer why the experimentally  <cit.> confirmed hadrons are composed of three valence quarks or a pair of quark and an anti-quark.


QCD is a gauge theory and physical observables should be gauge invariant objects. Gauge invariance implies
that only certain combinations of quarks and/or gluons can lead to observables particles. If one applies blindly such 
a simple rule, the observed hadrons are necessarily composite states involving multi-quarks and multi-gluon 
configurations. 
There is a priori no reason why states with other valence composition than mesons or baryons, called in general exotic states, should not be observed. 
Exotic states can be pure glue states (glueballs), multi-quark states (tetraquark, pentaquarks, etc) or hybrid states (mesons with a non-vanishing valence gluon content). Besides the hadron states compatible with the quark model, 
the particle data book  <cit.> also reports candidates for the different types of exotic states, see e.g. the reviews on pentaquarks and 
non-q q mesons. The masses of the experimental states listed as candidates to multi-quark/gluon hadrons cover the 
full range of energies of the particle spectrum. In particular the exotics with most observations are the tetraquarks.


In what concerns the experimental observation of exotic tetraquarks, the quarkonium sector of double-heavy tetraquarks including a QQ̅ pair is the most explored experimentally, see e.g. the recent reviews <cit.>. In particular, the charged Z_c^± and Z_b^± are crypto-exotic, but technically they can be regarded as essentially exotic tetraquarks if we neglect c c̅ or b b̅ annihilation.
There are two Z_b^± observed only by the collaboration BELLE at KEK
<cit.>, 
slightly above B  B̅^* and B^*  B̅^* thresholds,
 the Z_b(10610)^+ and Z_b(10650)^+.
Their nature is possibly different from the  two Z_c(3940)^± and  Z_c(4430)^±, whose mass is well above DD threshold <cit.>.
The Z_c^± has been observed with very high statistical significance
and has received a series of experimental observations by the BELLE collaboration <cit.>, the Cleo-C collaboration <cit.>, the
BESIII collaboration <cit.> and the LHCb
collaboration <cit.>. This family is possibly related to the closed-charm pentaquark recently observed at LHCb <cit.>. 
Notice that, using naïve Resonant Group Method calculations, in 2008, some of us predicted
<cit.>
a partial decay width to π J/ψ of the Z_c(4430)^- consistent with the recently observed experimental value<cit.>.

On the other hand, in what concerns lattice QCD simulations, the most promising exotic tetraquark sector is also double-heavy, but it has a pair of heavy quarks Q Q or antiquarks Q̅Q̅, and thus it differs from the quarkonium sector. Note that in lattice QCD, the study of exotics is presently even harder than in the laboratory, 
since the techniques and computer facilities necessary to study of resonances with many decay channels remain to be developed.

Lattice QCD searched for evidence of a large tetraquark component in the closed-charm Z_c(3940)^- candidate but this resonance is well above threshold, and Ref. 
<cit.>  concluded there is no robust lattice QCD evidence of a Z_c^± tetraquark resonance.

Lattice QCD also searched for the expected boundstate in light-light-antiheavy-antiheavy channels <cit.>.
Using dynamical quarks, the only heavy quark presently accessible to Lattice QCD simulations is the charm quark. No evidence for boundstates in this possible family of tetraquarks, say for a u d c̅c̅ was found.

Moreover the potentials between two mesons, each composed of a light quark and a static (or infinitely heavy) antiquark , have been computed in lattice QCD <cit.>.  A static antiquark constitutes a good approximation to a spin-averaged b̅ bottom antiquark. The potential between the two light-static mesons can then be used, with the Born-Oppenheimer approximation <cit.>, as a B-B potential,
where the higher order 1/m_b terms including the spin-tensor terms are neglected.
From the potential of the channel with larger attraction, which occurs in the Isospin=0 and Spin=0 quark-quark system, 
the possible boundstates of the heavy antiquarks have been investigated with quantum mechanics techniques.
Recently, this approach indeed found evidence for a tetraquark u d b̅b̅ boundstate <cit.>, while no boundstates have been found  for states where the heavy quarks are c̅b̅ or c̅c̅ (consistent with full lattice QCD computations <cit.>) or where the light quarks are s̅s̅ or c̅c̅ <cit.>. The b̅b̅ probability density in the only binding channel has also been computed in Ref.  <cit.>.





The quark models for tetraquarks with the most sophisticated description of confinement are the string flip-flop models. Clearly, tetraquarks are always coupled to meson-meson systems, and we must be able to address correctly the meson-meson interactions. The first quark models had confining two-body potentials proportional to the SU(3) colour Casimir invariant 
λ⃗_i ·λ⃗_j   V(r_ij) suggested by the One-Gluon-Exchange type of potential. However this would lead to an additional Van der Waals potential V_Van der Waals  =  V'(r)  r× T,
where T is a polarization tensor. The resulting Van der Waals 
<cit.>
force between mesons, or baryons would be extremely large and this is clearly not compatible with observations.
The string flip-flop potential for the meson-meson interaction was developed in Refs. 
<cit.>,
to solve the problem of the Van der Waals forces produced by the two-body
confining potentials.

The first considered string flip-flop potential was the one
minimizing the energy of the possible two different meson-meson configurations, say M_13  M_24 or M_14  M_23.
This removes the inter-meson potential, and thus solves the problem of the Van der Waals force.

An upgrade of the string flip-flop potential includes a third possible
configuration <cit.>, in the tetraquark channel, say T_12   ,   34,
where the four constituents are linked
by a connected string
<cit.>. 
The three confining string configurations differ in the strings linking the quarks and antiquarks, 
this is illustrated in 
Fig. <ref>. 
When the diquarks qq and q̅q̅ distances are small, the tetraquark configuration minimizes the string energy.
When the quark-antiquark pairs q q̅ and q q̅ are close, the meson-meson configuration minimizes
the string energy. 

With a triple string flip-flop potential, 
boundstates below the threshold for hadronic coupled channels have been found
<cit.>.

On the other hand, the string flip-flop potentials allow fully unitarized studies of resonances
<cit.>.
Analytical calculations with a double flip-flop harmonic oscillator potential,
<cit.>,
using the resonating group method again with a double flip-flop confining harmonic oscillator potential,
<cit.>,
anfd with the triple string flip-flop potential <cit.>
have already predicted resonances and boundstates.










So far. the theoretical and experimental interpretations of the observed states that can possibly be exotics is not clear crystal and, certainly, a better understanding of the colour force helps to elucidate our present view of the hadronic spectrum. 
For heavy quark systems its dynamics can be represented by a potential which, in general, is a function of the geometry of the hadrons, 
of the spin orientation of its components and of the quark flavours. In the limit of infinite quark mass one can compute the so-called
static potential using first principle lattice QCD techniques via the evaluation of Wilson loops. The static potential provides an important
input to the modelling of hadrons and it gives a simple realisation of the confinement mechanism. Moreover it can be applied to study tetraquarks Q  Q Q̅Q̅ with two heavy quarks and two heavy antiquarks, see for instance a Dyson-Schwinger study in Ref. 
<cit.>, at the intersection of the two sectors most studied experimentally and theoretically. 

The static potential has been computed using lattice QCD for mesons, tetraquark, pentaquarks and hybrid systems see 
<cit.>.
For a quark and an anti-quark system, the static potential V_QQ is a landmark calculation in lattice QCD and it is 
used to set the scale of the simulations. V_QQ has been computed both in the quenched theory and 
in full QCD with the lattice data being well described by a one-gluon exchange potential (a Coulomb like potential) at short distances and a linear rising function of the quark distances at large separations. The behaviour at large interquark distances
provides a nice explanation of the confinement mechanism. 
Moreover, for other hadronic systems and for large separations of its constituents a similar pattern of the corresponding
static potentials has been observed in lattice simulations, i.e. a linear rising potential which, once more, is a simple realisation of
quark confinement.

In the current work we revisit the static potential for tetraquarks using lattice simulations. The static potential for tetraquarks was 
computed for the gauge group SU(3) and in the quenched approximation in <cit.>. 
The hybrid potential defined and measured in <cit.> can also be viewed as a particular limit of the tetraquark potential. 
Herein, of all the possible geometries for the QQ Q̅Q̅ system we consider the case where quarks and anti-quarks are
at the corners of a rectangle, see Fig. <ref>, and recompute the static potential of the system both in the quenched 
approximation and in full QCD. We focus our analysis in the comparison of the quenched and full QCD and also in the transition between
a tetraquark system and a two meson system. Thus we go beyond the triple string flip-flop paradigm of Fig. <ref> and analyse, in the transition region, the mixing between the meson-meson and tetraquark string configurations. Moreover we explore not only the groundstate but also the first excited state.




The current work is organised as follows. In Sec. <ref> we discuss the possible colour structures for a 
QQ Q̅Q̅ system and introduce the Q Q̅ potentials used to compare the results of the static potentials
for the tetraquark. In Sec. <ref> we revisit the geometries used to compute the
static potentials and discuss the expected configurations at large separations. In Sec. <ref> the method used to 
evaluate the static potentials is described. In Sec. <ref> we report on the parameters used in the lattice simulations and how we set the 
scale of the simulations. The results for the static tetraquark potential for the two geometries are described in Sec <ref>.
In Sec. <ref> we resume and conclude. In the Appendix, the reader can find various tables with all our numerical results.







§ THE COLOR STRUCTURE OF A QQQ̅Q̅ SYSTEM 


The colour-spin-spatial wave function of a QQQ̅Q̅ system has multiple combinations, relevant for the computation of static potentials. In this section, we analyse the possible colour wave functions associated with a tetraquark system.

The quarks belong to the fundamental 3 representation of SU(3), while anti-quarks are in a 3 representation of the group. 
The  space built from the direct product 3⊗3⊗3⊗3 includes two independent colour singlet states. 

In a QQQ̅Q̅ system, quarks and anti-quarks can combine into colour singlet meson-like states, 
leading naturally to the two meson states,

    |1_131_24⟩   =   1/3δ_ikδ_jl|Q_iQ_jQ̅_kQ̅_l⟩ ,
    
    
    |1_141_23⟩   =   1/3δ_ilδ_jk|Q_iQ_jQ̅_kQ̅_l⟩ ,

where only the colour indices are written explicitly and 1_ij refers to the meson-like colour singlet state built
combining quark i and anti-quark j. The two colour singlet states in Eq. (<ref>) are not orthogonal to each other and a straightforward algebra gives,

    ⟨1_131_24|1_141_23⟩=1/3 .


Moreover, a quark and anti-quark pair, besides a colour singlet state, can also form a colour octet state. With two colour octets it is again possible to build a colour singlet state.
For the QQQ̅Q̅  system the colour singlet states
built from the octets read,

    |8_138_24⟩  =  1/4√(2)λ_ik^aλ_jl^a|Q_iQ_jQ̅_kQ̅_l⟩ ,
    
    
    |8_148_23⟩  =  1/4√(2)λ_il^aλ_jk^a|Q_iQ_jQ̅_kQ̅_l⟩ ,

where the factors comply with the normalization condition,

    ⟨8_138_24  | 8_138_24⟩ = 
    ⟨8_148_23  |  8_148_23⟩ = 1  .

The colour octet-octet states in Eq. (<ref>) can be written in terms of the meson-meson states defined in Eq. (<ref>),

    |8_138_24⟩   =   3|1_141_23⟩-|1_131_24⟩/2√(2) ,
    
    
    |8_148_23⟩   =   3|1_131_24⟩-|1_141_23⟩/2√(2) .

A simple calculation shows that the colour octet states (<ref>) are not orthogonal (in colour space) to each other.
However, each of the octet-octet states is orthogonal to the corresponding meson-meson state, i.e.

    ⟨1_131_24|8_138_24⟩   =     0  ,
    
    ⟨1_141_23|8_148_23⟩   =     0  .


The states in Eqs. (<ref>)  and (<ref>) do not represent all the possible colour singlet states that can be associated to a 
QQQ̅Q̅ system.   We can also consider diquark-antidiquark configurations.
For the group SU(3) it follows that  3⊗3=3̅⊕6, 
3̅⊗3̅=3⊕6̅ and the two colour singlet states belong to the space spanned by
3⊗3̅ and 6⊗6̅,

    |3̅_123_34⟩   =   1/2√(3) ϵ_ijm ϵ_klm  |Q_iQ_jQ̅_kQ̅_l⟩
       =   √(3/4)(|1_131_24⟩-|1_141_23⟩)  ,
      
    |6_126̅_34⟩   =   √(3/8)(|1_131_24⟩+|1_141_23⟩)
     .

The states in Eqs. (<ref>) and (<ref>) are orthogonal to each other in colour space,  i.e.
⟨3̅_123_34|6_126̅_34⟩  =  0. Furthermore,
they are eigenstates of the exchange operators of quarks or anti-quarks, and verify the following relations,

    P_12|3̅_123_34⟩   = P_34|3̅_123_34⟩  =    -  
    |3̅_123_34⟩ ,
    
    
    P_12|6_126̅_34⟩   = P_34|6_126̅_34⟩  =     +  
    |6_126̅_34⟩  ,

where P_ij is the exchange operator of (anti)quark i with (anti)quark j. 
Eqs.  (<ref>) and (<ref>) can be inverted, giving,

    |1_131_24⟩  =  √(2/3)|6_126̅_34⟩+1/√(3)|3̅_123_34⟩ ,
    
    
    |1_141_23⟩  =  √(2/3)|6_126̅_34⟩-1/√(3)|3̅_123_34⟩ ,

which shows that the meson-meson states of Eq. (<ref>) are not eigenstates of the quark and of the anti-quark
exchange operators P_12 and P_34.


The static potential V for a QQQ̅Q̅ system is a complicated object which may involve, two, three and four body interactions.
In general, V also depends on the allowed quantum numbers of the constituents of the multiquark state. The static potential should allow, when combined with quantum mechanics, for the groundstates to be the ones of Fig. <ref>.
For example, the static potential should allow for the formation of two-meson states when the quark-anti-quark distances are small compared to the quark-quark and anti-quark-anti-quark distances, or possibly for the formation of a tetraquark at other particular distances.

As an approximate model to understand the results of the lattice simulations for the static potential in terms of overlaps with the various colour singlets, one can consider the two-body potential given by the Casimir scaling,

    V_CS = ∑_i<j C_ij  V_M ,

where V_M is the mesonic static Q Q̅ potential
with C_ij = λ_i^a·λ_j^a/-16/3 and compare the results of the simulations with the one of any of the colour singlet states and the Casimir potential given by,

    V_Ψ = ⟨Ψ | V_CS | Ψ⟩ .

Note, for a two body system, the one gluon exchange predicts a static potential proportional to λ^a_i ·λ^a_j.













The expectation values ⟨Ψ | C_ij | Ψ⟩ for the possible colour singlet states associated to the QQQ̅Q̅ system
are reported in Table <ref>. These numbers are important to obtain a qualitative insight into the result of the simulations. 
For instance, if for a given state C_ij < 0, we don't expect that the lattice result would give us a strong attraction between the 
particles i and j and, therefore, one can expected significant deviations of the static potential relative to the potential V_CS associated to
the corresponding colour singlet state.


Moreover we consider as well the first excitation of the QQQ̅Q̅, which also depends in the particular distances of the system.
Based in the orthogonality conditions and in a crude Casimir scaling where V_M would be a spatial independent potential,
we would expect the pairs of colour singlet states,
(|1_13 1_24⟩, |8_13 8_24⟩) ,
(|1_14 1_23⟩, |8_14 8_23⟩) and 
(| 3_12 3_34⟩, |6_136_24⟩)
to form possible (groundstate, first excited state) pairs. This already goes beyond the simple paradigm of Fig. <ref>.


Nevertheless, Eq. (<ref>) is clearly an approximation, and our aim is to compute more rigorous potentials.
Previous lattice studies 
<cit.>
show that the static potential for a tetraquark system is not described entirely by a function proportional
to this potential. 
An example of such a kind of potentials is the two-meson potential,

    V_33   =   ⟨1_131_24 | V_CS | 1_131_24⟩ ,
    
       =   
    V_M(r_13) + V_M(r_24)  ,

which we expect to saturate the ground state when the quark-quark and anti-quark-anti-quark distances are large.

















§ GEOMETRICAL SETUP 











We aim to measure the static potential for the QQQ̅Q̅ system but also to investigate the
transition between the  tetraquark and a two meson state,  and the transition between the two two-meson states. This computation within lattice QCD simulations requires choosing a particular geometrical setup of the 
quark system under investigation. In principle, one could choose any of the available geometrical configurations allowed 
by the hypercubic lattice. In order to study in detail the transitions between the different states, in the current work we opt for restricting our study to the case where the four particles are at the corners of a 
rectangle and look at 

















two particular alignments.  In the so-called parallel alignment, see Fig. <ref> (left), the two quarks (anti-quarks) are at adjacent
corners of the rectangle. In the anti-parallel alignment, see Fig. <ref> (right), the quarks (anti-quarks) are at the opposite corners of the 
rectangle. 




















 §.§ Parallel Alignment of Quarks


For this geometry, where the two quarks are at neighbour corners of the rectangle, we can describe the system via
the intra-diquark distances,

    r_12=|𝐱_1-𝐱_2|=|𝐱_3-𝐱_4|  ,

and the inter-diquark distances,

    r_13=|𝐱_1-𝐱_3|=|𝐱_2-𝐱_4|  .

Note that for both cases the second equality  holds only due to the particular geometrical configuration considered.

If one assumes that quarks are confined within colourless states, this geometrical setup has two limits which allow to study the transition
between a tetraquark state and a two meson system. Indeed, when r_12≪ r_13 one expects the ground state of the 
QQQ̅Q̅  system to be that of a tetraquark,  while for the opposite case, i.e. for r_13≪ r_12, one expects the system, i.e.
its potential, to behave as a two meson system.





For this geometrical setup, in the evaluation of the static potential we consider the basis of operators shown in Fig. <ref>. 
They are associated with a tetraquark operator (left in the figure) and a two-meson operator (right in the figure), the two ground state configurations 
expected for this particular geometry.








 §.§ Anti-parallel Alignment of Quarks


For the anti-parallel alignment of quarks described in Fig. <ref> (right), we take as distance variables,

    r_13   =   |𝐱_1-𝐱_3|=|𝐱_2-𝐱_4|  ,
    
    
    r_14   =   |𝐱_1-𝐱_4|=|𝐱_2-𝐱_3|   ,

where, again, the second equalities are valid due to the particular characteristics of the geometrical distribution
of quarks and anti-quarks.




For this geometrical setup, one expects the ground state of the system when r_13≪ r_14 and r_14≪ r_13
to be dominated by the two possible independent two-meson states. For the computation of the static
potential we use the basis of operators shown in Fig. <ref> that are associated with the
two two-meson operators.






§ COMPUTING THE STATIC POTENTIAL 









For the computation of the static potential, including the groundstate and the first excited state, we rely on a basis of two operators 𝒪_i for each of the geometrical setups 
discussed in Sec. <ref>.  Defining the correlation matrix,

    M_ij   =   ⟨𝒪_i(0)^† 𝒪_j(t)⟩
       =   ∑_nc_in^*  c_jn  e^-V_nt ,

where ⟨⋯⟩ stands for vacuum expectation value, c_in=⟨ n|𝒪_i|0⟩
and |n ⟩ are the eigenstates of the Hamiltonian of the system,
the determination of the potential requires the knowledge of the solutions of the generalized  eigenvalue problem 

    M_ij(t)  a_j(t) = λ_k(t)   M_ij(t_0)   a_j(t) .

In our calculation, we assume that the creation of an excited state out of the vacuum occurs at t = 0. 
From the generalized eigenvalues λ_k, the energy levels of system V_k can be estimated from the plateaux on
the effective mass given by,

    M_eff(t)    =   logλ_k(t)/λ_k(t+1)
       =    V_k+𝒪(e^-(V_k+1-V_k)t)  .

In practice, the effective mass plateaux are identified fitting to a constant  both generalized eigenvalues.
In this way, one is able to compute both the static potential for the ground state and the first excited state of the system.




As described above, the basis of operators chosen to compute V depends on the geometry of the system 
and on the expected ground states. For the anti-parallel alignment, we use two meson-meson operators, while
for the parallel alignment a meson-meson operator and a diquark-antidiquark 
operator, i.e. a 3̅_123_34 colour configuration, are used to compute the correlation matrix.



In the case where the quarks are in the anti-parallel alignment  the operators used to compute the potential are,

    𝒪_13,24   =   1/3Q_1^iL_13^ijQ̅_3^j Q_2^kL_24^klQ̅_4^l ,
    
    𝒪_14,23   =   1/3Q_1^iL_14^ijQ̅_4^j Q_2^kL_23^klQ̅_3^l ,

where L are Wilson lines connecting the quark. Its representation in terms of closed Wilson loops is given in 
Fig. <ref>. The corresponding correlation matrix reads,

    M = 
    ([  W_13W_24 1/3W_1324; 1/3W_1423  W_14W_23 ])  ,

where W_i are normalized mesonic Wilson loops W=[U].



On the other hand, for the parallel alignment the two operators we consider are, 

    𝒪_YY  =   1/2√(3)Q_1^iQ_2^jϵ_i'j'kL_1a^ii'L_2a^jj'L_ab^kk'ϵ_k'l'm'L_b3^l'lL_b4^m'mQ̅_3^lQ̅_4^m , 
       𝒪_13,24     =  1/3Q_1^iL_13^ijQ̅_3^j Q_2^kL_24^klQ̅_4^l .

The closed Wilson loops associated to 𝒪_YY and 𝒪_13,24 are represented in Fig. <ref> and
the corresponding correlation matrix is given by,

    M =
    ([             W_YY 1/2√(3)W_YY,1324; 1/2√(3)W_1324,YY         W_13W_24 ])  .

















§ LATTICE SETUP 


From the static potential we aim to understand the transition between possible configurations of a
QQQ̅Q̅ system. Furthermore, we also want to glimpse any possible differences due to the quark dynamics.
Therefore, for the computation of V_k we consider two different simulations.  

Our quenched simulation uses an ensemble of 1199 configurations provided by the PtQCD collaboration <cit.>,  generated using the Wilson action in a
24^3 × 48 lattice for a value of β = 6.2. The quenched configurations were generated using GPU's and a combination of  Cabbibo-Marinari, pseudo-heatbath and over-relaxation algorithms, and computed in the GPU servers of the PtQCD collaboration.




Our full QCD simulation uses a Wilson fermion dynamical ensemble of 156 configurations generated in a 24^3 × 48 lattice and a β = 5.6. 
In the dynamical ensemble we take κ=0.15825 for hopping parameter, which corresponds to a pion mass of
m_π = 383 MeV. 
For Wilson fermions the deviations from continuum physics are of order 𝒪 (a)
in the lattice spacing and, therefore, one can expect relative large systematic errors. However, we expect that the static
potential as measured from the full QCD simulation away from the physical point to be more realistic when compared
to the quenched simulation.
The full QCD configuration generation has been performed in the Centaurus cluster <cit.> using the Chroma library <cit.>. The Hybrid Monte Carlo integrator scheme has been tuned using the methods described in 
<cit.>.


Then, with both the quenched and full QCD ensembles of configurations, we perform our correlation matrix computations at the PC cluster ANIMAL of the PtQCD collaboration.







The Wilson loops at large Euclidean time are decaying exponential functions of the static potential times the Euclidean time and, therefore, 
for large Euclidean times the Wilson loops are dominated by the statistical noise of the Monte Carlo. 
A reliable measurement of the static potential requires techniques which reduce the contribution of the noise to the correlation functions
used in the evaluation of V.

The quality of the measurement of the effective masses depends strongly on the overlap with the ground state of the system.
In order to improve the ground state overlap we applied  50 iterations of APE smearing <cit.> 
with w = 0.2 to the spatial links in both configuration ensembles. 
Furthermore, for the quenched ensemble, to further improve the signal to noise ratio, we used the extended multihit 
technique <cit.>. This procedure generalizes the multihit as described in <cit.>
by fixing the n^th neighbouring links instead of the first ones when performing the averages of the links. 
However, this technique has the inconvenient of changing the short distance behaviour of the correlators and, therefore,
one should not consider the points with r<r_min. In previous studies with the multihit, r_min=2 was sufficient, but in our study we consider r_min=4.
For the dynamical configurations the multihit technique can not be applied and, therefore, we resorted on
hypercubic blocking <cit.> with the parameters α_1 = 0.75, α_2 = 0.60 and
α_3 = 0.30 to improve the signal to noise ratio.











For the conversion into physical units we first evaluate Wilson loops to access the ground state meson static potential on a single
axis. In this calculation, we use a variational basis built using four different smearing levels to access the ground state meson
static potential. The lattice data for the static meson potential is then fitted to the Cornell potential functional form,

    V_M(r)=K-γ/r+σ r  .

The fits for different fitting ranges are reported in Tables <ref> and <ref> for
the quenched and the dynamical ensembles, respectively. 
The fits allows for the evaluation of the physical scale associated to the two ensembles through the Sommer method <cit.>.
Indeed, by demanding that,

    r_0^2dV_M/dr(r_0) = 1.65  ,

where r_0=0.5 fm, the lattice spacing a is measured and we present it in 
Tables <ref> and <ref> for various fitting ranges. The results show that a is fairly independent of
the fitting intervals and, in the following,  we take a ≃ 0.0681 fm for the quenched data ensemble and a ≃ 0.0775 fm
for the dynamical data set. 
Our QCD lattice spacing is essentially similar to the one obtained with different techniques.
It follows that the lattice volumes used in the simulation are
( 1.63fm)^3 × 3.27fm for the quenched case and 
( 1.86fm)^3 × 3.72fm for the dynamical simulation.
For completeness, in  Fig. <ref> we show the ground state meson potentials for the two ensembles in physical units.








§ RESULTS FOR THE STATIC QQQ̅Q̅ POTENTIAL 


In this section, we report on the results for the static potential with the two different geometries mentioned  in Sec. <ref>, and we apply fits with ansatze bases in the string flip-flop potential and in the Casimir scaling.












In Fig. <ref>, as an example, we show effective mass plots for the pure gauge simulation (left), full QCD simulation (right) and for the ground
state (top) and first excited state (bottom) for a QQ Q̅Q̅ system in the antiparallel geometry. 
The red curves are the results of fitting the lattice data to measure the static potential. See the appendix for further details on the numerics. 
We consider the maximum number of points aligned in a horizontal line with acceptable χ^2 / d. o. f.. 
Because the noise reduction technique in the quenched simulation rejects the cases with source distances smaller than 4a, we end up by accepting a few more results in the full QCD case than in the quenched case.




 §.§ The anti-parallel alignment








We start by analyzing the simpler case of the anti-parallel geometry, where the meson-meson systems are expected to have lower energies than the tetraquark system.  Our results are plotted in Figs. <ref> and <ref>. Clearly there are two different trends for r_13 < r_14 and for r_13>r_14 and a transition, with mixing, at the point r_13=r_24. Moreover we compare in detail our results with different ansatze.


From the string flip-flop paradigm of Fig. <ref> we expect  the ground state of the system to be that of a two meson system 
when the distance between a quark and an anti-quark, i. e. r_13 or r_14, is much smaller than the quark-quark distance, i.e. r_12.
Then, for sufficiently small r_13 and/or r_14 the potential of the ground state of the QQQ̅Q̅ should reproduce the string flip-flop potential,

    V_0   ≃    V_ff = min[ V_MM,V_MM'] ,
 
where the two different meson -meson potentials are

    V_MM   =    2 V_M(r_13)  ,
    
    
    V_MM'   =    2  V_M(r_14)  ,

and V_M is the ground state potential of a meson in Eq. (<ref>).
Previous lattice simulations <cit.> confirm that V_0 is compatible with
such a result. Deviations from Eq. (<ref>) are expected at intermediate distances together with a smooth transition
from one picture to the other, i.e. from the two meson state with valence content Q_1Q̅_3 and
Q_2 Q̅_4 to the two meson with valence content Q_1Q̅_4 and
Q_2 Q̅_3.



On the other hand, for the excited state, we have two possible scenarios. From the string-flip-flop, we would again expect, when the distance between quark and anti-quark  a quark and an anti-quark, i. e. r_13 or r_14, is much smaller than the quark-quark distance, i.e. r_12, the system to be that of the next two meson system,

    V_1?≃max[ V_MM,V_MM']  .

However, given that the colour wavefunctions of the two mesonic states are not orthogonal, see Eq. (<ref>), and Section II,
possibly the excited state is not another mesonic state and, but instead is an octet state, 

    V_1?≃max[ V_88,V_88') ]  ,

where we estimate the colour octet potential assuming Casimir scaling, i.e. using the decomposition  in Eq. (<ref>) and the values reported on 
Tab. <ref>,

    V_88     =    1/2V_M(√(r_13^2+r_14^2))+7/4V_M(r_14)-1/4V_M(r_13)  
     ,
    
    
    V_88'    =    1/2V_M(√(r_13^2+r_14^2))+7/4V_M(r_13)-1/4V_M(r_14)  ,




















Thus we have two different simple anzatse to interpret our results. 
The ground state potential V_0 and the first excited state potential V_1 for the quenched and dynamical ensembles are 
reported in Figs. <ref> and <ref>, respectively, together with V_MM,
V_MM' and the octet potentials V_88, V_88'.

As the figures show, the ground state static potential V_0 as a function of r_14 is compatible with two two meson potentials
for small and large values of r_14. 
Indeed, for all r_13, at small values of r_14 the static potential is compatible with V_MM, while for large r_14 
V_0 becomes compatible with V_MM'. 
We show in Table <ref> good fits with the meson-meson potentials. 
This type of behaviour is well described by the string flip-flop potential V_ff.






In the transition region r_13∼ r_14 where also V_MM∼ V_MM', deviations of V_0 from V_MM or V_MM' can be seen. 
The difference between the ground state potential and the sum of the two meson potentials 
in physical units is detailed in Fig. <ref>, and in particular the transition point 
r_12 = r_13 is analysed in Fig. <ref>.
The results for the quenched simulation are well described assuming an off diagonal term Δ in the correlation matrix, leading to the functional form,

    V_0(r_13,r_14)   =   V_MM+V_MM'/2
          
    -√((V_MM-V_MM'/2)^2+Δ^2) ,

where we may have either,

    Δ(r_1,r_2)=Δ_0e^-λ(r_1+r_2)/1+c(r_1-r_2)^2

or,

    Δ(r_1,r_2)=Δ_0/1+c(r_1-r_2)^2+d(r_1+r_2)^2 .

Eq. (<ref>) interpolates between the two potentials in flip-flop picture of a meson-meson.





The fits for the functional forms in Eqs. (<ref>) and (<ref>) are reported in Tables <ref> and <ref>. 
In order to quantify the deviation from the two limits where the system behaves as a two meson system, 
we refer that the fits give a Δ(0.5 ,0.5 )≃60, a number to be compared
with typical values for the meson potential which are of the order of GeV (see Fig. <ref>).
This results shows that the corrections due to Δ to the flip-flop picture are small when the quarks and anti-quarks are in an
anti-parallel geometry. 

























The full QCD simulation shows similar results to the quenched QCD simulation. 
However, the results for V_0 for the full QCD configurations are not described by the same
type of functional form  given in Eq. (<ref>)
which reproduces the flip-flop potential at large distances. We found no window where the fits are stable and, therefore, conclude that the dynamical V_0 is not reproduced by 
Eq. (<ref>) with the deviations parametrised by either Eq. (<ref>) or Eq. (<ref>). 


In what concerns the excited state potential V_1 there are clearly two different regimes for r_13 very different from r_14, but we are not able to find an
analytic form compatible with the lattice data, neither for the quenched simulations nor for the full QCD simulations. 
In both  Figs. <ref> 
and <ref>, it is clear the static potential V_1 lies between the functional forms of Eq. (<ref>) and Eq. (<ref>).
There are subtle differences between Fig. <ref> 
and <ref>. In general, the full QCD case is closer to the octet expression of Eq. (<ref>) than the quenched QCD case.

A fortiori, we are not able as well to find a good ansatze to fit V_1 in the transition region. For a detailed view of the differences for the quenched simulation in this region, see Fig. <ref>.


This observed behaviour for V_1 can be understood in terms of adjoint strings.
When the quark-anti-quark inside the octets are close to each other, they can be seen externally
as a gluon. Therefore, we have a single adjoint string with a tension of σ_A=9/4σ. 
On the other hand, when the quark and the anti-quark are pulled apart, the adjoint string tends to split into 
two fundamental strings, with a total string tension of 2σ. The splitting of the adjoint string, gives a 
repulsive interaction between the quark-anti-quark pairs that form octets in the excited state. This is qualitatively
consistent with the behaviour predicted by Casimir scaling, where the potential for a quark and an anti-quark in an octet corresponds
to a repulsive interaction.





  §.§.§ Mixing angle





For the anti-parallel geometry and for the ground state potential the lattice results show that the tetraquark is essentially a two meson state. Therefore,
one can write the most general ket describing the ground state |u_0⟩ of a QQQ̅Q̅ system as a linear combination of the available
colourless states

    |u_0⟩    =    cosθ  |6_126̅_34⟩+
                  sinθ  |3̅_123_34⟩
        =    √(3/4){( cosθ/√(2) + sinθ)  |1_131_24⟩
               + ( cosθ/√(2) - sinθ)  |1_141_23⟩} .

For a pure two-meson state, the mixing angle is either θ = θ_0,
for |1_131_24⟩, or θ = -θ_0, for  |1_141_23⟩,
with θ_0 = tan^-1 ( 1 / √(2) ). For the general case, the angle θ can be estimated using the generalized
eigenvectors obtained solving Eq. (<ref>) with the following operators,

    𝒪_S   =   √(3/8)(𝒪_13,24+𝒪_14,23)   ,
    
    𝒪_A   =   √(3/4)(𝒪_13,24-𝒪_14,23)  .

The results for θ for the quenched simulation can be seen in Fig. <ref>.
From the lattice data one can estimate a typical length, or broadness, associated to the transition between the two two-meson states. In the region when |r_13-r_14| ≲ d_trans, the transition occurs and the groundstate is a mixing of the MM and MM' states.  We estimate the typical transition length from,

    d^-1_trans∼dθ(r_13,r_14)/dr_14|_r_14=r_13 .

For the quenched data, see Fig. <ref>, the derivative stays within 0.36/a and 0.42/a
and, therefore, d_trans∼ 0.16 - 0.19 fm. 
For the dynamical simulation, see Fig. <ref>, the typical transition length is essentially the same and we find
d_trans∼ 0.16 -  0.20 fm.










The lattice data for the mixing angle gives a vanishing angle for r_13=r_14.
This means that the ground state for the anti-parallel alignment is given only by |6_126̅_34⟩
and has no |3̅_123_34⟩ component.

The results reported in Figs. <ref> and <ref> show that, in general, a QQQ̅Q̅ system
is in a mixture of two possible colour meson states and it approaches meson states as the distance between the pairs of quark-anti-quark is much
smaller than the distance between quarks or anti-quarks.









 §.§ The parallel alignment


For this particular geometry, the static potential was investigated with lattice methods in <cit.>. 
For the ground state and in the limit where r_12≪ r_13, the authors found that the lattice data is compatible with the double-Y (or butterfly) potential,



    V_YY   =    2K - γ(
        1/2r_12 + 1/2r_34+1/4r_13
          +1/4r_24+1/4r_14+1/4r_23)+σ L_min ,

where γ and K are the estimates of the static meson potential and σ is the fundamental string tension.
For the geometry described on the right hand side of Fig. <ref> and for r_13 > r_12 / √(3) the butterfly potential simplifies into,

    V_YY   =   2K-γ(1/r_12+1/2r_13
           +1/2√(r_12^2+r_13^2))+σ(√(3)r_12+r_13)  ,











Moreover, from the expression for the Casimir scaling potential given in (<ref>) and using the results reported on Tab. <ref>
it is possible to define various types of potentials to be compared with the static potential computed from the lattice simulations.

The potential associated to the state where the quarks and anti-quarks are in triplet states leads to the so-called triplet-antitriplet or diquark-antidiquark potential,

    V_33=∑_i<j⟨3̅_123_34|C_ij|3̅_123_34⟩ V_M(r_ij)  ,


or in a form similar to (<ref>),


    V_33   =   2K-γ(1/r_12+1/2r_13+1/2√(r_12^2+r_13^2))
    
           +σ(r_12+1/2r_13+1/2√(r_12^2+r_13^2)) .


Similarly, the anti-sextet-sextet potential is given by

    V_66   =   ∑_i<j⟨6_126̅_34|C_ij|6_126̅_34⟩ V_M(r_ij)
    
       =   5/4V_M(r_13)+5/4V_M(√(r_12^2+r_13^2))-1/2V_M(r_12)  ,

and the octet-octet potential reads

    V_88=1/2V_M(r_12)-1/4V_M(r_13)+7/4V(√(r_12^2+r_13^2))  .






The lattice estimates for the ground state and first excited (whenever possible) potentials can been in Figs. <ref> 
and <ref> for the quenched and for the dynamical simulation, respectively. 
The data shows that for large quark-anti-quark distances, i.e. for large r_13, the static potentials are compatible with a linearly 
rising function of r_13.  
This result can be viewed has an indication that the fermions on a tetraquark system are confined particles.





For both the pure gauge and dynamical simulations and for small quark-anti-quark distances, i.e. for small r_13,
and up to r_13≤  r_12 the ground state potential reproduces that of a two meson state V_MM.
In this sense, one can claim that for sufficiently small quark-anti-quark distances the ground state of a QQQ̅Q̅ system 
is a two meson state.
For the excited potential, the pure gauge results are among the double-Y potential (<ref>)
and the octet-octet potential (<ref>). However, for the dynamical results, the static
potential seems to be closer to V_88 at smaller and large r_13 and closer to V_YY as r_13 approaches r_12.


On the other hand for sufficiently large r_13, the ground state potential is essentially that of a diquark-antidiquark system
V_33 and the system enters its tetraquark phase. 
Indeed, the ground potential is given by 2 V_M for quark-anti-quarks distances up to r_13 = r_12 and
is just above V_YY for distances r_13≥ r_12 + 1 in lattice units. These results suggests that, for this geometrical setup,
the transition of a two meson state towards a tetraquark state occurs at r_13∼ r_12 + 1 (in lattice units).

In what concerns the dependence of V_0 on r_12, the lattice data suggests that the potential increases with the 
quark-quark distance and favours a V_0 ∼ V_YY for sufficiently large r_12 as was also observed in <cit.>. 


For the quark models with four-body tetraquark potentials, in particular the string flip-flop potential illustrated in Fig. <ref> it is very important to quantify the deviation of V_0 from the V_YY ansatz; and we have studied several ansatze for this difference. Clearly V_0 is more attractive than the tetraquark potential V_YY of Eq. (<ref>) reported by previous authors, and this favours the existence of tetraquarks.

Adding a negative constant (attractive) to the double-Y potential is not sufficient for a good fit of the lattice data for any of the sets of configurations.
Adding a correction to the double-Y potential which is linear in the
quark-quark distance,

    V_YY^B=V_YY(r_12,r_34) + δ K + δσ_12  r_12 ,

describes quite well the dynamical simulation data and a fit gives δ K = -0.12(3)  √(σ), δσ_12=-0.34(5)  σ,
where σ is the fundamental string tension, for a χ^2 / d.o.f =0.46 (see the tables on the appendix for details on the fits). 
The dynamical data for the deviations from V_YY are also compatible with a Coulomb like correction

    V_YY^C = V_YY(r_12,r_34) + δ K + δγ_12/r_12 ,

for a δ K=-0.67(4)  √(σ), δγ_12=0.22(3) 
with a χ^2/d.o.f. = 0.62 (see appendix for details).
Such a functional form is not compatible with the lattice data for the pure gauge case. A possible explanation could come from
the difference in the statistics of both ensembles. Recall that the number of configurations for the pure gauge
ensemble is about ten times larger than for the dynamical simulation and, therefore, the associated statistical errors are much smaller.


In what concerns the first excited potential V_1, the data for the pure gauge and for the dynamical fermion simulation follows slightly 
different patterns. In the quenched simulation and for r_13 < r_12, the potential is close to V_T and the behaviour for
larger values of r_13 does not reproduces any of the potentials considered here. On the other hand, in the dynamical simulation
V_1 for small and large values of r_13 is just below the data for anti-sextet-sextet potential

    V_66=∑_i<j⟨6_126̅_34|C_ij|6_126̅_34⟩   V_M(r_ij)

which, for this geometry, is given by

    V_66=5/4V_M(r_13)+5/4V_M(√(r_12^2+r_13^2))-1/2V_M(r_12)

and, at intermediate distances where r_13∼ r_12, is compatible with the octet-octet potential, 

    V_88=1/2V_M(r_12)-1/4V_M(r_13)+7/4V(√(r_12^2+r_13^2))  .

Further, at very small distances the potential seems to flatten for full QCD and the data also suggests a flattening or a small repulsive core. 
Note, for the quenched simulation smaller distances than 4 are not accessible, this short distance effect is not visible. 






§ SUMMARY AND DISCUSSION 


In this work the static potential for a Q Q Q̅Q̅ system was investigated using both quenched and Wilson Fermion full QCD simulations for two
different geometric setups. The two geometries are designed to investigate sectors where dominantly meson-meson or tetraquark static potentials are expected. 


The simulations show that whenever one distance is much larger than the other, the ground state potential and the first excited state potential are compatible
with a linearly rising function of the distance between constituents, suggesting that quarks and anti-quarks are confined particles. For the distances studied, the quenched and full QCD results are qualitatively similar, and their subtle differences only become clearer when we compare the lattice data with anzatse inspired in the string flip-flop potential and in Casimir scaling. 


For the anti-parallel geometry setup, the groundstate potential V_0 is approximately described by a sum of two two meson potentials, i.e. it is compatible with the string flip-flop type of potential. We take this result as an indication that the Q Q Q̅Q̅ wave function is given by a superposition of two meson
states and we compute the mixing angle, as a function of the quark-anti-quark distances, which caracterize such a quantum state. 
The mixing angle shows that the tetraquark system undergoes a transition from one of the meson states to the other configuration as the
quark-antiquark distance increases, and the broadness of this transition has a typical length scale of 0.16 - 0.20 fm.
Moreover, for the quenched simulation, we found an analytical expression which describes well the lattice groundstate.
The analytical expression is essentially a flip-flop type of potential with corrections, parametrized by Δ (r_1, r_2), which are typically 
≲ 10% than the sum of two two mesons potentials. 

In what concerns the first excited potential V_1 in the anti-parallel geometry, the results show that for small enough
quark-anti-quark distances the potential is just below one of the possible octet-octet potentials and approaches a two meson potential from above
from large quark-anti-quark distances. This results for the excited potential can be interpreted in terms of and excited state including a combination of meson-meson and octet-octet states.

For the parallel geometry setup, the groundstate potential V_0 is compatible with  a diquark-antidiquark potential for large quark-antiquark distances and a sum of two meson potentials for small separations. Moreover, the lattice data for the full QCD simulation is 
compatible with a butterfly type of potential with corrections that we are able to parametrize. For the quenched simulation we found no analytical
expressions that are able to describe the data, but the trend is the same. 

The interpretation of the first excited potential V_1 for the parallet geometry, in terms of possible colour configurations is not as compliant with models as in the anti-parallel geometry. It seems that V_1 for the full QCD simulation
is just below the octet-octet from small quark-anti-quark distances and approaches again the octet-octet potential for at large distances. For the quenched simulation, the interpretation of V_1 in terms of colour components is not so clear, as the lattice data seems to point for a combination of different colour potentials.


Importantly for quark models with four-body tetraquark potentials, in particular for the string flip-flop potential illustrated in Fig. <ref>, we obtain a groundstate potential V_0 more attractive, by a difference of -300 to -500 MeV, than the butterfly potential reported by previous authors
<cit.>, and this favours the existence of tetraquarks.

As an outlook, it would be interesting to measure the static Q Q Q̅Q̅ potentials for larger distances and for different geometries. We leave this for future studies.



The authors are extremely grateful to Nuno Cardoso <cit.> for generating the ensemble of quenched configurations utilized in this work. The authors also acknowledge both the use of CPU and GPU servers of the collaboration PtQCD <cit.>, supported by NVIDIA, CFTP and FCT grant UID/FIS/00777/2013, and the Laboratory for Advanced Computing at University of Coimbra  <cit.> for providing HPC resources that have contributed to the research results reported within this paper. 
M. C. is supported by FCT under the contract SFRH/BPD/73140/2010.
P. J. S. acknowledges support
by FCT under contracts SFRH/BPD/40998/2007 and SFRH/BPD/109971/2015.



ieeetr











§ TABLES OF RESULTS 





















