Disentangling different structures in heavy-light four-quark states
Paul C. Wallbott Gernot Eichmann Christian S. Fischer
hep-ph hep-lat









[e-mail: ]paul.wallbott@physik.uni-giessen.de
Institut für Theoretische Physik, Justus-Liebig Universität Gießen, 35392 Gießen, Germany


[e-mail: ]gernot.eichmann@tecnico.ulisboa.pt
LIP Lisboa, Av. Prof. Gama Pinto 2, 1649-003 Lisboa, Portugal
Departamento de Física, Instituto Superior Técnico, 1049-001 Lisboa, Portugal


[e-mail: ]christian.fischer@physik.uni-giessen.de
Institut für Theoretische Physik, Justus-Liebig Universität Gießen, 35392 Gießen, Germany
Helmholtz Research Academy Hesse for FAIR (HFHF), Campus Gießen, 35392 Gießen, Germany




Models proposed to explain recently discovered heavy-light four-quark states already assume certain internal
structures, i.e. the (anti)quark constituents are grouped into diquark/antidiquark clusters, heavy-meson/light-meson
clusters (hadrocharmonium) or heavy-light meson molecules. We propose and use an approach to
four-quark states based on Dyson-Schwinger and Bethe-Salpeter equations that has the potential to discriminate
between these models. We study the masses of heavy-light cqq̅c̅ and ccq̅q̅ four-quark states 
with q=u,d,s and quantum numbers I(J^PC)=0(1^++),1(1^+-),0(0^++) and 1(0^+),0(1^+),1(1^+). 
We identify the dominant components of the ground states with these quantum numbers and suggest candidates for corresponding
experimental states. Most notably, we find strong heavy-light meson-meson and negligible diquark-antidiquark
components in all cqq̅c̅ states, whereas for ccq̅q̅ states diquarks are present. 
A potential caveat in the I=0 channels is the necessary but costly inclusion of cc̅ components 
which is relegated to future work. 


Disentangling different structures in heavy-light four-quark states
    Christian S. Fischer
    December 30, 2023
===================================================================



In the past two decades a number of highly interesting states have been identified in the charmonium and bottomonium
energy regions that cannot be accommodated for in the conventional quark model for mesons
made of a quark and an antiquark. Since the quark model is otherwise extremely successful in predicting
spectra of heavy QQ̅ states (with Q=c,b), these exceptional states are considered to be exotic hadrons.
Some of them carry electromagnetic charge and thus may be naturally explained as four-quark states QQ̅qq̅ (q=u,d,s)
with a light charged quark-antiquark pair in addition to the overall neutral QQ̅ component.
Thus, four-quark states are considered as promising candidates to explain the properties of these exotic hadrons, see e.g. <cit.> for reviews.

There is, however, no agreement on the internal structure of these four-quark states.
Model approaches usually assume some kind of internal clustering from the start. One possibility, the
hadroquarkonium picture <cit.>, suggests a heavy quark and antiquark grouped together
in a tight core surrounded by the light qq̅ pair. This is motivated by the experimental observation
of final states with a specific charmonium state and light hadrons. The
second possibility is the clustering of constituents in diquark-antidiquark ()
components which interact via colored forces, see e.g. <cit.> for a review. A third
possibility, especially relevant for states close to open-charm thresholds, is the
meson-molecule picture of arrangements into pairs of D^(*)D^(*)  mesons
that interact with each other by short- and/or long-range forces <cit.>.

It is important to note that these possibilities are not mutually exclusive: In general, every experimental
state may be a superposition of components with a different structure and the `leading' component may be different
on a case-by-case basis. It is therefore important to develop theoretical approaches to QCD that can
deal with all these possibilities. Lattice QCD is one such approach and has made
interesting progress so far, see
<cit.> and references therein.
Nevertheless most simulations have been performed at an exploratory level using light quarks with unphysical
large masses. Functional methods, on the other hand, have
been restricted to four-quark states with equal masses <cit.> or specific
quantum numbers <cit.>.

In this work we present a generalization of the functional approach to four-quark states that
has the potential to systematically address and compare heavy-light states in different flavor combinations
and with different J^PC quantum numbers. Based on a well-studied and understood truncation of the
underlying quark-gluon interaction, we work with an approximated version of the four-body Faddeev-Yakubovsky
equations that takes into account the two-body correlations that lead to the internal clustering described
above. We apply the resulting formalism to the experimentally interesting
cqq̅c̅ hidden-charm states with quantum numbers
J^PC=0(1^++), 1(1^+-) and 0(0^++), which are carried by the
X(3872) <cit.>, the neutral Z(3900) <cit.>
and (likely) the X(3915) <cit.>, respectively. Furthermore, we discuss four-quark states
with open charm ccq̅q̅ in the channels 1(0^+),0(1^+) and 1(1^+).
Currently there are no experimental candidates for these states but searches are underway.
Corresponding states in the heavier bottom-quark region received a lot of attention in recent years
since they are promising candidates for deeply bound and narrow states, see e.g.
<cit.>. It is certainly interesting
to see whether this is still the case for the experimentally more easily accessible open-charm states.




Four-body equation.
The homogeneous Bethe-Salpeter equation (BSE) shown in Fig. <ref> has the form
  
    Γ = K G_0  Γ  ,

  where Γ is the BS amplitude, K is the four-quark interaction kernel that contains all possible  three- and four-body interactions,
  and G_0 is the product of four dressed (anti)quark propagators; see <cit.> for details.
  Each multiplication represents an integration over all loop momenta.
  Eq. (<ref>) holds at a given pole position of the offshell qqq̅q̅ scattering matrix T, which satisfies the scattering equation T = K + K G_0 T.
  Poles on the real axis of the total squared momentum P^2 correspond to bound states,
  whereas resonances appear as poles in the complex plane on higher Riemann sheets.


 




In this work we focus entirely on the two-body correlations in K since these generate the internal two-body clusters discussed above.
This leads to

    KG_0 = ∑_aa' K_aa' ,    K_aa' = K_a + K_a' - K_a K_a'

where a, a' stand for qq, q̅q̅ or qq̅ pairs and aa' is either (12)(34), (13)(24) or (14)(23).
The subtraction is necessary to avoid overcounting <cit.>.
Irreducible three- and four-body interactions are not (yet) taken into account for three reasons: first, this would complicate an already
tremendous numerical task further beyond the resources currently available to us; second, a similar strategy has been employed with
great success in the baryon sector, where strong two-body correlations naturally lead to a diquark-quark picture, which in turn leads
to a spectrum in one-to-one agreement with experiment <cit.>; third and most important, the pictures
of internal structures that we like to discriminate (i.e. diquark/antidiquark vs. hadro-charmonium vs. meson molecule) all rely on
strong two-body clustering. Thus for the purpose of this work it is indeed sufficient to focus on two-body interactions. Nevertheless,
of course, the effects of irreducible three- and four-body forces have to be explored in future work.


For the two-body kernels we employ the same rainbow-ladder interaction that is used in the Dyson-Schwinger equation (DSE)
for the quark propagator. This truncation has recently been reviewed in <cit.>, where the DSE for the quark
propagator is discussed around Eq. (3.18) and the effective interaction
in Eqs. (3.95–3.96). We use Λ=0.72 GeV for the scale parameter, adjusted to reproduce the pion decay
constant f_π, and η=1.8 ± 0.2. Together with the current-quark masses, these are the only input parameters in all equations.
The construction satisfies chiral constraints such as the Gell-Mann-Oakes-Renner relation,
ensures the (pseudo-) Goldstone-boson nature of the pion and has been extensively applied to meson and baryon phenomenology.
As discussed in <cit.>, the truncation is well-known
to reliably reproduce  many properties of pseudoscalar and vector mesons (and, correspondingly, scalar 
and axialvector diquarks). Since we focus on two-body clusters inside tetraquarks in these channels only,
we may expect qualitatively reasonable results.

The quantitative reliability of the approximation of the two-body kernel may be judged from the results for meson masses in Table <ref>.
We work in the isospin symmetric limit where m_D^+ = m_D^- = m_D^0. The u/d current-quark mass is fixed by m_π, the strange quark
mass is chosen such that the sum m_D_s+m_D_s^* equals the sum of the experimental values <cit.> and analogously
for the charm quark mass in m_D+m_D^*. The deviations between the theoretical and experimental meson masses are then
below 7% in all cases.


Four-quark amplitude.
  The main challenge in solving Eq. (<ref>) for given J^PC is the structure of the BS amplitude.
  Its general decomposition can be written as
  
    Γ^(μ)_αβγδ(p_1 … p_4) = ∑_i f_i(…) τ_i^(μ)(p_1 … p_4)_αβγδ ,

  where the Lorentz-invariant dressing functions f_i(…) depend on the ten Lorentz invariant momentum variables
  that can be constructed from four independent momenta. The tensors τ_i are the direct products of Dirac, color and flavor parts.
  A J=0 state has 256 linearly independent Dirac tensors and a J=1 state  768,
  which are collected in Ref. <cit.> and the Appendix of <cit.>.
  The color part of the amplitude consists of two independent
  color-singlet tensors and the flavor wave functions depend on the particular system, cf. Appendix <ref> for details.

  To extract physical content from the BS amplitude Γ^(μ)(p_1 … p_n), we observe that the amplitude develops
  internal two-body clusters, which for heavy-light systems occur in the three different channels corresponding to
  hadroquarkonia, heavy-light meson-meson components and dq-d̅q̅ clusters.
  These clusters may go on-shell provided that the sum of their masses is smaller than the mass of the
  four-body state. If this occurs in color-singlet channels, the four-quark state becomes a
  resonance in the two-body hadronic system of the corresponding clusters. But even if the masses of the
  two-body clusters are large enough such that the probed momenta only come close to the corresponding
  singularities, this will influence the four-body system. Thus the guiding idea is to represent Γ^(μ)(p_1 … p_n)
  in terms of these two-body clusters.

  Since we are interested in specific quantum numbers with experimental candidates for four-quark
  states, we draw on existing information on the decay channels of these states and construct our
  representation along the content displayed in Table <ref>. For example, in the heavy-light meson
  sector we took into account combinations of the I(J^P)=1/2(0^-) multiplets (D,D^*) and their anti-particles (D̅,D̅^*),
  omitting the heavy combination D^*D^*. Note that since we work in the isospin-symmetric
  limit the charged and neutral states are mass-degenerate and both taken into account.

  

  The construction of the τ_i^(μ)(p_1 … p_4) for 
  the configurations qqqq, cqqc and ccqq is detailed in App. <ref>.
  We construct the Dirac parts according to the dominant two-body clusters 
  and combine them with appropriate color and flavor wave functions such that the charge-conjugation and Pauli exchange symmetries are respected. 
  As a result, one populates a physically motivated subset of all possible basis elements. 
  Whereas in this work we restrict ourselves to the combinations displayed in Table <ref>,
  in principle one could construct a complete basis for the four-body amplitude with entangled Dirac, color and flavor tensors
  including all possible meson and diquark channels (i.e., also those with higher total angular momentum).

With this setup, we are in a position to solve the BSEs as an eigenvalue problem with structure 
λ Γ = K G_0  Γ and general eigenvalue λ. All elements of this equation depend
explicitly on the total momentum P. By varying P^2 such that λ=1 one finds the mass of the bound 
state/resonance in the four-body system via P^2=-M^2. 
One problem that appears in this search is the potential appearance of singularities
in the plane of the complex total momentum due to the internal meson and diquark correlations. Typically, this does not happen for large masses of the
lighter quark pair, where the resulting bound state is in general (well) below the meson-meson thresholds, cf. Figs. <ref> and <ref>.
However, this situation changes
for smaller masses and in particular close to the physical point of light quark masses, where the mass of the bound state/resonance comes
close to or is even larger than the meson-meson threshold. For these cases we determine the eigenvalue curve λ(M^2)
in the singularity free region and extrapolate it further into the time-like momentum domain using rational functions.
This procedure can only pick up the real part of potentially complex masses, i.e. it is not possible to extract decay widths. 
To this end one would need to use the much more involved approach described in <cit.>, 
which is not (yet) at our disposal for heavy-light four-quark states. For many further details and first results in the 
light quark sector we refer the reader to Ref. <cit.>. 


Hidden-charm states.
We first discuss our results for the hidden-charm cqqc four-quark states in the I(J^PC)=0(0^++) and 1(1^+-) channels.
In this context we wish to emphasise that the results in the I=0-channel must be seen as preliminary since the necessary but costly inclusion 
of cc̅ components in our approach is relegated to future work (see however <cit.> for 
corresponding results in the light quark sector). 

The results for the 0(1^++) heavy-light state can be found in <cit.>, where we also described our procedure
to estimate (part of) the error of the calculation. In that case we found a dominant
heavy-light meson component, whereas the hadrocharmonium component is rather weak and the diquark component
has almost no effect at all.




A similar pattern arises in the scalar 0(0^++) case displayed in the top panel of Fig. <ref>.
We show the mass evolution of the four-quark
state when the mass of the cc̅ pair is fixed and the 
mass of the other qq̅ pair is varied from the charm mass (rightmost vertical dashed line) to the strange and
light quark masses (other two vertical lines). We compare calculations with hadrocharmonium content
only (J/Ψ ω; squares), heavy-light meson components only (DD; circles) and a combination of the two (triangles).
The full calculation including also the diquark/antidiquark component (DD+J/Ψ ω+SS) is marked with crosses.
Not contained in the figure is our result for using diquark components only (SS). For the scalar channel and in fact for
all channels that we studied we obtain very large masses with diquarks only. These are typically
of the order of the diquark/antidiquark thresholds around 5 GeV (cf. Table <ref>) and therefore unphysical.
Comparing the results for 0(0^++) shown in Fig. <ref> we find that for large masses, the  component
has almost no effect on the results, whereas the hadrocharmonium component gives only mild corrections
to the leading heavy-light meson components.
The diquark corrections become somewhat more prominent for small masses, however without changing the general picture.
Due to the sizeable error bar we cannot discriminate between a bound state and a resonance.
The masses quoted in Table <ref> are obtained from
a linear fit to the mass evolution at larger quark masses;
only for 0(0^++) the direct calculation yields a mass below any threshold.
 









The general picture changes somewhat for the 1(1^+-) state shown at the bottom of Fig. <ref>.
Once again we find negligible diquark components and a strong heavy-light meson component, but also non-negligible
contributions from the hadrocharmonium component. It is interesting to compare this behavior to the
1^++ channel discussed in <cit.>. From Table <ref> we observe that
the only difference between the two states are the hadrocharmonium components.
The J/ψ π component in the 1^+- channel is significantly lighter
than the J/ψ ω component in the 1^++ channel, which lifts the degeneracy between the two states and
leads to a lighter mass of the 1^+- (which is, however, opposite to the current experimental situation).
We have also tested further components, which can contribute to the axialvector 1(1^+-) state. We
found the D^∗ D^∗ components to be negligible; however, the η_c ρ component is sizeable
and enhances the mass splitting once included. Further studies in this direction are necessary.

The resulting masses are collected in Table <ref>. It is interesting to compare
our results with expectations from the literature. In Ref. <cit.> heavy-quark symmetry has been
used to predict patterns for molecular states. This led to the identification of the X(3872) with a molecular
state in the 1^++ channel and the neutral Z_c(3900) with a molecular state in the 1^+- channel. Our
results agree with this identification: for the X(3872) we expect an almost pure heavy-light meson state which
is then natural to expect to sit very close to the DD̅^∗ threshold. For the Z_c(3900), however, we find
non-negligible corrections from other components, which may shift the physical state away from the
threshold. For the scalar channel no predictions
have been made in <cit.>, see however <cit.>
for a detailed discussion. The lightest scalar molecule, if it exists, would be expected
at the DD̅ threshold. This is indeed the case for our scalar state, which sits in the region of the
threshold of our D mesons (cf. Table <ref>). Thus the mass pattern emerging in Table <ref>
is in line with our observation of heavy-light meson dominance in all cqq̅c̅ states studied so far.

In Table <ref> we also list the masses of the charm-strange css̅c̅ states
extracted from the mass evolution.
For I=0 these correspond to observable states, whereas for
I=1 they are unphysical. 
In  the 1^++ channel there is an experimental candidate, namely the X(4140) with a mass only
slightly above the upper range of our error bar. Provided this identification holds, we predict a strong
heavy-light meson component of the X(4140), even though it is not overly close to the D_sD̅_s^∗ threshold.
In addition, we find a corresponding state in the 0^++ channel, although the large error bars in this
case make a prediction of its mass rather imprecise.





Open-charm states.
Our results for the open-charm states with flavour content ccqq are shown in Fig. <ref>.
The basis construction in the open-charm case is significantly different from hidden charm
since charge-conjugation symmetry is replaced by Pauli antisymmetry.
The two heavy-light meson combinations (cq̅)(cq̅) are identical,
whereas the  component (cc)(q̅q̅) with a heavy diquark and a light antidiquark
inherits the role of the hadrocharmonium component. For I=0,
 the light qq antidiquark must be scalar (S) due to symmetry, whereas for I=1 it is
axialvector (A). The heavy cc pair is always an axialvector diquark.



This change of internal
dynamics is also reflected in the results. The heavy-light meson component alone produces a state that is well
below the DD^* threshold and moves only slightly above threshold for decreasing quark masses.
Whereas the AS diquark contribution in the I=0 case is negligible compared to the
DD^* contribution, the AA diquark component for I=1 has a significant impact in pushing the mass evolution up
 above threshold.
 In the graph for 0(1^+) we also included the D^∗D̅^∗ component explicitly, which is negligible as in all other channels.
 The resulting mass hierarchy between the isosinglet and isotriplet states
is as expected from heavy quark symmetry <cit.>. Our extrapolated values for the masses are in
the ballpark expected from other approaches, see e.g. <cit.> and Refs.
therein.

In the open-charm case, the ccss states (with I=0) must be read off from the ccqq curves with I=1
since those have the same wave-function components cf. Appendix <ref>.
As a consequence, several slots in Table <ref> are empty because they do not support physical states.
Moreover, there is a large gap between the light and strange state in the 0(1^+) channel but
it comes again with a sizeable mass uncertainty.



Conclusions.
In this work we have studied and compared the masses of heavy-light four-quark states in the charm energy region.
We developed a dynamical framework that takes into account all possible combinations of internal two-body
clusters and is therefore able to decide dynamically whether pictures from effective field theory and models
(meson molecule, hadrocharmonium, diquark-antidiquark) are realized. For hidden charm, in all cases considered
we do not find a sizeable diquark-antidiquark component. Instead, the heavy-light meson component
is favoured, with channel-dependent negligible (for 0(1^++)) or small but significant (0(0^++) and 1(1^+-))
contributions from the hadrocharmonium component. The situation for open charm is similar; the dominant
contribution is the heavy-light meson component although significant corrections from the diquark-antidiquark
component arise. Although the masses of the charm quarks are far from static, mass
patterns expected from heavy-quark symmetry are visible: In the open-charm sector we observe the expected
mass hierarchy and difference between the I=0 and I=1 axialvector channels <cit.>.
The observed mass pattern in the hidden charm sector resembles the one expected of the lowest-lying multiplet of
states in the hadronic molecular approach <cit.>. In order to make further contact
with heavy-quark symmetry and lattice QCD, it would be interesting to further increase the masses of the heavy
quarks and explore the bottomonium sector. This is technically challenging and therefore left for future studies.
A potential caveat of the present formalism is that it does not take into account potentially important effects
from mixing with ordinary cc̅ states in the I=0 channels <cit.>.
Again, this is left for future work.

Acknowledgments.
We are grateful to Christoph Hanhart, Soeren Lange, Sasa Prelovsek and Marc Wagner for discussions.
This work was supported by the DFG grant FI 970/11-1 and by the FCT Investigator Grant IF/00898/2015.



apsrev4-1





   


 




§ CONSTRUCTION OF THE FOUR-QUARK BETHE-SALPETER AMPLITUDE
 


  In this appendix we describe the construction of the tetraquark amplitude Γ^(μ)_αβγδ(p_1, p_2, p_3, p_4) in detail.
  It is the direct product of Dirac, color and flavor parts; here, {α, β, γ, δ}
  either represent the Dirac indices or absorb all Dirac, color and flavor indices to avoid clutter.
  If we consider the color parts alone, we work with latin indices {A,B,C,D}.
  We also suppress the Lorentz index μ that distinguishes between J=0 and J=1 states.

  We consider the three different physical systems shown in Fig. <ref>:
  a state with four identical quarks qqqq,
  a hidden-charm system cqqc,
  and an open-charm state ccqq.
  Depending on how we connect the indices, each system
  can come in a diquark-antidiquark () configuration (12)(34) and two meson-meson configurations (13)(24) and (14)(23).
  In terms of diquark and meson clusters,
  the hidden-charm system cqqc contains
  two heavy-light diquarks in (12)(34), whereas (13)(24) describes a `molecule'
  of two heavy-light mesons and (14)(23) a `hadrocharmonium' with a heavy and a light meson.
  For the open-charm system ccqq, (12)(34) contains a heavy diquark and a light antidiquark
  and the two meson-meson configurations are both molecular.


  Depending on which system we study, the amplitude is subject to certain symmetry constraints, namely Pauli antisymmetry in (12) or (34)
  and charge-conjugation symmetry in (13)(24) or (14)(23). Pauli antisymmetry in (12) or (34) means
  
    Γ(p_2, p_1, p_3, p_4)_βαγδ!=   - Γ(p_1, p_2, p_3, p_4)_αβγδ , 
    Γ(p_1, p_2, p_4, p_3)_αβδγ!=   - Γ(p_1, p_2, p_3, p_4)_αβγδ ,





  where a permutation of all (Dirac, color and flavor) indices is understood. Likewise,
  charge-conjugation symmetry in (13)(24) amounts to
  
    C_αα' C_ββ' C_γγ' C_δδ' Γ(-p_3,-p_4,-p_1,-p_2)_γ'δ'α'β'
               !=±Γ(p_1, p_2, p_3, p_4)_αβγδ







  and in (14)(23) to
  
    C_αα' C_ββ' C_γγ' C_δδ' Γ(-p_4,-p_3,-p_2,-p_1)_δ'γ'β'α'
               !=±Γ(p_1, p_2, p_3, p_4)_αβγδ ,







  where C = γ_4 γ_2 is the charge-conjugation matrix and the signs ± determine the C parity of the state.
  The amplitude for a qqqq system satisfies all four relations, whereas
  a hidden-charm amplitude cqqc is only subject to Eq. (<ref>) and
  an open-charm configuration ccqq only satisfies Eqs. (<ref>–<ref>).

  To proceed, we write the most general tensor basis decomposition of the amplitude as
  
    Γ(p_1 … p_4)_αβγδ = ∑_i f_i(…) τ_i(p_1 … p_4)_αβγδ ,

  where the Lorentz-invariant dressing functions f_i(…) depend on the ten Lorentz invariant momentum variables
  that can be constructed from four independent momenta. The tensors τ_i are the direct products of Dirac, color and flavor parts.

  Let us first work out the color tensors.
  From 3⊗3⊗3̅⊗3̅ =
       ( 3̅⊕6) ⊗ (3⊕6̅) = 1⊕1⊕ ..., the
  color part of the amplitude
  consists of two independent color singlet tensors, which can be taken from
  the  (3̅⊗3, 6⊗6̅)
  or either of the meson-meson configurations (1⊗1, 8⊗8),
  for example:
  
    (_11)_ABCD   =  1/3 δ_AC δ_BD
    
         (_11')_ABCD   =  1/3 δ_AD δ_BC .

  The two tensors in the  decomposition are linear combinations of these,
  
    _3̅3 = -√(3)/2 (_11-_11') ,   _66̅ = √(3/8) (_11+_11') ,

  as well as the remaining octet-octet tensors:
  
    _88 = _11 - 3 _11'/2√(2) ,   _88' = _11' - 3 _11/2√(2) .

  The tensors {_11, _88}, {_11', _88' } and {_3̅3, _66̅}
  are mutually orthogonal.
  The color tensors are invariant under both charge-conjugation operations (<ref>–<ref>), i.e., they carry a positive sign, whereas
  either of the two Pauli symmetries (<ref>) or (<ref>) transforms them into each other: _11↔_11'
  and therefore _3̅3↔ -_3̅3.

  Next, we write down the Dirac-color tensors. A scalar tetraquark has 256 linearly independent Dirac tensors
  (see Appendix A in Ref. <cit.> for the full basis construction),
  which together with the two independent color structures amounts to 512 tensors in total:
  
    ϕ_1    = γ^5_αγ γ^5_βδ _11 , 
    ϕ_2    = γ^5_αδ γ^5_βγ _11' , 
    ϕ_3    = γ^μ_αγ γ^μ_βδ _11 , 
    ϕ_4    = γ^μ_αδ γ^μ_βγ _11' , 
    ϕ_5    = (γ_5 C)_αβ (C^T γ_5)_γδ _3̅3 , 
    ϕ_6    = (γ^μ C)_αβ (C^T γ^μ)_γδ _3̅3 , 
         ⋮

  Here the Greek subscripts are Dirac indices only and we suppressed the color indices as well as those for the full Dirac-color amplitudes ϕ_i;
  they are always understood to be in the standard order αβγδ or ABCD.


  The tensors in (<ref>) correspond to the dominant `physical' two-body clusters
  in terms of mesons and diquarks:
  ϕ_1 and ϕ_2 describe configurations with two pseudoscalar mesons in (13)(24) and (14)(23),
  ϕ_3 and ϕ_4 represent vector-vector meson configurations,
  ϕ_5 a (12)(34) diquark-antidiquark system with two scalar diquarks (SS)
  and ϕ_6 one with two axialvector diquarks (AA).
  They are part of a Fierz-complete subset of 16 Dirac tensors that do not depend on any relative momentum,
  i.e., the `s waves' which do not carry orbital angular momentum.
  In principle one could continue the list until the basis is complete.



  Similarly, an axialvector tetraquark with J^P=1^+ has 3× 256 = 768 linearly independent Dirac tensors (48 of which are s waves),
  which amounts to 1536 Dirac-color tensors in total <cit.>:
  
    ψ_1^±   = (γ^5_αγ γ^μ_βδ±γ^μ_αγ γ^5_βδ) _11 , 
    ψ_2^±   = (γ^5_αδ γ^μ_βγ±γ^μ_αδ γ^5_βγ) _11' , 
    ψ_3    = ε^μνρσ P̂^ν γ^ρ_αγ γ^σ_βδ _11 , 
    ψ_4    = ε^μνρσ P̂^ν γ^ρ_αδ γ^σ_βγ _11' , 
    ψ_5    = (γ_5 C)_αβ (C^T γ^μ)_γδ _3̅3 , 
    ψ_6    = (γ^μ C)_αβ (C^T γ_5)_γδ _3̅3 , 
    ψ_7    = ε^μνρσ P̂^ν (γ^ρ C)_αβ (C^T γ^σ)_γδ _3̅3 , 
         ⋮

  We suppressed again the indices of the ψ_i and color tensors.

  ψ_1^±, ψ_2^± are pseudoscalar-vector and
  ψ_3, ψ_4 vector-vector meson configurations,
  ψ_5, ψ_6 describe scalar-axialvector diquark configurations (SA)
  and ψ_7 is one with two axialvector diquarks (AA).




  In Table <ref> we collect the combinations of the Dirac-color tensors ϕ_i and ψ_i that carry definite Pauli symmetry.
  For example, a combination of type Φ_ss is fully symmetric under both (12) and (34),
  Φ_sa is symmetric under (12) but antisymmetric under (34), etc.
  Table <ref> lists the transformation properties of the ϕ_i and ψ_i under charge conjugation.

  The symmetries of the Dirac-color tensors dictate the symmetries of the remaining flavor tensors.
  As an example, consider the leading Dirac-color-flavor tensors (γ_5 C) _3̅3 _0
  for a scalar diquark and (γ^μ C) _3̅3 _1 for an axialvector diquark.




  Here, _0=[ud] is the I=0 flavor wave function and _1 denotes the isospin triplet consisting of uu, {ud} and dd, where
  […] means antisymmetrization and {…} symmetrization.
  The color tensor _3̅3 is Pauli-antisymmetric and the Dirac and flavor parts are either both antisymmetric or both symmetric,
  which ensures the antisymmetry of the total wave function. As a consequence, scalar diquarks carry isospin I=0 and axialvector diquarks I=1.
  One can also construct diquarks with reversed isospin by
  including Pauli-antisymmetric momentum prefactors p_1^2-p_2^2 (in analogy to quantum-number exotic mesons), but those
  prefactors would induce a strong suppression of the amplitudes.


  Following this argument, we restrict ourselves to dressing functions
  f_i(…) in Eq. (<ref>) that are even under all symmetries, so that
  the symmetry properties are carried by the basis tensors τ_i alone.

  The flavor wave functions depend on the explicit quark content (qqqq, ccqq or cqqc) and determine the isospin of the state.
  They also transform under the respective symmetries.
   Then, for a system which has Pauli symmetry, the Dirac-color amplitudes in Table <ref> must be combined
  with appropriate flavor wave functions  to form a totally antisymmetric tensor:
  
    [ Φ_ss _aa ,; Φ_aa _ss , ]    [ Ψ_ss _aa ,; Ψ_aa _ss , ]    [ Ψ_sa _as ,; Ψ_as _sa . ]

  If C parity does not apply (such as for ccqq), these are  the final Dirac-color-flavor wave functions;
  otherwise (qqqq) one must find appropriate combinations which also have definite C parity.
  If the system has no Pauli symmetry and is only subject to charge-conjugation invariance in one channel
  (cqqc), one can combine the ϕ_i and ψ_i directly with the flavor tensors to obtain total wave functions with definite C parity.











  



  



  The resulting Dirac-color-flavor tensors that can be constructed in this way are listed 
  in Table <ref>, together with the isospin and C parities they carry.
  In the following we discuss the construction for different quark content in detail.




  


 

  ▪ nnnn: A system made of four light quarks n ∈{u,d} satisfies all symmetries (<ref>–<ref>).
  From 2⊗2⊗2̅⊗2̅= (1_a ⊕3_s) ⊗ (1_a ⊕3_s)
  one obtains 16 flavor wave functions with definite Pauli symmetry, which can be arranged into isospin multiplets:
  
    _aa     …   I=0 , 
    _ss     …   I = 0, 1, 2 , 
    _as     …   I=1 , 
    _sa     …   I=1 .

  For example, _aa = [ud][u̅d̅] is the I=0 flavor wave function typically associated with the σ meson.
  _aa and the neutral members of the multiplets for _ss have positive C parity under any of the two charge-conjugation operations,
  whereas the remaining ones transform into each other: _as↔_sa under (13)(24) and _as↔ -_sa
  under (14)(23).



  The C parities of the full wave functions Φ_ss _aa and Φ_aa _ss from Eq. (<ref>) can be read off from Table <ref>.
  Since all ϕ_i as well as _aa and _ss have positive C parity,
  the same is true for their combinations and thus they describe I(0^++) states, where the isospin depends on the flavor wave functions.
  Similarly, ψ_1,2^+ and ψ_7 have negative C parity so that the combinations Ψ_ss _aa, Ψ_aa _ss describe I(1^+-) states.
  The remaining wave functions Ψ_sa _as and Ψ_as _sa carry I=1 due to their flavor parts;
  their combinations with definite C parity are given by
  
    Ψ_± = {[ (ψ_1^–ψ_2^-) F_as∓ (ψ_1^-+ψ_2^-) F_sa ,;        (ψ_3-ψ_4) F_as± (ψ_3+ψ_4) F_sa ,;                     ψ_5 F_as±ψ_6 F_sa , ].

  where Ψ_± corresponds to I(J^PC)=1(1^+±).

  The final Dirac-color-flavor wave functions satisfy all four constraints in Eqs. (<ref>–<ref>).
  Each has an interpretation in terms of a meson-meson or  configuration. The corresponding physical
  two-body systems are listed in the right column of Table <ref>.






  

  ▪ cccc: For an all-charm state there is only one flavor wave function _ss = cccc
  which is symmetric under any quark exchange, invariant under both charge-conjugation operations and carries I=0.
  As a consequence, only the combinations Φ_aa _ss and Ψ_aa _ss in Eq. (<ref>) survive
  and produce states with quantum numbers 0(0^++) and 0(1^+-), respectively.
  Note that due to the symmetry of _ss
  there is no overlap with the `σ-like' components Φ_ss in the nnnn system.

  

  ▪ ccnn: For an open-charm system, the light quarks can be arranged into
  an isosinglet [u̅d̅] and an isotriplet u̅u̅, {u̅d̅}, d̅d̅.
  Thus there are four flavor wave functions:
  
    _sa     …   I=0 , 
    _ss     …   I=1 .

  In this case charge conjugation symmetry does not apply, so that the final wave functions are
  Φ_aa _ss, Ψ_aa _ss, Ψ_as _sa with quantum numbers 1(0^+), 1(1^+) and 0(1^+), respectively.
  Note that this does not produce a 0(0^+) state since there is no entry Φ_as in Table <ref>.

  

  ▪ ccss: Here the only flavor wave function is _ss, however with I=0,
  so that the combinations Φ_aa _ss and Ψ_aa _ss produce the quantum numbers 0(0^+) and 0(1^+).
  Compared to the previous case, going from light to strange quarks for the same quantum numbers therefore leads to quite different wave-function components.


  

  ▪ cnnc: In the hidden-charm system Bose symmetry does not apply.
  Instead, the flavor wave functions
  
    _0    = c (uu̅+dd̅) c̅        …   I=0 , 
    _1    = {[         cud̅c̅; c (uu̅-dd̅) c̅;        c du̅c̅ ]}  …   I=1

  have positive C parity in (14)(23).
  The C parities of the full wave functions {ϕ_i, ψ_i}×{_0, _1} can then be read off
  from the bottom line of Table <ref>: The Dirac-color amplitudes
  
    Φ_1    = {ϕ_1,  ϕ_2,  ϕ_3,  ϕ_4,  ϕ_5, ϕ_6 } , 
    Ψ_1    = {ψ_1^-, ψ_4,  ψ_5-ψ_6 } , 
    Ψ_2    = {ψ_1^+,  ψ_2^+ ±ψ_2^-,  ψ_3,  ψ_5+ψ_6,  ψ_7 }

  correspond to the quantum numbers I(0^++), I(1^++) and I(1^+-), respectively,
  where the isospin I=0,1 depends on the flavor wave function (<ref>).

  

  ▪ cssc: Here the situation is analogous except there is only one
  flavor wave function _0 = css̅c̅ with I=0.

  

  Up to this point we have made two approximations in the construction of the Bethe-Salpeter amplitude.
  Instead of the complete tensor basis, we have restricted ourselves to the dominant tensors in Eqs. (<ref>–<ref>),
  and we only considered dressing functions f_i that are even under the relevant symmetry operations.
  Those dressing functions, however, still depend on ten independent variables: the total momentum P^2 = (p_1 + p_2 + p_3 + p_4)^2 = -M^2,
  where M is the mass of the tetraquark, plus nine variables that are dynamical.
  This is what makes a full numerical solution practically impossible.

  In the final step we therefore assume that the momentum dependence of the f_i is mainly carried by the symmetric variable
  
    _0 = 1/4(p_1^2 + p_2^2 + p_3^2 + p_4^2 + M^2/4),

  and that the remaining angular dependencies can be well described by two-body intermediate poles
  corresponding to the physical components listed in Table <ref>. Such a behavior has been found
  directly in the dynamical solution of the four-body equation for light scalar tetraquarks <cit.>
  and it is also the approximation employed in Ref. <cit.>.
  We therefore assume
  
    f_i(…) ≈ f_i(_0) (m_1, m_2)_ab,cd ,

  where the pole factors depend on the respective tensor component:
  
    (m_1,m_2)_ab,cd = 1/(p_a+p_b)^2 + m_1^2 1/(p_c+p_d)^2 + m_2^2 .

  For example, for a cnnc system with I=0, the pole factors corresponding to the ϕ_i in Eq. (<ref>) are
  
    (m_D,m_D)_13,24 , 
       (m_η_c,m_η)_14,23 , 
       (m_D^∗,m_D^∗)_13,24 , 
       (m_J/ψ,m_ω)_14,23 , 
       (m_S,m_S)_12,34 , 
       (m_A,m_A)_12,34 ,

  where we compute all meson and diquark masses from their two-body Bethe-Salpeter equations.
  Note that these pole factors are invariant under the symmetry operations (<ref>–<ref>)
  that are relevant for the particular system (in the cnnc example
  only charge conjugation in (14)(23) applies), so we can also attach them directly to the tensors (<ref>–<ref>).
  As a result, only the dependence on the symmetric variable _0 remains dynamical.

  The f_i(_0) are therefore the amplitudes
  that are determined in the solution of the four-body equation shown in Fig. <ref>, together with the tetraquark mass M.
  By switching off individual tensor components we can determine their strengths and isolate the dominant contributions.
  This completes our construction of the tetraquark Bethe-Salpeter amplitude.

 

















 




