Color halo scenario of charmonium-like hybrids

The internal structures of $J^{PC}=1^{--}, (0,1,2)^{-+}$ charmonium-like hybrids are investigated under lattice QCD in the quenched approximation. We define the Bethe-Salpeter wave function $\Phi_n(r)$ in the Coulomb gauge as the matrix element of a spatially extended hybrid-like operator $\bar{c}{c}g$ between the vacuum and $n$-th state for each $J^{PC}$ with $r$ being the spatial separation between a localized $\bar{c}c$ component and the chromomagnetic strength tensor. These wave functions exhibit some similarities for states with the aforementioned different quantum numbers, and their $r$-behaviors (no node for the ground states and one node for the first excited states) imply that $r$ can be a meaningful dynamical variable for these states. Additionally, the mass splittings of the ground states and first excited states of charmonium-like hybrids in these channels are obtained for the first time to be approximately 1.2-1.4 GeV. These results do not support the flux-tube description of heavy-quarkonium-like hybrids in the Born-Oppenheimer approximation. In contrast, a charmonium-like hybrid can be viewed as a"color halo"charmonium for which a relatively localized color octet $\bar{c}c$ is surrounded by gluonic degrees of freedom, which can readily decay into a charmonium state along with one or more light hadrons. The color halo picture is compatible with the decay properties of $Y(4260)$ and suggests LHCb and BelleII to search for $(0,1,2)^{-+}$ charmonium-like hybrids in $\chi_{c0,1,2}\eta$ and $J/\psi \omega (\phi)$ final states.


I. INTRODUCTION
In the constituent quark model, mesons are interpreted as quark-antiquark (qq) bound states. However, since gluons are also fundamental degrees of freedom in quantum chromodynamics (QCD), they are expected to build hadrons either by themselves, such as glueballs, or make create hybrids with quarks. Glueballs and hybrids are exotic hadron states that have been searched for a long time through experiments. For the hybrid mesons (denoted by QQg) involving a heavy quark-antiquark pair QQ, an interesting phenomenological description is the flux-tube picture in the leading Born-Oppenheimer approximation [1,2], in which the gluonic excitations are considered as fast degrees of freedom, which concentrate along the QQ axis and result in an instantaneous confining potential that obeys the cylindrical symmetry along the QQ axis and the refection symmetry with respect to the QQ midpoint. These gluonic excitations provide effective potentials between the QQ pair, which are frequently labeled as Λ η in analogy with diatomic molecules, where Λ = 0, 1, 2, · · · is the magnitude of the angular momentum of the gluon state projecting to the axis and is conventionally denoted by Σ, Π, ∆, respectively; η = ± is the CP quantum number of the gluon state, and = ± is the reflection quantum number of the gluon state with respect to a plane perpendicular to the axis at the midpoint of the axis. Based on this description, many phenomenological studies have been conducted on QQg hybrids through gluon-excitation pictures [3,4] or by solving the non-relativistic Schrödinger equations of QQ systems with this type of potential [5,6].
QQg hybrids, particularly the charmonium-like hybrids ccg, have also been extensively investigated in lattice QCD [7][8][9][10][11][12][13][14][15]. A recent lattice calculation [16] demonstrated the existence of a {1 −− , (0, 1, 2) −+ } charmoniumlike supermultiplet with nearly degenerate masses of approximately 4. 2-4.4 GeV, which overlaps strongly to ccg type operators. These states may have similar internal dynamics, while the different quantum numbers are due to the different couplings of the spin states of thecc and chromomagnetic excitation (J P = 1 + ). In the Born-Oppenheimer potential picture, this supermultiplet can be assigned to Π + u (1P ) states since they are the lowest in mass. The properties of this supermultiplet have practical interest for current experimental and theoretical studies of XY Z particles (see Ref. [17] for a review). Among the XY Z particles, Y (4260) (or ψ(4230) as named by PDG 2018 [18]) can have the possible assignment of a vector charmonium-like hybrid [19] and can be a member of the aforementioned supermultiplet, owing to its strange production and decay properties, as well as its mass adjacent to that of the 1 −+ charmonium-like state predicted by previous phenomenological and lattice QCD studies. In addition to its mass, a quenched lattice calculation also predicted the leptonic decay width of the vector charmonium-like hybrid to be approximately smaller than 40 eV [20], which explains to some extent the absence of Y (4260) (if a hybrid) in the R-value scan of e + e − annihilation processes; this was compatible with the estimate from its isospin symmetric decays [21,22]. Therefore, a joint study of this multiplet is necessary to understand the experimental observations relevant to Y (4260) and predict the properties of other members to be searched, particularly the state of the exotic 1 −+ quantum number.
In this work, we investigated the internal structure of the {1 −− , (0, 1, 2) −+ } charmonium-like hybrid supermultiplet under lattice QCD. As an exploratory study, we adopt the quenched approximation, in which at least the 1 −+ hybrid is well-defined. By constructing spatially extended operators for hybrids, we extracted their Bethe-Salpeter (BS) wave functions through the corresponding correlation functions calculated in the Coulomb gauge.

II. FORMALISM
We generated gauge configurations on two anisotropic lattices using the tadpole improved gauge action [23,24]. The aspect ratio is set to be ξ = a s /a t = 5, where a s and a t are the spatial and temporal lattice spacings, respectively. The much finer lattice in the temporal direction enables us to address heavy particles on relatively coarse lattices. The configuration parameters are listed in Table I, where the values of a s are determined from r −1 0 = 410 (20) MeV. For the charm quark, we use the tadpole improved clover action for anisotropic lattices [25], and the bare charm quark mass is determined by the physical mass of J/ψ, m J/ψ = 3.097 GeV. As will be addressed in the following sections, we use spatially extended operators to calculate the relevant correlation functions; therefore, the configurations are first fixed to the Coulomb gauge through the standard gauge fixing procedure [26] in lattice QCD studies before the quark propagators are computed.
In this study, we extracted the spectrum and the BS wave functions of charmonium-like ccg hybrids of quantum numbers {1 −− , (0, 1, 2) −+ }. For a particular quantum number, the state-of-art technique in extracting excited states in lattice QCD is the variational method (VM) based on an operator set {O α (t), α = 1, 2, . . . , M } of M different operators. The procedure of VM is outlined as follows: First, the correlation matrix C αβ (t) = O α (t)O † β (0) is calculated on a gauge ensemble. Second, the generalized eigenvalue problem can be solved to derive the eigenvector v (m) α corresponding to the m-th eigenvalue λ (m) (t). We can obtain an α O α (t) that couples mostly to the m-th state of the given quantum number.
The key aspect of VM is that the operators O α in the operator set are as different as possible from each other in the sense that they couple to each state in the spectrum differently, such that the linear combination α O α (t) produces almost distinct operators for different m values. Therefore, smeared quark and gluon fields are conventionally adopted in constructing operators O α using the conjecture that the operators with different smearing sizes satisfy the requirement. VM has been extensively and successfully applied in the precise determination of energy levels in lattice QCD studies. However, VM is not suitable for the derivation of hadron decay constants and form factors, which are usually defined through the matrix elements of local operators between physical states.
The BS wave function is defined by the matrix element of a spatially extended operator O(r) between the vacuum and a hadron state, i.e., Φ n (r) ∼ Ω|O(r)|n , where r reflects the internal spatial coordinates of a hadron. The effect of O(r) on the vacuum state |Ω is which implies that the BS wave function Φ n (r) actually controls the coupling of the operator O(r) to the state |n . Therefore, if the r behaviors of Φ n (r) are different for different state |n , then O(r) can be considered a different operator for different r values. On the other hand, the correlation function of O(r, t) and some other source operator O S (0) has the following spectral expression where m n is the mass of the n-th state and the rindependent factors are absorbed into Φ n (r) in the last equality. Since the m n spectrum is common for different r values, if N r correlation functions C(r, t) with different r values are calculated, then in proper time windows t ∈ [t min , t max ], the m n and Φ n (r) values for n = 1, 2, . . . , N m can be determined by fitting these correlation functions simultaneously with the function form in Eq. (3) involving N m mass terms. Specifically, if the time window [t min , t max ] is uniform for all values of r and there are N t time points in this window, there will be N r ×N t data points for the N m +N r ×N m parameters to be fitted. Thus, the number of the degrees of freedom can be sufficiently large if N t is significantly larger than N m . We adopted the above strategy to study the spectrum and BS wave functions of {1 −− , (0, 1, 2) −+ } charmoniumlike hybrids. The operator prototypes for these states are selected to be the traditionalcΓc • B type:  R  A1  A2  E  T1  T2  d(R)  1  1  2  3  3  J  0, 4  3  2, 4, 5  1, 3, 4  2, 3, 4 O k where a, b = 1, 2, 3 are color indices, B k is the chromomagnetic field strength, and the summation over x guarantees the operator to couple to a state in its rest frame. We can easily observe that the quark-bilinear operators,cγ i c andcγ 5 c, are a spin triplet and singlet, respectively, in the non-relativistic approximation. Note that the spatial symmetry of the lattice is described by the octahedral group, O, whose irreducible representations are R = A 1 , A 2 , E, T 1 and T 2 , whose dimensions d(R) are 1, 1, 2, 3, 3 respectively. Group O is a subgroup of SU (2) such that the subduced representation of SU (2) with respect to O is generally reducible. Table II  of the lattice symmetry group, respectively, they have the same forms as the 0 −+ , 1 −+ , 2 −+ and 1 −− operators in the continuum limit. Therefore, we do not distinguish them from each other under the assumption that the lattice discretization effects are small. To investigate the inner structure of the charmonium-like states through BS wave functions, we introduce two types of spatially extended operators based on the operator prototypes in Eq.(4). Type-I operators are constructed by splitting thecΓc component from the chromomagnetic field strength operator B by a spatial separation r: where D(r) is the number of r values that satisfy |r| = r, Γ R is the gamma matrix appearing in the operator prototype O R in Eq.(4) with R = (0, 1, 2) −+ and 1 −− , and the symbol "•" indicates the proper color summation and specific combination of the spatial indexes. Type-II operators are obtained by assuming a spatial separation betweenc and c with B being at the same point as c: (6) where "c.c." is the complex conjugate of the first term that is required to guarantee the correct parity and charge conjugate quantum numbers.
These two types of operators are not gauge invariant; therefore, the correlation functions involving these operators should be calculated in a fixed gauge. In practice, we work in the Coulomb gauge and calculate the following correlation functions for each quantum number R: where O (W ) R (τ ) is a wall-source operator defined on the time-slice τ for R as In the data analysis stage, these correlation functions are parameterized to the function form of Eq.(3), based on which masses m n and BS wave functions Φ n (r) are derived simultaneously through a correlated joint fit to the lattice data of C R (r, t) at different values of r and t.

III. RESULTS USING TYPE-I OPERATORS
In this section, we present the results of the spectrum and the BS wave functions of {1 −− , (0, 1, 2) −+ } states using the type-I operators O (I) R (r, t), where r is the spatial separation between the cc and B components. In each channel, the Coulomb wall-source correlation functions C R (r, t) are modeled as where τ is the source time slice and N τ = 10 is the number of sources on each configuration, and the spatial indexes of operators are summed implicitly in the quenched approximation. Based on the aforementioned numerical strategy, we performed the data fitting using the model in Eq. (9) to C R (r, t) functions with r ∈ [0, r max ] (n = 3, r max = 0.89 fm for β = 2.4 and n = 3, r max = 0.74 fm for β = 2.8). In the simultaneous fitting procedure, the correlations of data points at different r and t values were considered by calculating their jackknife covariance matrix. The key point is that the sink operator with different r values couples differently to different states; thus, it has different spectral weights Φ (I) n (r). We fixed the upper limit of the fitting window (t max /a t = 19 for both β = 2.4 and β = 2.8) and let the lower bound t min vary. The fit results of m n at different t min values on the two lattices are listed in Table III, where the masses are converted to the values in physical units using the lattice spacings in Table I; the table also lists the χ 2 values per degree of freedom (χ 2 /dof). As shown, the fitted masses are very stable and insensitive to t min . Additionally, the masses of the lowest two states exhibit slight finite a s effects. Figure 1 shows the fitted masses m n (in physical units) of the lowest three states versus t min , where the colored bands are the averaged mass values weighted with the reciprocals of the error squared at different values of t min .
The Φ (I) n (r) functions can be derived simultaneously for the three lowest states. Table IV shows the results at t min /a t = 2 for β = 2.4 and t min /a t = 5 for β = 2.8, where the fitted parameters have smaller statistical errors and acceptable values of χ 2 /dof (see Table III). After the normalization using Φ  belongs to the 1D multiplet {2 −+ , (1, 2, 3) −− } whose 1 −− member is ψ(3770) and 2 −− member ψ 2 (3823) was observed by Belle [27] and BESIII [28]. Recently, LHCb also observed a candidate ψ 3 (3842) for the 3 −− charmonium [29]. Their masses comply with the prediction of lattice QCD studies [16]. In addition to η c2 , lattice QCD studies find a higher 2 −+ state with a mass around 4.4 GeV, which couples almost exclusively to ccg operators [16,30]. This is also the case in this paper. The 2 −+ correlation functions C R (r, t) using the operators O k (4), and the data analysis is the same as for 1 −+ . The mass spectrum and BS wave functions ( Fig. 2 and Table V ) are very similar to those of 1 −+ states. n (r), n = 1, 2, 3 of 1 −+ states at tmin/at = 2 for β = 2.4 and tmin/at = 5 for β = 2.8.   spin triplet, respectively. In contrast, recalling that thecc component of the hybrid operators defined in Eq. (4) for 0 −+ and 1 −− are in spin triplet and spin singlet, respectively, we expect that the couplings of these operators to conventionalcc states will be suppressed to an extent owing to the spin-flipping of charm quarks. Despite this kind of suppression, we observe that the lowest lying conventionalcc states contribute significantly to the twopoint functions in Eq. (7). We perform three-mass-term and four-mass-term fits to the correlation functions of 0 −+ and 1 −− channels and address the fit results in the following.
For the three-mass-term fits in both channels, we fix the upper bound t max of the time window [t min , t max ] and vary the t min to monitor the stability of the fit. For β = 2.4, t max is set to t max /a t = 24, and t min /a t varies from 9 to 13. For β = 2.8 with t max being fixed at t max /a t = 39, and t min /a t varying from 15 to 19. The fitted masses m n and Φ (I) n (0) (n = 1, 2, 3) are listed in Table VI, where the minimal-χ 2 per degree of freedom   (χ 2 /dof) values are also given for all the t min . The masses are converted into values with physical units using the a s in Table I. We observe that the fits are all acceptable and stable in these fit ranges with reasonable χ 2 /dof values. The masses of the lowest two states, i.e., m 1 and m 2 , are in good agreement with those of 1S charmonia (η c and J/ψ) and 2S charmonia (η c (2S) and ψ(2S)), respectively, while m 3 is close to the ground state masses in 1 −+ and 2 −+ channels. In contrast, the spectral weight magnitude of the third state, namely Φ (I) 3 (0), is seemingly larger than those of the lowest two states. If the third state is dominated by the charmonium-like hybrid, this type of differences in the spectral weight magnitudes complies with our previous expectation that the couplings of hybrid-like operators to conventional charmonia are suppressed due to the spin-flipping effect. Here, we provide a qualitative interpretation of the negative sign of Φ (I) 2 (0). As introduced in Sec. II, the correlation functions C R (r, t) are calculated using the wall-source operator in Eq. (8), which can be viewed as X,sc (X)Γ R c (X + s) with c being a "dressed" charm quark field Bc, such that its matrix between the vacuum and a charmonium(-like) state H in its rest frame can be qualitatively expressed as where V 3 is the spatial volume and φ H (s) is the BS wave function of H with respect to the distance s between the charm and anti-charm quark. If H is a 2S charmonium state, we expect that φ H (s) has a radial node (this is actually observed;, see below) such that the integral in the above equation will result in a negative sign (note that the integrand is enlarged by the factor s 2 for a large s). We now discuss the results of the four-mass-term fits. To incorporate the data points in the small t region, we add the fourth mass term in the fit model. Table VII shows the fit results. We observe that the four-massterm model closely fits the data for t min /a t starting even from t min /a t = 1 for β = 2.4 and t min /a t = 2 for β = 2.8. The ground state masses are compatible with the 1S charmonia (η c for 0 −+ and J/ψ for 1 −− ) but are slightly lower than the results of the three-mass-term fits. The masses of the second lowest states (m 2 ) are approximately 4.4 GeV and significantly higher than the expected 2S charmonium masses. The masses of the third and the fourth states (m 3 and m 4 ) are close to those of the second and third states in the 1 −+ and 2 −+ channels (see Table III and V). For the BS wave functions at the origin r = 0, Φ (I) 1 (0) is very stable with respect to the varying of t min /a t but slightly smaller than the threemass-term fit results (this may be due to the smaller fitted masses m 1 ) whereas Φ (I) 2 (0) changes dramatically. Because there are a few conventional charmonium states at approximately 4 GeV, e.g., ψ(3686), ψ(4040), ψ(4415) etc. and a would-be charmonium-like hybrid state in the 1 −− channel, the second state may be an admixture of these states, such that its fitted spectral weight Φ 2 (r) is very sensitive to the fit window. In contrast to Φ   Until now, we have provided a detailed description of the fitting procedures and fitted results of the correlation functions involving the cc − g operators in the 1 −− and (0, 1, 2) −+ channels. The results of the 1 −+ and 2 −+ are solid and have slight ambiguities. The situation for 0 −+ and 1 −− channels is a little more complicated since both conventional charmonia and possible charmoniumlike hybrids contribute to correlation functions. In the three-mass-term fits for both channels, at larger t min values of the fit window, the ground states are stable with masses consistent with those of J/ψ and η c , the second states have masses of approximately 3.7 GeV, while the third states have masses of roughly 4.5 GeV with larger spectral weights. We consider the lowest two states to be the 1S and 2S conventional charmonia, and tentatively assign the third state in each channel to be the corresponding ground state charmonium-like hybrid. In the four-mass-term fits with smaller t min values, for each channel, while the lowest state does not change much, the mass of the second state varies slightly, and the corresponding wave function changes drastically. This can be attributed to the interplay of conventional nS charmonia and the ground state hybrid. The third state is relatively stable with varying t min and has a mass close to those of the second states in (1, 2) −+ channels; therefore, we consider it to be the first excited hybrid state in the corresponding channel. The r-dependence of the BS wave functions reinforces the above assignments: the wave functions of the three states in the three-massterm fit have no nodes in the r-direction, and those of the third and fourth states in the four-mass-term fit have one and two nodes, respectively. The fitted mass spectrum is shown in Table VIII, in which the values are averages of the two lattices. After singling out the states that correspond to the conventionalcc states in 0 −+ and 1 −− channels, the masses of the states whose BS wave functions have the same number of nodes are listed in the same row for all the four channels. We observe that the masses of the 0-node and 1-node states are nearly degenerate (we do not consider the masses of the 2-node states to be significant since they may have substantial contamination from higher states). The BS wave functions of the 0-node states (except for the states corresponding the conventional charmonia in 0 −+ and 1 −− channels) and 1-node states in the four channels are plotted in  3 (r) of (1, 2) −+ states (see Fig. 1 and Fig. 2).

IV. RESULTS FROM TYPE-II OPERATORS
As a self-consistent check, another type of BS wave function Φ The data fitting strategy is the same as that for type-I operator, and will not be repeated here. We directly present the results (this part of the calculation was carried out only on the β = 2.4 lattice).
A. 0 −+ and 1 −− states Let us begin with the results of the 0 −+ and 1 −− channels. The results of m n and Φ (II) n from the two-massterm, three-mass-term and four-mass-term fits in the 0 −+ channel, for which the masses are converted into values with physical units, are presented; the χ 2 /dof of each fit is also presented to indicate the fit quality. The results in the 1 −− channel are similar, which can be checked by comparing Table IX and Table X; therefore, we omit these to save space.
When t min /a t > 15, the C R (r, t) functions can be fitted using two mass terms with reasonable χ 2 /dof values, as shown in the first five rows in Table IX  VIII. Mass spectrum of the 1 −− and (0, 1, 2) −+ states. The lowest three masses in the 1 −− and 0 −+ channels are from the three-mass-term fits. The other two masses in the two channels are the masses of the third and fourth states, respectively, from the four-mass-term fits. Except for masses of the lowest two states in 1 −− and 0 −+ channels which correspond to conventional charmonia, other masses are arranged in each row by the number of nodes (#node) of the BS wave functions of related states. The values are averages of the two lattices and are converted into the physical units using the lattice spacings listed in Tab. I. (r) has no node in the r-direction, whereas Φ (II) 2 (r) has one radial node. This type of r-behaviors is in qualitative agreement with the expectation of the non-relativistic quark model and can be understood as follows. First, we assume the non-relativistic approximation is reasonable for charm quark systems to a certain extent. Since operator O (II) R (r, t) has a spatial structure in that the charm fieldc is spatially separated from block cB by distance r, and the BS wave function is defined by Φ n (r) can be interpreted as the probability amplitude of annihilating an anti-charm quark at the origin byc and annihilating a charm quark at r by cB in the state |n . In other words, block cB can be considered to be a dressed charm quark field belong to the fundamental representation of the color SU (3) group. Thus Φ (II) n (r) qualitatively reflects the non-relativistic wave functions of the corresponding cc state and thereby has the expected nodal structure.
When t min decreases to the range 7a t ≤ t min ≤ 11a t , the third mass term is required in the fit function to describe the lattice data of C R (r, t). The results of the three-mass-term fits are tabulated in the middle part of Table IX, Table VI). Note that the values of Φ (r) has no radial nodes. In this sense, the lowest two states can be identified to be η c (1S) and η c (2S), while the third state is significantly higher and can not be the purely η c (3S) state but is likely dominated by the would-be lowest hybrid charmonium state.
If t min decreases further to t min /a t = 2, 3, 4, then the correlation functions in 0 −+ and 1 −− channels can be fitted by four mass terms with χ 2 /dof being 2.27, 0.99 and 1.06 for 0 −+ , and 1.91, 1.04, and 1.05 for 1 −− . The TABLE IX. Fitted masses of 0 −+ states at different tmin values through the two-mass-term (top part), the three-mass-term ( middle part) and the four-mass-term fits to the correlation functions involving type-II operators. The mass values are converted into physical units using the lattice spacings listed in Tab. I. The χ 2 /dof of each fit is provided to indicate the fitting quality. Φ (II) n (r) values at r = 0 are also listed.   tmin/at  Table IX. Although the mass m 1 of the lowest state is compatible with (actually slightly smaller than) the η c (1S) state, m 2 and m 3 are clearly different from those of the two-mass-term and three-mass-term fits. Φ (II) n (r) functions obtained at t min /a t = 3 are shown in Fig.8. The r-behavior of these states are very strange and have no radial nodes; therefore, no physical information can be inferred yet. The reason for this observation may be that, in the small t range many states significantly contribute to the correlation functions including the conventional and possible charmonium-like hybrids, such that the fitted second and third state might be the admixture of them. This situation is similar to that of the four-mass-term fits for the type-I operator.

B. 1 −+ and 2 −+ states
Since the quantum number 1 −+ is exotic, the spectrum of this channel is much simpler than those of 0 −+ and 1 −+ channels. As mentioned earlier, the hybridlike operators might couple with hybrid states predominantly in the 2 −+ channel although it is permitted for conventional D-wave charmonia. Therefore, similar to the case of type-I operators, the C R (r, t) functions can be well described by the function form of three mass terms in the time window beginning from the very early time slices. The fitted results for different t min /a t values (t max /a t is fixed at 19) in these two channels are tabulated in Table XI. The χ 2 /dof values of these fits are approximately one or smaller and illustrate the goodness of the fits. We observe that the masses of the three states in each channel are in very close agreement with those of type-I operators (see Table III and Table V). This is exactly expected since the spectrum is independent of the operators for the same quantum number. The Φ (II) n (r) functions at r = 0 are approximately twice as large as those of type-I operators owing to the factor of two in the definition of the type-II operator at r = 0.
The BS wave functions Φ (II) n (r) can be derived precisely and the results at t min /a t = 3 are plotted in Fig. 9 (the upper panel for 1 −+ and the lower one for 2 −+ ). The r-behaviors are strikingly different from those of the type-I operator and there are no nodes in the r direction at all. A tentative interpretation of this difference is that, if a charmonium-like hybrid can be considered to be ā c − c − g three-body system, then its internal motion can be described through the Jacobi variables (ρ, λ), where ρ is the relative displacement betweenc and c, and λ is the displacement between the gluonic degree of freedom from the center of mass of cc. Thus the BS wave function from type-I operators reflect the internal motion with respect to λ, and Φ (II) n (r) signal the projection of the full wave function φ(ρ, λ) onto the r-direction. In contrast to the BS wave function derived from the type-I operators, the r-behavior of Φ (II) n (r) does not exhibit the typical feature of radial excitations and the separation r ofc and cB components is a less meaningful variable for the pattern of excited hybrid states.

V. DISCUSSION
In the previous sections, we presented the spectrum and BS wave functions for 1 −− , (0, 1, 2) −+ states by fitting the correlation functions of type-I and type-II hybrid-like operators. The results are summarized as follows: • The spectrum is independent of the operators, as expected. If we exclude the lowest two states in 0 −+ and 1 −− channels, which correspond to the 1S and 2S conventional charmonia, all the four channels contain states that are nearly degenerate in mass of approximately 4.4 GeV and around 5.6 GeV, respectively, which can be tentatively assigned to be the ground and first excited charmonium-like hybrid supermultiplet, respectively. Note that the mass splitting between the ground state and first excited state is approximately 1.2 GeV. This is in contrast to the 1S − 2S mass splitting of conventional charmonia, which is approximately 0.6 GeV. In other words, this r is a meaningful dynamical variable for charmonium-like hybrids.
• For the states in these four channels, we also obtain the BS wave functions Φ (r, t), with r here being the spatial distance between the charm quark fieldc and cB components. In 0 −+ and 1 −− channels, the lowest two states in each channel have the masses compatible with those of the 1S and 2S charmonia, and the r-behaviors of their BS wave functions O (II)) R (r, t) have the similar feature of the non-relativistic Schrödinger wave functions of cc systems in nS states. However, for higher states in the four channels, no radial nodes are observed in Φ (II) n (r). If these higher states are tentatively assigned to be hybrids, this observation might imply that the r here is less meaningful to describe the internal motion of hybrid than the distance between the cc component and chromomagnetic field strength B.
These results can be interpreted as follows. Since a charm quark is heavy and if the relativistic effect is not important, the Φ (I) n (r) function defined through a type-I operator may be considered the approximation of the radial wave function of a charmonium-like hybrid state to some extent. Thus, Φ (I) n (r) implies that the color octet cc pair resumes a center-of-mass motion recoiling against some additional degrees of freedom, which are necessarily gluonic in the quenched approximation. We consider the gluonic degree of freedom to be a 'constituent' gluon in the chromomagnetic mode, which functions as a color octet source and provides a potential in the non-relativistic picture. A more conceptually reasonable picture is considering the charmonium-like hybrids to be a color octet cc pair dressed by a color halo composed of gluons. In both scenarios, the binding mechanism is the strong interaction between color octets. Previous lattice studies demonstrated that for a pair of static color charge and anti-charge, their interacting potential is Cornell type [31], where the subscript D indicates the color SU (3)  eigenvalue C D of the second order Casimir operator in the D representation. This is called "Casimir scaling". The static potential of a heavy quark-antiquark pair is usually expressed as V QQ (r) = V 0 − 4α s /3r + σr, where the string tension σ is proportional to C F = 4/3. Therefore, for the color octet charge in the adjoint representation A of the color SU (3) with C A = 3, we obtain σ A = C A /C F σ = 9σ/4. If the excitation of a ground state hybrid is primarily along the direction from the centerof-mass of the cc to the gluonic degree of freedom, as is manifested by the wave function Φ n (r), then the quantum number J P C will not change and the excitation energy will be larger than that of the charmonium system since σ A is 2.25 times larger. Although we cannot yet derive the precise relation of the excitation energy to σ A , this may explain qualitatively that the mass splitting (approximately 1.2 GeV) of the ground state and first excited state hybrids is larger than the 1S − 2S mass splitting (about 0.6 GeV) of charmonium. We observe that the "color halo" picture can be also applied to the strangeonium-like hybrids [32].
Obviously, the "color halo" picture we propose is clearly conceptually different from the flux-tube picture of the hybrid involving a heavy quark-anti-quark pair QQ based on the Born-Oppenheimer approximation [1,2,5,6,33]. In the Born-Oppenheimer approximation, the gluonic degrees of freedom are considered to be light and fast and distributed along the displacement of the quark-antiquark pair, whose non-trivial representation along with the spin of the gluonic degrees of freedom indicates the specific quantum number of the potential and subsequently the J P C of the hybrid. On the other hand, the QQ is bound by an effective potential induced by the excitation of gluonic degrees of freedom. If the QQ pair can be considered static color sources, the excited gluonic degrees of freedom obey the cylinder symmetry along the QQ axis, which result in an excited static potential denoted by Λ η , where Λ = 0, 1, 2, . . . is the projected total angular momentum of gluons with respect to the QQ axis and is labeled as Σ, Π, ∆ for Λ = 0, 1, 2 etx., η represents the combined parity (P ) and the charge conjugate (C) of gluon excitations with η = g, u for P ⊗ C = ±, respectively, and is the P parity of the glue state. Therefore, the quantum number of a QQ state with this type of potential is whereL =L QQ +Ĵ g withL QQ is the orbital angular momentum of QQ with respect to the midpoint of the QQ axis, andĴ g is the total angular momentum of gluons. The conventional Cornell-type potential between the QQ pair is indicated by the ground Σ + g potential thus the P and C quantum number reproduce the conventional quantum number. The lowest 1 −− and (0, 1, 2) −+ hybrid supermultiplet is associated with the Π + u (L = 1) potential such that the radial Shrödinger equation can be expressed explicitly as where r is the distance between Q andQ, µ is the reduced mass of the QQ pair, and u(r) is related to the radial wave function φ(r) by u(r) = rφ(r). The effective potential V eff is where L 2 QQ = L(L + 1) − 2Λ 2 + Ĵ 2 g and Ĵ 2 g = 2. The eigenvalue E is independent of the total spin S of the QQ pair and thereby results in the 1 −− and (0, 1, 2) −+ supermultiplet. When the parameters of V QQ are set phenomenologically or from the lattice QCD results, the Shrödinger equation in Eq. (14) can be solved and the mass spectrum can be obtained. In Ref. [5], the mass splitting between the lowest Π + u multiplet and its first radial excitation is determined to be approximately 350 MeV for charmonium-like hybrids and 207 MeV for bottomium-like hybrids, which are significantly smaller than 1.2 GeV obtained in this paper. Moreover, the wave function with respect to the distance between the QQ pair should behave like a P -wave one, which is clearly different from the behaviors of the wave functions Φ (II) n (r). In other words, our results do not support the flux-tube description of heavy-quarkonium-like hybrids in the Born-Oppenheimer approximation. It should be emphasized that even though the interpretation of the wave functions can be debatable, the pattern of the spectrum should be solid and model-independent since it is derived directly from the lattice QCD calculation.
The "color halo" picture can have physical consequences. Based on the discussion above, the cc pair and gluonic component are bound through the potential of color octets. This binding can be easily broken by the excitation of gluons such that the cc component is neutralized in color and is emitted as a conventional charmonium, while the gluonic component is hadronized into light hadrons. That is, the decay modes of a charmonium state plus light hadrons can be important for the decay of a charmonium-like hybrid. This decay property seems compatible with the decay pattern of Y (4260) (now named ψ(4230) aka Y (4230) by the PDG), which is occasionally assigned to be a candidate for the 1 −− charmonium-like hybrid. Y (4260) was first observed in the invariant mass spectrum of J/ψπ + π − . Over the past several years, BESIII has observed structures at the center-of-mass of 4.22-4.23 GeV of the final states J/ψππ, χ c0 ω, h c ππ, ψ(3686)ππ in e + e − annihilations [22], which can be the previous Y (4260) if they come from the same state. The cross sections of these processes are comparable at the peak positions and can be understood by the hybrid assignment of Y (4260). Note that thecc pair in the 1 −− hybrid is a spin singlet, the decay modes involving h c and η c should be preferable to those involving J/ψ, χ c0 and ψ(3686) owing to the suppression of the spin flipping of charm quarks. However, the h c π + π − decay mode, in which h c and π + π − are in relative P -wave, is suppressed by the centrifugal potential barrier in contrast to the other channels in which the charmonium and light hadrons are in S-wave. The two effects compensate and result in comparable cross sections. On the other hand, the 1 −− states Y (4360) and Y (4660), if they do exist, are disfavored to be the radial excitations of 1 −− charmonium hybrids, since their masses are much lower than the predicted value in this work.
The (0, 1, 2) −+ charmonium-like hybrids can be also searched in the decay modes involving a charmonium state, among which the final state χ c0,1,2 η may be important, since no suppression occurs from the spinflipping of heavy quarks and the centrifugal potential barriers. The disadvantage of these modes is that η only has a small fraction of the flavor singlet component, but the QCD anomaly may enhance its production to some extent if it couples to two gluons. Other modes, such as J/ψω and J/ψφ (in P -wave), are worth considering, which are similar to the J/ψπ + π − and χ c0 ω decay modes of the 1 −− case. Since these states are heavy and cannot be observed directly in the e + e − annihilation, BelleII and LHCb may implement a mission to search for them in B meson decays.

VI. SUMMARY
The internal structures of the J P C = 1 −− , (0, 1, 2) −+ charmonium-like hybrids are investigated for the first time through their BS wave functions from lattice QCD in the quenched approximation, where the wave functions Φ (I) n (r) are defined by the state-to-vacuum matrix elements of spatially extended hybrid-like operators (type-I operators) with the color octetcc component separated from the chromomagnetic operator by a spatial distance r. After singling out the conventionalcc states in the 0 −+ and 1 −− channels, we confirmed the existence of a 1 −− and (0, 1, 2) −+ supermultiplet of nearly degenerate masses of approximately 4.3 − 4.6 GeV and similar BS wave functions Φ 1 (r) without nodes in the r direction. The first excited hybrid states also compose a supermultiplet with masses of approximately 5.7 GeV, and their BS wave functions Φ 2 (r) are almost the same and have one node. We checked these results by using type-II operators that the charm quark fieldc is separated from the cB component by a spatial distance. While the spectra from the two types of operators are consistent with each other, the wave functions Φ (II) n (r) do not have a nodal structure with respect to the distance betweenc and cB. These observations imply that r can be a significant dynamical variable for charmonium-like hybrids. The spectrum and information from the wave functions obtained in this study do not support the fluxtube description of heavy quarkonium-like hybrids in the Born-Oppenheimer approximation. Instead, we propose a "color-halo" scenario for the internal structure of the charmonium-like hybrids in which a relatively compact color octetcc pair is surrounded by gluonic degrees of freedom. Thus, the decay modes of a charmonium plus light hadrons are important for charmonium-like hybrids.