Hidden-flavour four-quark states in the charm and bottom region

We discuss the spectrum and the internal composition of ground and excited four-quark states in the charm and bottom energy region. To this end we extend previous calculations within the framework of the relativistic four-body Faddeev-Yakubovsky equation to include quantum numbers with $J^{P C} = 0^{++} , 0^{-+} , 1^{--} , 1^{+-}$ and $1^{+ +}$ and study their internal composition in terms of heavy-light meson pairs, hadroquarkonia and diquark-antidiquark clusters. We observe similar patterns in the charm and bottom energy region with different compositions of the four-quark states depending on $J^{P C}$ quantum numbers. Most notably, we find that all states with $C \cdot P = +1$ are dominated by heavy-light meson contributions, whereas for axialvector states with $J^{P C} = 1^{+-}$ including the $Z_c (3900)$ we find a much more complicated picture depending on the flavour content. We systematically compare our results for the spectrum with existing experimental results and provide predictions for future analyses.


I. INTRODUCTION
Starting out with the unexpected detection of a narrow state in the J/ψπ + π − invariant mass spectrum by the Belle collaboration in 2003 [1], i.e., the χ c1 (3872), an ever increasing number of exotic states has been identified in the charmonium and bottomonium mass region by Belle, BaBar, BSE III and the LHC experiments in the last two decades.Many of these 'XY Z states' cannot be accommodated for in the conventional quark model for Q Q mesons (with Q = c, b) and are therefore generally referred to as exotic hadrons.Some of these states carry electric charge, which can only be explained assuming four valence (anti-)quarks with either hidden or open heavy flavour configurations, i.e.Q Qq q or QQq q, Q Qqq (with q = u, d, s).Thus, four-quark states are considered to be good candidates to study the properties of these exotic hadrons, see, e.g., [2][3][4][5][6][7][8][9] for recent reviews.
A highly debated and unsettled property of four-quark states is their internal structure.In most effective field theory and model approaches, one can generally distinguish between three different a priori assumptions regarding a possible internal clustering.The first possibility is motivated by the experimental observation of final states with a specific charmonium/bottomonium state and light hadrons.This is known as the hadroquarkonium picture [10] and features a tight heavy quark and antiquark (Q Q) core which is surrounded by the light q q pair.The second prominent possibility is the clustering of the constituents into diquark-antidiquark (dq − dq) pairs which interact via colour forces, see, e.g., [2] for a review.Finally, there is the meson-molecule picture, which is especially relevant for states close to open-flavour thresholds.In this picture the constituents arrange in pairs of D ( * ) D( * ) or B ( * ) B( * ) mesons which interact via short-and/or longrange forces [6].
On a general note, these three possibilities of internal clustering are not mutually exclusive, i.e., experimental states may be a superposition of all three different clusters with the 'leading' component possibly different on a case-by-case basis.To thoroughly study this behaviour, it is important to develop theoretical approaches to QCD that can deal with all three possibilities.A prominent approach is lattice QCD, where interesting progress has been made in recent years, see [11][12][13][14][15][16][17][18][19] and references therein.Recently, the framework of functional methods has been generalized to systematically investigate four-quark states with any J P C and flavour combination, see [20] and references therein for a review on results in the light and charm mass region.
In this work we reanalyse and confirm the findings for the hidden-charm four-quark states in [21], i.e., states with quantum numbers I(J P C ) = 0(1 ++ ), 1(1 +− ) and 0(0 ++ ).Furthermore, we investigate the experimentally very interesting vector (1 −− ) channel in the hidden-charm region, and discuss our findings for possible pseudoscalar (0 −+ ) four-quark states in that region, where there are currently no experimental exotic candidates.We also present our novel results for the bq qb hidden-bottom mass spectrum with the aforementioned quantum numbers.Last but not least, we discuss a new method to investigate the internal structure of four-quark states by considering how much each of the three internal clusters contributes to the normalization of a given state [22,23].
The paper is organized as follows: In Section II we briefly introduce the four-body Bethe-Salpeter equation (BSE) and discuss the technical details, i.e., truncation of the two-body interaction, construction of the physical basis and the internal structure.In Sec.III we present and discuss our results for the hidden-charm and hiddenbottom four-quark states, before we conclude in Sec.IV.

II. SETUP A. Four-quark Bethe-Salpeter equation
In this work we focus on heavy-light and heavy-heavy four-quark states with hidden flavour.We denote their quark content by Qq q Q with q ∈ {u, d, s, c, b} and Q ∈ {c, b}.Any bound state or resonance in QCD that has overlap with Qq q Q appears as a pole in the quark eightpoint correlation function or scattering matrix T (4) , which satisfies the scattering equation T (4) = K (4) + K (4) G (4) 0 T (4) . (1) Here, K (4) is the four-body interaction kernel and G (4) 0 denotes a product of four dressed (anti-)quark propagators; see [24] for details.From the pole residue of Eq. ( 1) one obtains the homogeneous four-quark BSE, written in compact notation as 0 Γ (4) . ( Each multiplication implies an integration over all loop momenta, and Γ (4) is the four-quark Bethe-Salpeter amplitude of a given state.Eq. ( 2) is an eigenvalue equation for K (4) G (4) 0 with eigenvalues λ i (P 2 ), which depend on the total hadron momentum squared P 2 ∈ C. If the condition λ i (P 2 i ) = 1 is satisfied, this corresponds to a pole in the scattering matrix at The index i indicates whether the eigenvalue satisfying the condition corresponds to the ground state (i = 0), the first radial excited state (i = 1), etc. Below a given meson-meson threshold, M i is real and we have found a bound state.For a resonance, on the other hand, this condition is only satisfied in the complex plane on a higher Riemann sheet.In principle the homogeneous BSE is able to detect both bound states and resonances, where contour deformations are required to calculate λ i (P 2 ) above the lowest threshold [25][26][27].
The scattering kernel K (4) consists of irreducible two-, three-and four-body correlations.In this work we primarily want to study the internal two-body clusters, hence we neglect the three-and four-body forces.The resulting kernel is then the sum of the two-body interactions where a and a ′ denote interactions between quark-(anti-) quark pairs and aa ′ is one of three combinations (12)(34), (13)( 24) and ( 14) (23).These correspond to the two-body interaction topologies for quark content Qq q Q: Diquarkantidiquark (Qq)(q Q), meson-meson (Qq)(q Q) and hadroquarkonium (Q Q)(q q).Note that the last term on the r.h.s. of Eq. ( 3) is necessary to avoid overcounting [28][29][30][31].
The resulting four-quark BSE with the kernel in Eq. ( 3) is shown in Fig. 1.
For the two-body kernels, we employ the rainbow-ladder truncation in combination with the effective Maris-Tandy FIG. 1. Four-quark BSE for a generic hidden-flavour Qq q Q system in the (12)(34) topology; the permutations (13)( 24) and ( 14) (23) are not shown here.The green half-circles denote the Bethe-Salpeter amplitudes, blue boxes represent the twobody interaction kernels and the blobs denote fully-dressed quark propagators.
(MT) interaction [32,33].This setup models the combined effect of the dressed gluon propagator and quarkgluon vertex and has been extensively applied to meson, baryon and four-quark phenomenology, see [20,24] for recent reviews.The explicit form of the interaction can be found in Eq. (3.96) of [24]; we use the scale parameter Λ = 0.72 GeV tuned to reproduce the pion decay constant and fix the shape parameter to η MT = 1.8.The Dyson-Schwinger equation (DSE) for the quark propagator, which is needed as an input for the BSE (cf.Fig. 1), is solved using the same interaction.The MT interaction is known to describe the phenomenology of light mesons in the pseudoscalar and vector channel reasonably well, whereas in the scalar and axialvector channels it does not provide satisfactory results.Its qualitative reliability can be judged from the meson masses in Table I obtained via the quark-antiquark BSE.
We work in the isospin symmetric limit, i.e., m s and m D + m D * have < 0.5% deviation from the sums of the respective experimental values [34].The bottom quark mass is fixed to m b = 3750 MeV (again same scheme) such that the pseudoscalar and vector b b masses and the sum m B + m B * matches experiment within 0.8% relative error.With these constraints set, the deviations between the theoretical and experimental meson masses are then below 7% in most channels.

B. Physically motivated four-quark amplitude
In general, the four-quark Bethe-Salpeter amplitude is a direct product of Dirac (D), colour (C) and flavour (F) parts.For a given J P C it can be written as with p 1 , . . ., p 4 denoting the quark momenta.The µν . . .occur as Lorentz indices for states with higher spin J.We show the experimental candidates [34], the masses mRL obtained in our rainbow-ladder calculation, and the relative error of these to the masses given in the PDG (if the experimental state has been identified).In the last two columns we also show the obtained rainbow-ladder masses for the corresponding Qq diquarks with quantum numbers J P = {0 + , 1 + }.All values are given in MeV.† : The π and the η are mass degenerate in this work, since we neglect the strange component in the η and the indirect effect of the topological mass via octet-singlet mixing.‡ : We do not consider the lightest scalar meson nonet as potential internal components, since they themselves are of four-quark nature [35,36].Instead, we resort to the scalar nonet with masses above 1 GeV.
The colour and flavour part of the amplitude are straightforward to work out and discussed in detail in the supplemental material of [21].
The Dirac part can be written as where the f i (Ω) are Lorentz-invariant dressing functions depending on ten Lorenz invariant momentum variables Ω (see [37] for details).The τ (µ) i are the corresponding Dirac structures (with Dirac indices α, β, γ , δ) and N is the number of Dirac basis elements.The full Dirac bases for J = 0 (N = 256) and J = 1 (N = 768) states are collected in [36] and in the Appendix of [37].
Following the arguments in [21,36], we note that the amplitude dynamically develops two-body clusters in the three different topologies mentioned in Section II A. For heavy-light systems, this is a heavy-light meson-meson component (M 1 ), a hadro-quarkonium component (M 2 ) and a dq − dq cluster (D).These clusters were found to heavily influence the four-body system, thus the guiding idea is to express Γ (µ) (p 1 , p 2 , p 3 , p 4 ) in terms of these internal two-body clusters.
Applying this idea to Eq. ( 4), we construct a physically motivated basis for the Bethe-Salpeter amplitude by projecting onto a subset of the full basis corresponding to the dominant two-body clusters.For given quantum numbers J P C , we draw on existing information on the decay channels of experimental four-quark candidates to identify the possible two-body clusters.This in turn fixes the Dirac tensors that enter in the amplitude.A collection of the chosen internal configurations used in this work is given in Table II.With all of the above, the amplitude in Eq. ( 4) reduces to where the sum includes the dominant physical components for the different interaction topologies in Table II.The resulting structure is visualized in Fig. 2. We emphasize that this representation is different from the two-body framework used in [38], as the sub-amplitudes depicted in Fig. 2 are not actual two-body amplitudes but rather specific Dirac-colour-flavour basis elements of the fourbody amplitude that reflect its internal clustering.In Eq. ( 6), τ C i is the attractive colour singlet structure corresponding to each interaction topology, with 1 ⊗ 1 tensors for each of the two meson-meson topologies and a 3 ⊗ 3 tensor for the diquark-antidiquark topology (see Supplemental Material of [21] for details).We relegate the inclusion of repulsive colour-singlet channels to future work.The index i on the flavour part τ F i is only used to construct the Dirac and colour part of the wavefunction but is otherwise irrelevant for the calculation, as the interaction kernel described in Section II A is flavourblind.

C. Two-body poles
Let us now take a closer look at the Dirac part of the amplitude.One can express the quark momenta p 1 , . . ., p 4 by three relative momenta k, q, p and the total hadron  and A c/b , respectively, where the subscript characterizes the heavy quark that is paired with the light quark.The f0 here denotes the f0(1370).The notation f0, f1, f2 for the physical components corresponds to Eq. ( 8) and is used again to display results in subsection III D.
momentum P using the relations Here, σ 1 , . . ., σ 4 are (quark-)momentum partitioning parameters which can be used to maximize the calculable domain for the bound-state mass by optimally distributing the total hadron momentum among the quarks.In the rest frame of the four-quark state, P µ = (0, 0, 0, iM ) is imaginary whereas the relative momenta p, q, k are real.Therefore, the p 2 i obtained from Eq. ( 7) describe different parabolas in the complex plane which are limited by the nearest quark singularities, which translates into an upper limit for M .An optimized choice of momentum partitionings then maximizes the mass range M for which the BSE can be solved.
Using the relations in [39], one can group the Lorentzinvariant momentum variables Ω = {q 2 , p 2 , k 2 , . ..} in multiplets of the permutation group S 4 .This yields a singlet variable S 0 = (k 2 + q 2 + p 2 )/4 carrying the momentum scale, a doublet containing the internal two-body clusters, and two triplets, thus totalling to 3+3+2+1 = 9 momentum variables, plus P 2 = −M 2 .Following [36,37], the leading momentum dependence of the dressing functions f i beyond the singlet variable S 0 comes from the two-body clusters, whose poles are dynamically generated when solving the equation.We therefore pull out these poles from the f i and, to reduce computational effort, assume that the remainder only depends on S 0 .The resulting Dirac part of the amplitude then reads where the two-body poles of the amplitude are given by Here, p + ab = p a + p b is the momentum of the meson or diquark with mass m ab in a given topology (ab)(cd) = (13)( 24), ( 14) (23) or ( 12) (34).The sum in Eq. ( 8) runs over the physical components given in Tab.II.These occur in the diagrammatic topologies discussed in Fig. 2. For example, for a four-quark state with I(J P C ) = 0(1 ++ ) the M 1 topology has clusters m 13 = m D and m 24 = m D * , the M 2 cluster m 14 = m J/ψ and m 23 = m ω , and the D topology m 12 = m Scq and m 34 = m Aqc .On the other hand, the first two states of Tab.II would have a single contribution (f 0 ) in M 1 topology, two contributions (f 1,2 ) in M 2 and none in D. In general, the respective meson and diquark masses are calculated from their rainbowladder BSEs as described above and compiled in Table I.
Eq. ( 9) introduces two-body poles in the integration domain of the r.h.s. of Eq. ( 2), thus restricting the range of P 2 = −M 2 where the BSE can be solved directly.This can be somewhat remedied by optimizing the quark momentum partitioning parameters in Eq. ( 7) which split the hadron momentum P amongst the quarks.On the other hand, the momenta in the denominator of Eq. ( 9), which are linear combinations of Eq. ( 7), also form complex parabolas which are limited by the meson and diquark poles.To this end, it is advantageous to relate the σ i to new parameters η, ζ and χ, where η corresponds to the meson-meson topology (M 1 ), ζ to the hadro-quarkonium (M 2 ) and χ to the dq − dq (D) cluster: The choice Current-mass evolution of the cq qc (crosses) and bq qb ground states (dots) in the 0(1 ++ ) and 1(1 +− ) channels.The grey vertical dashed lines mark the position of the q = n, s, c, b current-quark masses, where n = u/d.The dotted, dashed and dash-dotted curves represent the meson-meson and diquark-antidiquark thresholds for the charm and bottom four-quark system, respectively.The colours of the thresholds comply with Fig. 2, i.e., blue for the M1, purple for the M2 and orange for the D cluster.The grey bands are the fits to the data.
then maximizes the value of M to be the lowest sum of masses of the individual physical components in Table II, e.g., m D + m D * for the χ c1 (3872).Here we also remedy a slight inconsistency in the previous works [21,37], where the quark momentum partitioning parameters were chosen as σ i = 1 4 in the pole terms (9) but set to their optimal values in the rest of the equation.As a result, the maximum value of M did not exhaust its full range, which resulted in the need for extrapolating the eigenvalue curves over a large momentum range causing a large extrapolation error.In the present work we overcome this limitation, thereby reducing the extrapolation error considerably.
To calculate the eigenvalues above the thresholds, in principle one needs to employ contour deformation techniques or elaborate analytic continuations [26,27].As the former would present an enormous technical challenge in the four-body approach, we relegate it to future work and analytically continue the eigenvalues on the real axis using the Schlessinger-point method [40], see Sec.B for details.Thus, if masses above thresholds are quoted they merely serve as a rough estimate for the real part of the corresponding resonance locations.In many cases, however, the resulting ground-state masses are below all thresholds and thus no analytic continuation is needed.

III. RESULTS
In the following we present results for hidden-charm and hidden-bottom four-quark states with quantum numbers 0(0 ++ ), 0(1 ++ ), 1(1 +− ), 0(1 −− ), 0(0 −+ ).The first three of these were already investigated in the hidden-charm sector in [21,37].Here we extend the calculations to the bottom region and also present novel results for the vector and pseudoscalar channels.We first present our mass evolution curves in section III A and then discuss our results for the physical spectrum in the charm and bottom energy region in sections III B and III C. The internal composition of our states is then the topic of section III D.

A. Mass evolution curves
We first discuss the charm-and bottom-like ground states in the I(J P C ) = 0(1 ++ ), 1(1 +− ) and 0(0 ++ ) channels.While the 0(1 ++ ) channel has established experimental four-quark candidates only in the charm region [34], e.g., χ c1 (3872) and χ c1 (4140), the 1(1 +− ) channel features four-quark candidates both in the charm region (e.g., the Z c (3900)) and the bottom region, namely the Z b (10610) and Z b (10650).The situation in the 0(0 ++ ) channel is not yet settled in the literature, but the existing exotic candidates also only occur in the charm region.
We display the results of our calculation for these three channels in Figs. 3 and 4, where we show the mass evolution curve (MEC) of the four-quark state with fixed heavy-quark pair Q Q = cc (lower group of curves) and b b (upper group of curves).The quark pair q q is varied from the bottom mass m b (rightmost vertical dashed line) to the light quark mass m n (leftmost vertical dashed line).The results for charm-like states are marked by crosses and those for bottom-like states by dots.The dotted, dashed and dash-dotted curves show the quark-mass evolution of the two-body thresholds (cf.Table II).We show the MECs for the first radial excited states in Appendix C.
We find that for increasing current-quark masses the four-quark states become more deeply bound with respect to the lightest meson-meson threshold.For quark masses m q ≥ m c , the MECs become approximately linear both in the charm and bottom region.However, below the charm mass we observe an upwards bending of the MEC when it approaches the lowest meson-meson threshold in the system, which is stronger for the bq qb compared to the cq qc states.We note that such a bending is also observed in the MECs for the two-body heavy-light states in our framework.In the pseudoscalar channel, the MECs bend downwards for m q ≃ m n , resembling the observed behaviour of the MECs for the two-body q q states in the pseudoscalar (and scalar) channel, see, e.g., right figure in Fig. 3.10 in [24].Since the masses from the aforementioned two-body MECs serve as input for the four-quark state calculations, the similarity in the behaviour of the MECs could be interpreted as a first indication that the states in the pseudoscalar channel have a strong hadroquarkonium component.A more detailed discussion about the internal structure can be found in Sec.III D.
In Tables III and IV we quote the masses of all ground and excited states calculated in this work.To this end, we employ fits of the form to the MECs at quark masses reasonably far away from the two-body thresholds.These fits are shown in Figs. 3 and 4 by the grey bands.The data not taken into account in the fits are depicted slightly opaque.The error of the masses, given in the brackets in Table III, is then determined by combining the width of the bands (i.e., the error of the fit function M (m q )), with the extrapolation error as described in Appendix B. To facilitate comparisons with the literature, we also list the resulting binding energies E B = M − M th with respect to our calculated lightest heavy-light meson-meson threshold M th in both tables.For the pseudoscalar channel, we instead use the dominant hadro-quarkonium threshold, see Sec.III D. Note that these thresholds do not necessarily coincide with the experimental ones given that our meson masses are calculated in the same framework as our four-quark masses.
We now move on to the novel results of this work, i.e., the vector 0(1 −− ) and pseudoscalar 0(0 −+ ) four-quark states.For the 0(1 −− ) channel there are many established exotic candidates in the charm and bottom mass region [34]: the ψ(4230), ψ(4360), ψ(4660), Υ(10753), Υ(10860), and Υ(11020).By contrast, as of today there are no experimental exotic candidates in the pseudoscalar channel.Thus, we can only compare our results with the experimental situation in the vector channel and crosscheck our obtained results with theoretical predictions for the pseudoscalar channel.
In contrast to the other three quantum numbers, our "physical basis" for the vector and pseudoscalar four-quark states does not contain a diquark topology but rather a second M 2 cluster (cf.Table II).The reason for this is that, on one hand, the diquark clusters always constitute the highest two-body thresholds in our system and are generally found to be almost negligible for hidden-flavour four-quark states.On the other hand, the construction of S-wave dq − dq pairs in the 0(1 −− ) and 0(0 −+ ) channels would require including pseudoscalar and vector diquarks, which are strongly suppressed compared to their ("good") scalar and ("bad") axialvector diquark counterparts (see Appendix A for details).
Our results for the ground states of the vector and pseudoscalar states are displayed in Fig. 5.In both channels, the behaviour of the MECs are very similar to the ones described before, i.e., the results are again affected by the meson-meson thresholds for light quark masses.
As mentioned in Sec.II A, we calculate both the fourquark ground states and radial excitations from Eq. ( 2).The resulting masses of the radial excitations are given in Table IV.Note that because the lowest two-body thresholds are identical for the ground and excited states, the eigenvalue curves need to be extrapolated much further in some cases and thus the errors for these states increase.Most of the excited states are therefore also unbound resonances.
As noted before, the binding energies in Tables III and IV are determined with respect to our calculated lightest heavy-light meson-meson thresholds.In particular, in the vector and pseudoscalar channels the thresholds depend on the calculated masses of the scalar and axialvector Qq states, which may differ from experiment by a couple of percent (cf.Sec.II A and Tab.III for details).

B. Charm spectrum
In the left panel of Fig. 6 we compare our results for the masses in Tables III and IV to the experimental spectrum in the charmonium region.We find that the cnnc and cssc states, depicted by blue and green boxes, respectively, lie quite close together in most channels.This closeness can be attributed to the plateau-like behaviour of the MECs in Figs.3-5 for small quark masses.The only exception is the pseudoscalar channel, where the MEC bends downwards at small quark masses.
Starting with the 1 ++ channel, we find that the cnnc ground state nicely agrees with the experimental state χ c1 (3872).On the other hand, the cssc ground state is too light compared to the χ c1 (4140), while the cnnc and cssc radial excitations are in the right mass region to be identified either with the χ c1 (4140) or χ c1 (4274).
In the 1 +− channel, the cnnc ground state is close to the Z c (3900) and the cssc ground state close to the X(4020) ± .Their first radial excitations might be candidates for the Z c (4200) + and Z c (4430), although their masses lie substantially above the respective thresholds and should thus be treated with caution.
In the 0 ++ channel, we find that the cnnc ground state agrees with the χ c0 (1P ), which in the literature, however, is considered as a cc ground state.In addition, also the cssc ground state appears in the same mass region.The excited cnnc state is in good agreement with the χ c0 (3915) and the excited cssc state matches very nicely with the recently observed X(3960) [41].
In the 1 −− vector channel we find that the cnnc ground state agrees with the ψ(4230), thus rendering it the lowestlying four-quark candidate in this channel with the dominant physical component being D D1 , followed by χ c0 ω and J/ψσ.The states below the ψ(4230) are not picked up by our analysis as they feature different decay channels.
As a caveat, we note that our calculated D 1 is substantially lighter than its experimental counterpart, which also lowers the D D1 threshold so that our state is far above the threshold whereas the experimental ψ(4230) is a shallow bound state.The corresponding cssc ground state is close to the ψ(4360), although an identification may be questionable as the dominant decays of the experimental ψ(4360) do not point towards ss components [34].The cssc excited state is however close to the ψ(4660), which because of its prominent decays to D s Ds1 (2536) and ψ(2S)ππ (ψ(2S)f 0 (980)) is assumed to be a hiddencharm, hidden-strange four-quark state.This leaves the cnnc excited state, which appears to be missing from the experimental spectrum.
The pseudoscalar 0 −+ channel features the physical components χ c0 η, η c f 0 and D * D1 .It should be kept in mind that in our present truncation the η only features an nn component and is thus mass-degenerate with the pion (cf.Table I).We therefore expect this component to be possibly too dominant in our current calculation as compared to a more complete approach.In the η c hadro-charmonium component we chose the f 0 (1370) as companion state, since the σ is itself a four-quark state [35,36] and too broad to act as companion.Finally, for the heavy-light meson components we did not consider the DD 0 combination, because the experimentally measured D 0 is again much too broad to form a molecular bound state as already argued in Refs.[9,42].Instead, we consider the combination D * D1 .
As a result, we obtain the cnnc ground state at a mass of about 3.37 GeV and a corresponding cssc ground state mass at 3.68 GeV.We also compared our result for the cnnc ground and excited state masses to the one using the masses for the η and χ c0 given in the PDG.We found that both masses increase about 200 MeV with respect to the masses quoted in Table III.However, using the PDG  III.Ground-state masses for the hidden-charm (cq qc) and hidden-bottom (bq qb ) states in GeV.For completeness we also display the binding energies EB with respect to the lightest (calculated) heavy-light meson-meson threshold in each channel except for the 0 −+ , where we take the χc0η in the charm and the η b f0 in the bottom region, which are the most relevant thresholds in the system; the "binding energies" for resonant particles above the threshold are shown in grey.The error given in the brackets is the combination of the extrapolation error and the error of the fit from Eq. (11).masses leads to a significantly higher threshold in the χ c0 η channel of about 3.96 GeV, rendering both states deeply bound.The masses for the corresponding excited cnnc and cssc states are 3.61 GeV and 4.00 GeV, respectively.Considering the lowest-lying heavy-light meson-meson S-wave thresholds for each channel we investigated, one can identify the following threshold hierarchy.The threshold in the scalar channel is the lightest, followed by the axialvectors, the vector and finally the pseudoscalar channel thresholds.Our states, even including other components than only heavy-light meson-meson pairings, follow this pattern except for the pseudoscalar channel, which is fully hadro-quarkonium dominated.We come back to this point in Section III D below when we discuss the internal structure of our states.

C. Bottom spectrum
Moving on to the bottomonium spectrum in the right panel of Fig. 6, we observe similar features as in the charmonium spectrum.Due to the plateau-like behaviour of the MECs for small quark masses, the cnnc and cssc states (green and blue boxes, respectively) are even closer to each other compared to the charmonium spectrum and overlap in most cases.
In the bottomonium spectrum there are only two experimentally well-established four-quark candidates, namely the Z b (10610) and Z b (10650) with quantum numbers 1 +− .These are potential members of an isospin triplet and therefore indistinguishable in our isospin symmetric framework.We find a bnn b ground state with a slightly lower mass than the experimental Z b (10610) which seems to match reasonably well.However, we also find a corresponding state with bss b flavour content close by, which has not yet been detected in experiment.
The predicted bottomonium partners of the charmonium-like 0 ++ , 1 ++ and 2 ++ four-quark states are in the literature referred to as W bJ states.In the 1 ++ channel, our bnn b ground state coincides with the experimental state χ b1 (3P ), which is a radial excitation of the χ b1 (1P ).We therefore predict a W b1 state with a mass of about m W b1 = 10.52(2)GeV.In the scalar channel, our bnn b ground state is close to the χ b0 (1P ), as was the case in the charmonium spectrum.If the pattern that the first excited state is more in line with the experimental four-quark candidates in this particular channel can be carried over from the charm to the bottomonium spectrum, we predict a W b0 with a mass of m W b0 = 10.38(2)GeV.We note that the scalar 0 ++ and axialvector 1 ++ W bJ masses using heavy-quark spin symmetry (HQSS) and effective field theories are predicted in a similar region [43], except that they are resonances above the respective B B and B B * thresholds while our states are below their thresholds.
Considering the spectrum in the vector channel, we find the ground state with quark content bnn b in the vicinity of three states, i.e., ψ(10753), ψ(10860) and ψ(11020).FIG. 6. Hidden-charm (left) and hidden-bottom spectrum (right) for the ground and first excited four-quark states compared to experiment [34].The coloured boxes are our results, where the height of the boxes stands for the error of the extracted masses.The gray and black boxes are the PDG masses (real parts of the pole positions) for conventional and exotic hadrons, respectively.The pale gray coloured states are not yet well established.
The latter two are considered to be radial excitations of the Υ(1S) often termed Υ(5S) and Υ(6S) in the literature.Therefore, despite the higher mass of our bnn b ground state an identification with the experimental ψ(10753) seems to be an option.
Finally, for the pseudoscalar channel we find the bnn b ground state at a mass of 9.9 GeV being slightly heavier than our obtained bnn b scalar state mass.Also here, substituting our pseudoscalar nn (η) mass with the η mass from the PDG yields a bnn b ground state which is about 200 MeV heavier.
The bss b and bcc b states in each channel do not currently have any experimental candidates.We can, however, compare our findings to predictions from the literature, especially regarding the bcc b states.These have been investigated using various methods such as augmented QCD sum rules [44], diquark-antidiquark models [45] and lattice-QCD inspired quark models [46].The results for the investigated channels with these methods lie mostly above 12 GeV and are very close or above the respective meson-meson thresholds.From the MECs in Figs.3-5 and the resulting binding energies in Tables III and IV one can clearly see that our states get more deeply bound if we increase the mass of the q q pair, i.e., when we go from bnn b to bss b and bcc b.Thus, as expected, we find the bcc b ground states (and even their first excited states) to be deeply bound in every channel.

D. Internal structure
One of the most interesting questions on four-quark states concerns their internal structure.In the preceding works using functional methods [21,27,37] such information was extracted from the MECs only.The mass spectrum of four-quark states with quark content cq qc was calculated by keeping the mass of the cc pair fixed and varying the mass of the light current-quark pair q q from up/down to charm.By changing the physical components entering in the calculation, one can then observe how well the MEC for a single sub-cluster (or combinations of sub-clusters) agrees with the MEC of the full state.For example, for the χ c1 (3872) the D D * component alone agrees reasonably well with the full result across a wide range of current-quark masses, while the J/ψ ω cluster contributes marginally and the effect of the SA diquark cluster is almost negligible (cf.Fig. 2 in [37]).
Here we employ a different strategy to obtain information about the internal structure of four-quark states.To this end, we investigate the dressing functions f i (S 0 ) in Eq. ( 8) and in particular their norm contributions, similarly to Refs.[22,23,47] where the orbital angular momentum composition of baryons and the strengths of different internal diquark clusters were quantified along the same lines.We follow the canonical normalization procedure for BSEs [48,49] by contracting the amplitude in Eq. ( 6) with its charge conjugate Γ(µ) .Because the amplitude sums over the three dominant internal structures (cf.Fig. 2), the product Γ(µ) G (4) 0 Γ (µ) is a sum of nine terms as illustrated in Fig. 8.The diagonal terms represent the contributions coming from the three topologies Current-quark mass evolution of the norm contributions for the cq qc ground states in the 0(1 ++ ) channel.The colour coding is the same as in Fig. 8.
M 1 , M 2 and D, and the off-diagonal terms arise from the mixing of these topologies.
To make statements about the internal structure of the physical four-quark states with quark content cnnc and bnn b, we need to calculate the norm contributions with the on-shell Bethe-Salpeter amplitude for that quark configuration.However, for some channels it is not possible to use the on-shell amplitude directly due to the two-body thresholds.We can, however, calculate the norm contributions for a four-quark state with quark content Qq q Q with the on-shell BSAs at quark masses for q where the two-body thresholds do not affect the state, cf.Figs. 3, 4 and 5.This yields a quark mass evolution of the norm contributions, in analogy to the MECs, which is then extrapolated to the physical quark content.An example for the 0(1 ++ ) channel is shown in Fig. 7.One can see that the evolution follows a clear trend and does not change drastically when varying the quark masses.The opaque datapoints at small current-quark masses are the norm contributions not considered in the fit and correspond to the opaque points in Fig. 3.The curves for the other quantum numbers and the excited states behave in a similar way.Fig. 9 shows the results for the contributions of this correlation matrix elements to each state.The dressing functions f 1 , f 2 and f 3 correspond to the first, second and third physical component in Table II for each quantum number, e.g., D D, J/ψω and S c S c for the 0(0 ++ ) state.The plot on the left illustrates the arrangement and colour scheme.For the 0 −+ , 1 −− , 1 +− and 1 ++ channels we only show the contributions for the cnnc, bnn b and bcc b ground states, since they hardly change for the corresponding excited states or states with hidden strangeness (cssc and bss b).In the scalar channel we plot the results for the respective excited states as they are more in line with the exotic candidates.8. Graphical illustration of the norm contribution matrix.Each entry in the matrix is an overlap integral that contributes the normalization of the four-quark state shown in the denominator.The diagonal terms correspond to the norm contribution coming from the M1, M2 and D topologies and the off-diagonal terms arise from the mixing of the topologies.Note that the matrix is symmetric, so that the contributions with the same background colour are summed up.
Starting with the 1 ++ channel in the charmonium sector, we find that this state has an overwhelming D D * component (blue) which contributes about 88% to the state.The J/ψω component (red ) and its mixing with D D * (orange) are almost negligible, but these are still bigger than the S c A c diquark component (brown) and its mixing with D D * (green).This nicely reproduces the hierarchy found in [27,37].It also mirrors the experimentally known decays: the dominant hadronic decay channel for the χ c1 (3872) is D D * with ∼ 86% when combining the D 0 D0 π 0 and D 0 D * 0 channels, followed by J/ψω with ∼ 8% [34].Furthermore, as the χ c1 (3872) is very close to the D D * threshold, a strong D D * component in its wave function is expected.The same behaviour is also found for the cssc ground and excited state, which are dominated by the D s D * s component.In the bottom sector we find the (would-be W b1 ) bnn b state to be dominated by the B B * component with about 77%; here the other components are still all below 10% but the mixing of B B * − Υω and B B * − S b A b becomes more prominent.
In the 1 +− channel we find the J/ψπ component to be dominant (55%), however with a substantial D D * admixture (29%) and a D D * − J/ψπ mixing component (14%).In the literature the internal structure of the Z c (3900) is debated as its mass is close to, but still above, the D D * threshold.Furthermore, the D D * decay channel is preferred over J/ψπ by a factor ∼ 6 [50].This led to the conclusion that the J/ψπ component might be suppressed and the Z c (3900) could be explained as a D D * molecule.The HAL-QCD lattice collaboration studied the internal structure of the the Z c (3900) using a D D * − J/ψπ coupled-channel analysis [51,52].They found a strong D D * − J/ψπ mixing potential and thus evidence against the state having only a meson-meson or hadrocharmonium structure.Instead, they concluded that this strong potential leads to a formation of the Z c (3900) as a threshold cusp.An analysis of the experimental data with effective field theories [53,54] [34] and the closeness to the B B * threshold [9].
In the 0 ++ channel we find a dominant D D component for the cnnc state which contributes about 89% to the state.The mixing components of D D − J/ψω and D D − S c S c amount to about 9%, while the remaining contributions are negligible.This again nicely confirms the findings of [21,27] where this state was found to be predominantly D D. The corresponding cssc state has a contribution of D s Ds of about 89% which fits with the observed decay channel of the X(3960) [41].In the bottom region, we observe again a B B dominance of 80% with the mixings B B − Υω and B B − S b S b becoming more prominent with a total contribution of 16%.
Concerning the cnnc state in the 1 −− vector channel, we find that it is almost exclusively dominated by the D D1 component (93%) with a small contribution coming from the χ c0 ω (4%) and D D1 −χ c0 ω mixing (3%).This is in line with Refs.[55,56] concluding that a description of the ψ(4230) as a D 1 D molecular state agrees with the experimental data.As the D D1 threshold is also the lowest S-wave meson-meson threshold in the system [6,9,57], the closeness of the state to that threshold also points to a strong meson-molecule component in the wave function.The internal structure of the ψ(4660), which we identified with our cssc excited state, is a little more elaborate.Here, motivated by the observed decays, there are claims to describe this state as a hadro-charmonium ψ(2S)f 0 (980) state [58], a D ( * ) s Ds1 (2536) meson-molecule [59], or a Pwave tetraquark (dq − dq) state [60].In our analysis we find the D s Ds1 component to be dominant with 92%.The picture is quite similar in the bottom region: All states are dominated by the respective B B1 component with about 83%, followed by a contribution of about 10% coming from the pure χ b0 ω cluster.
Finally, we turn to the 0 −+ pseudoscalar channel.For the cnnc state we find that with a dominant χ c0 η component of 99% our obtained cq qc states are exclusively hadro-charmonium with no mixing components.As discussed, this picture might change when we include the strange flavour components in the η together with the octet-singlet mixing and therefore needs to be explored again in a more complete approach in the future.In the bottom sector the dominant substructure is still hadro-quarkonium, but the state is now almost exclusively dominated by the η b f 0 (1370) hadro-bottomonium component with about 95%, augmented by an almost negligible χ b0 η contribution of 5%.Presumably this result will not change in a more complete approach given that a heavier η will decrease the importance of this component rather than increase it.When going to the bcc b state, the dominant component is still η b χ c0 with 96%, with the B * B 1 and the χ b0 η c components contributing about 2% each.

IV. SUMMARY AND CONCLUSIONS
In this work we determined the spectrum and internal composition of heavy-light four-quark states with hidden flavour in the charm and bottom energy region.Using rainbow-ladder two-body interactions between quarks and (anti-)quarks, we solved the four-body Faddeev-Yakubovsky equation in a fully covariant framework and obtained spectra for states with quantum numbers J P C = 0 ++ , 0 −+ , 1 −− , 1 +− and 1 ++ .Our wave functions routinely incorporate contributions from internal heavy-light meson-meson, hadro-quarkonium and diquarkantiquark contributions.Whereas we find the latter ones to be subleading in most cases, it turns out that states with different quantum numbers correspond to different internal contributions, which in some cases also depend on the (hidden) flavour of the states.For all states with CP = +1 (i.e.J P C = 0 ++ , 1 −− and 1 ++ ) we find highly dominant heavy-light meson contributions, which almost exclusively determine the masses of the ground and excited states in the charm and bottom energy region.For axialvector states with J P C = 1 +− , however, we find a much more complicated picture.Our hidden-charm state corresponding to the experimental Z c (3900) comprises both heavy-light meson but also hadro-charmonium components in qualitative agreement with results from other approaches [51][52][53][54].The corresponding state in the bottom energy region is even dominated by the hadrobottomonium configuration.The same is observed for our pseudoscalar states regardless of flavour.Thus the most important message from our study is: the internal composition of XYZ states is by no means uniform but varies between different quantum numbers and flavours.It is certainly interesting to extend our findings to the open flavour case.Corresponding work is in progress and will be reported elsewhere.

Appendix A: Construction of amplitudes
In this section we describe the construction of the Bethe-Salpeter amplitudes for the 0(1 −− ) vector and 0(0 −+ ) pseudoscalar four-quark states based on their dominant physical clusters.We refer to the supplemental material of Ref. [21] for the analogous construction for the quantum numbers J P C = 0 ++ , 1 +− and 1 ++ .
To begin with, we consider the quantum numbers of the meson-meson and diquark-antidiquark components that can appear in a four-quark state with the desired quantum numbers.These components can then be assigned to one of the interaction topologies M 1 and M 2 (meson-meson) or D (diquark-antidiquark).We then compare with the PDG [34] for the experimentally dominant or realized decays of the four-quark state under consideration.If there is no experimental evidence regarding the dominant decays, we take the lowest-lying thresholds of the system as the dominant contributions spanning the physical basis.For the vector states this procedure is straightforward as there are well-established exotic candidates with vector quantum numbers such as the ψ(4230).At present, however, there is no evidence of a pseudoscalar four-quark candidate, so we have to determine its physical clusters based solely on the lowest-lying thresholds of the possible internal components.
As a first approximation, we neglect the diquark components in the basis for the vector and pseudoscalar fourquark states.The reasons for this are twofold.First, the possible diquark clusters feature either S-wave pairings of scalar/axialvector diquarks with pseudoscalar or vector diquarks (which are not only heavier than their scalar and axialvector counterparts, but also unreliable in a rainbow-ladder truncation), or P -wave scalar and axialvector diquark cluster pairings (which include higher orbital angular momentum and are therefore suppressed).The notion of S and P wave here refers to the combination of Dirac tensors needed to construct a basis element of a four-quark state with vector or pseudoscalar quantum numbers.The P -wave combination would induce angular momentum between the two diquarks, which is however not possible with the approximations explained in Sec.II C. Secondly, the diquark clusters always form the highest threshold in the system and were previously found to have a negligible influence on the mass of hidden-flavour four-quark states [21,[36][37][38].The resulting physical components chosen for the vector and pseudoscalar four-quark states are then only of meson-meson type and given in Table II.
To construct the basis describing the pseudoscalar and vector states, we write down the leading Dirac-colour tensors for I(0 −+ ) and I(1 −− ) in analogy to the supple- works as follows: Having obtained a set of eigenvalues A r = {λ(P 2 i )} r i=1 , we extrapolate the behaviour of all r eigenvalues as a function of P 2 to the value where the condition λ(P 2 i = −M 2 i ) = 1 is fulfilled.This gives us a base estimate of the mass, i.e., M base , and the first element in a set of extrapolated masses denoted by B. Next we choose a random subset of eigenvalues A m ⊂ A r , with m ∈ {r − 1, r − 2, r − 3, r − 4, r − 5, r − 6} and extrapolate the eigenvalues chosen in A m to fulfil the condition again.This procedure is repeated about 300 times for each m.The corresponding extrapolation results are added to the set B if they lie in a 5% region around the M base value.To obtain the masses shown in Figs.3-5 we average the values in the set B and the error is given by the standard deviation of the values in B. A nice example of the distribution of eigenvalues in B is shown in the histogram in Figure 10.This method has also been applied to states where the condition λ i = 1 can be read off from the eigenvalue curve directly.The corresponding histogram plot of eigenvalues then shows an extreme dense clustering around the directly determined value, which is the expected behaviour.
Throughout this paper we use the abbreviation n = u/d when referring to the light up and down quarks.The u/d-quark mass is fixed by m π to m n = 3.7 MeV at a renormalization point µ = 19 GeV in a MOM scheme.The strange and charm quark masses are chosen as m s = 85 MeV and m c = 795 MeV (in the same renormalisation scheme) such that the sums m Ds + m D *

FIG. 2 .
FIG. 2. Graphical representation of the Bethe-Salpeter amplitude in the physical basis.The diagrams on the r.h.s.represent the direct product of the internal clusters spanning the physical basis.The first diagram shows the meson-meson configuration (M1), with the heavy-light meson clusters depicted as blue half-circles.The second diagram is the hadro-quarkonium contribution (M2), with individual components shown as violet half-circles, and the last diagram is the diquark-antidiquark configuration (D), where each diquark is depicted by an orange half-circle.
of the Bethe-Salpeter amplitude for cnnc and bnn b configurations, where n stands for light u/d quarks.Scalar and axialvector diquarks are denoted by S c/b FIG.4.Current-quark mass evolution of the cq qc and bq qb ground states in the 0(0 ++ ) channel; see Fig.3for details.
bss b excited bss b ground bnn b excited bnn b ground bcc b excited bcc b ground Exp.(exotic) Exp.(conv.)

FIG. 9 .
FIG.9.Norm contributions for the cnnc, bnn b and bcc b states.In the scalar channel we show the results for the excited states and in all other channels those for the ground states.The f0 here corresponds to the f0(1370).The bars sum up to 100%.
suggests that the mesonmolecule D D * component is at least equally important than non-molecular structures.Interestingly, the contributions in the 1 +− channel change in the bottom sector.Here we find a very strong Υπ component (86%) and a Υπ − B B * mixing of about 10% for the bnn b state, with the other correlations being negligible.For the bcc b state the Υη c component is weaker (∼ 75%) and the B c B * c and Υη c −B c B * c mixing contribute about 10% to the state.This is somewhat at odds with the common picture in the literature, where the Z b (10610) is considered as a B B * molecule because of the dominant decay Z b (10610) → B + B * 0 + B * + B0 with about 86%

FIG. 10 .
FIG. 10.Example distribution of extrapolated eigenvalues for the I(J P C ) = 0(1 ++ ) state with quark configuration cq qc and mq = 2750 MeV.The extrapolation results are centred around a value of Mcqqc ≈ 9.4 GeV with a spread that resembles a Gaussian distribution.
FIG.11.Current-mass evolution of the cq qc (crosses) and bq qb first radial excited states (dots) for the investigated quantum numbers in this work; see Fig.3for details.