Identifying spin and parity of charmonia in flight with lattice QCD

The spectrum of charmonium resonances contains a number of unanticipated states along with several conventional quark-model excitations. The hadrons of different quantum numbers $J^P$ appear in a fairly narrow energy band, where $J^P$ refers to the spin-parity of a hadron at rest. This poses a challenge for Lattice QCD studies of (coupled-channel) meson-meson scattering aimed at the determination of scattering amplitudes and resonance pole positions. A wealth of information for this purpose can be obtained from the lattice spectra in frames with nonzero total momentum. These are particularly dense since hadrons with different $J^P$ contribute to any given lattice irreducible representation. This is because $J^P$ is not a good quantum number in flight, and also because the continuum symmetry is reduced on the lattice. In this paper we address the assignment of the underlying continuum $J^P$ quantum numbers to charmonia in flight using a $N_f = 2 + 1$ CLS ensemble. As a first step, we apply the single-hadron approach, where only interpolating fields of quark-antiquark type are used. The approach follows techniques previously applied to the light meson spectrum by the Hadron Spectrum Collaboration. The resulting spectra of charmonia with assigned $J^P$ will provide valuable information for the parameterization of (resonant) amplitudes in future determinations of resonance properties with lattice QCD.


I. INTRODUCTION
Over the past two decades many interesting resonancelike structures have been discovered in hadronic final states in the energy regime of heavy quarkonium (cc orbb). A large collection of these structures (generally referred to as XYZs) appears close to open flavor strong decay thresholds and does not fit into a simple non-relativistic quark-antiquark picture. Theoretically, resonance peaks are associated with pole singularities of scattering amplitudes in the complex energy plane. Hence their nature and properties have to be inferred from the respective scattering matrices and from their decays. Potential models and effective field theories, that approximate certain regimes of the strong interaction, provide insight into describing these resonance peaks in various different ways including as tetraquarks (typically compact diquark-antidiquark states), mesonic molecules, hadro-quarkonia or hybrid mesons with excited gluonic content. There are also interpretations in which these hadronic resonance structures arise from strongly coupled scattering channels. For a detailed review on the different theoretical investigations, see Refs. [1,2]. In the end, understanding these hadron resonances from first principles requires non-perturbative approaches such as lattice QCD.
So far, only one lattice simulation studied the charmonium resonances above open charm threshold, performing a finite volume analysis using Lüscher's method [9][10][11] and obtaining exploratory results for masses as well as decay widths in the vector and scalar channels [12]. In particular, all previous studies were performed in the rest frame, which provides only rather limited information on the scattering matrix. Lattice calculations of hadronic systems in moving frames (i.e. with non-zero total momentum) provide additional information on the relevant two-meson scattering matrices and have been quite successful in providing a comprehensive understanding of various hadron resonances (for a review, see Ref. [13]).
The first step towards a detailed finite volume analysis is a rigorous extraction of the excited hadron spectrum on the lattice in multiple inertial frames, followed by a reliable identification of the continuum quantum numbers. In this work, we investigate the charmonium spectrum on the lattice in the lowest three inertial frames with the square of the total spatial momenta, |p| 2 = 0, 1 and 2 (in units of (2π/L) 2 , where L is the spatial extent of the lattice). Compared to the rest frame, the determination of the spectrum on the lattice in moving frames poses additional challenges. This holds, in particular, for the charmonium spectrum, since there are a number of states with different spin (J) and parity (P ) in a narrow energy region and hence disentangling the J P information becomes non-trivial.
The main aim of this paper is to extract the charmo-arXiv:1811.04116v1 [hep-lat] 9 Nov 2018 nium spectrum on the lattice in the rest frame and in moving frames and to discuss our procedure for determining the underlying J P . Note that spin and parity refer to continuum quantum numbers of a hadron in its rest frame. Experimentally, J P for hadronic resonances are determined from partial-waves of their decay products transformed to their center of momentum frame. On the lattice, this information has to be inferred from the hadron energy spectrum and the overlaps with the lattice interpolating fields. This is rather challenging as lattice eigenstates transform according to irreducible representations (irreps) of the reduced symmetry group on the lattice. Hence infinitely many continuum states are constrained to appear in the spectrum of a finite number of lattice irreps. This means that continuum 1 eigenstates with different J contribute to a given irrep of the octahedral group (O h ) on a cubic lattice. For p = 0, the symmetry is further reduced to the little groups, implying that the continuum eigenstates not only with different J but also with different parity, P , contribute to a given irrep of the respective little group. For example, the only irrep of the Dic 2 little group (corresponding to p = (1, 1, 0)) that contains scalar charmonia with J P = 0 + , also contains states with J P = 1 − , 2 ± , 3 ± , ... . Thus the charmonium spectrum on the lattice can be quite dense -with several states of different J P contributing to a single lattice irrep. A rigorous spin-identification procedure for the charmonium spectrum on the lattice is therefore desirable.
To identify the spin and parity, we follow the approach proposed by the Hadron Spectrum Collaboration (HSC) and first applied to light iso-vector mesons at rest [14] and in flight [15]. We construct correlation functions from single meson interpolators on an N f = 2 + 1 ensemble generated by the "Coordinated Lattice Simulations" (CLS) consortium [16,17] with m π 280 MeV, L 2.1 fm and lattice spacing a 0.086 fm. Using these correlation functions, we illustrate the spin identification procedure that allows us to determine the underlying J P . For practical reasons we limit our discussion to J ≤ 3 (and consequently to helicities |λ| ≤ 3) in all irreducible representations for the three inertial frames considered. Some preliminary results have been presented in Ref. [18]. This spin identification procedure can be straightforwardly applied to excited charmonium studies on different lattices.
We restrict ourselves to the single-hadron approach employing only quark-antiquark interpolating fields and assuming that this procedure qualitatively captures the single-meson spectrum in the finite volume. In particular this approach does not aim at extracting the full spectrum of multi-hadron scattering states in finite volume. We briefly address indications for the presence of mesonmeson admixtures and postpone a rigorous extraction of scattering amplitudes from a basis of meson-meson and quark-antiquark interpolators to future calculations. We emphasize that realizing this technology for multiple inertial frames within the single-hadron approach is a valuable first step in a precise determination of the complete finite-volume spectrum, which will then be used in the extraction of resonance information [15].
We are currently performing an analysis of excited charmonia beyond the single-hadron approach on CLS lattices, including the ensemble used here; preliminary results were presented in Ref. [19]. This study also considers charmonium systems with non-zero momenta. The results from the present study can inform (model) parametrizations of scattering amplitudes with narrow resonances.
The layout of this paper is as follows. In Section II, we describe our spin identification procedure for charmonium at rest and in the moving frames. The construction of meson interpolators used in this study is detailed in Section III and our lattice setup is briefly outlined in Section IV. Section V illustrates our spin identification procedure with examples and our conclusions are presented in Section VI.

II. CHARMONIA AT REST AND IN FLIGHT ON THE LATTICE
Energy eigenvalues of charmonia, E n , are extracted on the lattice by computing two-point correlation functions constructed from annihilation (creation) operators O ( †) i,Λ C of typecc, inserted at initial and final times, t and t f , respectively: The interpolators (O † j,Λ C ) are chosen to transform according to the irreps Λ C of the lattice symmetry group. Subscripts i and j are indices referring to the operator label. The spectral decomposition contains all states that transform according to the lattice irrep Λ C . Charmonia are isosinglets and charge parity C is a good quantum number at rest as well as in flight, also on the lattice. The overlap of eigenstate n with an operator with label i is given by Z n i = O i,Λ C |n . A tower of low lying energy levels for a given channel can be determined by employing a large basis of interpolators to form correlation matrices and solving the generalized eigenvalue problem (GEVP) [20][21][22]. As states with different rest frame quantum numbers J P contribute to correlators corresponding to a given Λ C of the reduced symmetry on the lattice, assigning values of spin and parity to these levels is difficult. We overcome these difficulties by adopting an approach developed in Refs. [14,15], where interpolators well-suited for the GEVP, while also aiding the assignment of spin, are constructed.
In the following, we consider in detail the spin-parity assignments of levels with J ≤ 3, first for the more familiar case of p = 0 and then for the more challenging case of p = 0. In each case, we briefly discuss the symmetries and associated quantum numbers relevant in the continuum and on the lattice (which determine which J P can contribute to a given lattice irrep). Then we summarize our spin identification procedure based on energies and Z factors.

A. Mesons at rest
Spin J, its z-component M and parity P are good quantum numbers for a meson at rest in the infinite volume continuum. The corresponding eigenstates are labeled as |p = 0; J P C . On the lattice P remains a good quantum number. Lattice eigenstates transform according to one of the five irreps of the octahedral group O h , given in Table I. Group-theoretical subduction relates the infinite number of continuum spins to the finite number of lattice irreps as shown for J ≤ 3 in Table I. Since several J contribute to each lattice irrep Λ, the spin assignment of lattice energy levels within each Λ C is not immediate.
Following Ref. [14], we first construct interpolators in the continuum with good quantum numbers J, M and P at p = 0 (O J P C ,M i ). Using these continuum interpolators, one then builds the lattice interpolators O For details, see Section III. In this way these lattice interpolators are expected to have good overlaps with states of the respective continuum quantum numbers, aiding their clean determination. Below we outline the 'guidelines' we follow in our spin assignment procedure for mesons at rest: • Trivial assignments: According to Table I, con-tinuum states with J = 0 are expected to appear in the A 1 irrep. Assuming negligible contributions from any J ≥ 4 spin state in the energy range considered, all extracted levels in the A 1 irrep can be trivially assigned J = 0. Analogously, the levels in A 2 and E can only result from J = 3 and J = 2 continuum states, respectively.
• (Non-)degeneracy of energy levels: The appearance of a set of near degenerate energy levels (up to discretization and finite volume effects) or the absence of any near degenerate partner level across different lattice irreps immediately suggests possible spin assignments. For example, the appearance of a level in the T 1 irrep of O h with no near degenerate partner levels in any other irreps suggests a J = 1 spin assignment (c.f. Table I).
Similarly, three near degenerate energy levels, one each in the spectrum of the T 1 , T 2 and A 2 irreps of O h , suggest a J = 3 assignment.
• Enhanced Z factors: Those lattice eigenstates resembling continuum states with quantum numbers J P C are expected to have the largest overlap factors with the lattice operators, O i,Λ C |n suggest a spinparity assignment J P . To make such an assessment, we utilize the quantityZ n i , which is normalized to unity with respect to the largest Z n i for a given operator with label i across all lattice states If the effects of rotational symmetry breaking are sufficiently small, one expects for J = J.
• Degeneracy in Z factors for spin J > 1: One expects similar Z factors for a continuum state |p = 0; J P C with interpolators subduced from the same continuum operator For example, a J = 3 continuum state manifests as three almost degenerate energy levels, one each in the spectrum obtained using T 1 , T 2 and A 2 operators. Each of these levels is expected to have degenerate Z n i factors for lattice interpolators subduced from the same continuum interpolator with spin 3. i.e. O [3] i,T1 |n T1 O [3] i,T2 |n T2 O [3] i,A2 |n A2 , if |n T1 , |n T2 and |n A2 represent the same state in the continuum. In the above example we have suppressed the P and C indices for clarity.
3 ± TABLE I. Lattice irreps Λ for the symmetry groups corresponding to the momentum p = (0, 0, 0), (0, 0, 1) and (1, 1, 0), and continuum quantum numbers that can contribute to these lattice irreps. For the rest frame, the table presents the distribution of J across different lattice irreps. Only values of J ≤ 3 are listed, as we do not consider higher values of spin in this work. For p = 0, the table lists, which helicities (λ) and whichη = P (−1) J can appear in each lattice irrep;η is a good quantum number only for λ = 0. In the third column, the J P that can appear in the moving frame irreps are also shown.

B. Mesons in flight
For mesons with non-zero momentum in the infinite volume continuum, the O(3) symmetry group is broken to its subgroup, U (1). The eigenstates |p, J P C , λη are labeled by the magnitude of helicity λ (the projection of the spin component along the direction of momentum p/p) andη parity defined asη = P (−1) J .η is a good quantum number only for λ = 0 and is related to the reflection in the plane containing p. Unlike in the rest frame, J P are no longer good quantum numbers and the hadron spectrum in moving frames will be a mixture of different J P for any given λ [15]. Hence, discerning the J P information of the spectrum in moving frames requires more care than for hadrons at rest.
On the lattice, λ andη are no longer good quantum numbers. Lattice eigenstates transform according to a given lattice irrep Λ C of the reduced discrete lattice symmetry group associated with a particular spatial momentum. We consider p = (0, 0, 1) (in units of 2π N L ) with symmetry group Dic 4 and p = (1, 1, 0) with symmetry group Dic 2 . Table I presents all irreps for both cases and lists the values of λ andη (for λ = 0) that can contribute to each of these lattice irreps. Note that we only consider J ≤ 3 (|λ| ≤ 3) in this work, and that the spectrum of every irrep of Dic 2 and the E irrep of Dic 4 receive contributions from two different helicities.
Knowing that a meson with spin J can have helicities 0 ≤ |λ| ≤ J, we list the J P of states that can contribute to each irrep in the rightmost column of Table I. For example, the A 1 irrep of Dic 4 contains only |λ| = 0 with an allowed value ofη = +1. Thus mesons expected to appear in A 1 with |λ| = 0 can only have J P = 0 + , 1 − , 2 + , 3 − . For A 1 in Dic 2 , we arrive at the same pattern of J P for |λ| = 0, however, states with |λ| = 2 can also contribute. In this caseη is not a good quantum number and hence there is no restriction on the parity. States with J ≥ |λ| = 2 are possible, leading to the allowed combinations J P = 2 ± , 3 ± .
Clearly, identifying the underlying J P of energy levels extracted in each lattice irrep is non-trivial when studying systems in flight, as a number of states with different J P can appear in the spectrum for each lattice irrep. To aid this study, we follow the procedure in Ref. [15], which is detailed in Section III. Using the continuum interpolators O J P C ,M i (p = 0), we construct operators at p = 0 with good helicity λ (O J P C ,λ i (p)). The lattice interpolators respecting the reduced symmetry of the lattice (O [J P C ,|λ|] i,Λ C ) are built from these moving frame continuum interpolators and are therefore expected to have strongest overlap with continuum states with quantum numbers λ andη. Hence our spin identification proceeds in two steps : 1) identification of the helicity of lattice energy levels (analogous with the J identification for the spectrum at rest), 2) the spin and parity assignments. Below we summarize our 'guidelines' for spin-parity assignments in a moving frame :
• (Non-)degeneracy of energy levels: Appearance of near degenerate energy levels (up to cutoff and finite volume effects) or absence of a near degenerate partner level across the lattice spectrum according to Table I suggests possible helicity assignments.
• Z factor degeneracy for non-zero helicities: Similar Z factors are expected for a continuum state |p; J P C , λ with interpolators subduced from the same continuum operator For example, a |(1, 1, 0), 3 P C , 3 continuum state manifests as two near degenerate energy levels, one each in B 1 and B 2 of Dic 2 . Both these levels should have degenerate Z n i factors for lattice interpolators subduced from the same continuum interpolator with λ = 3. Hence if |n B C 1 and |n B C 2 represent the same λ = 3 state in the continuum.
2 Identifying J P from known |λ|η • Helicity components: The spin of the state with a given helicity is constrained to be J ≥ |λ|. Different |λ| components of a continuum J P state are constrained to have the same energy by Lorentz symmetry [15]. These constraints suggest J P assignments by matching the observed pattern of energy levels with the allowed distribution of helicity components of a continuum J P as summarized in Table I.
• Information inZ factors: At rest a continuum interpolator O J P ,|λ| i (p) can only overlap with states of spin-parity J P , whereas at non-zero momenta this operator can also overlap with states that have other spin and parity [15]. If the employed momentum and the effects of breaking of the continuum symmetry on the lattice are small, one expects to have dominantZ factors for O [J P ,|λ|] i,Λ (p) with states that have quantum numbers J P in the continuum at rest. We utilize this expectation not only to confirm the initial spin-parity assignments, but also to identify the quantum numbers for remaining ambiguous levels.

III. INTERPOLATORS
Only quark-antiquark (cc) interpolators are considered in this work. To construct interpolators with good overlaps to physical states, we follow the methods of Refs. [14] and [15] for the rest frame and moving frames, respectively. The approach involves first building interpolators in the continuum with angular momentum, J, its z-component M and parity P at p = 0. Using these rest frame interpolators, operators at p = 0 and helicity λ are constructed. These continuum interpolators are then subduced onto their lattice counterparts using the respective subduction coefficients. If the effects of a finite lattice spacing are small, these subduced lattice operators will have a good overlap with the states in the continuum. For completeness, we outline the main steps below and refer the reader to Refs. [14,15] for further details.
The starting point is a basis of continuum quark bilinear interpolators with spin structures built from gamma matrices and single and double covariant derivatives, Similarly using one-derivative bilinears in Eq. (8), one can construct J ≤ 2 operators and J ≤ 1 interpolators can be built using the local bilinears. At p = 0 the projection onto lattice interpolators for each lattice representation Λ proceeds via subduction co- where µ indicates the row of the lattice irrep and operators on the r.h.s are understood to be discrete versions of the continuum ones. See Appendix A of Ref. [14] for subduction coefficients for J ≤ 4. At finite p one first forms continuum operators of definite helicity λ: where D (J) * M,λ (R) is the Wigner-D matrix and R refers to the rotation from (0, 0, |p|) to p. For the trivial case of p = (0, 0, |p|), D As discussed in Ref. [15], in order to ensure consistency between different momentum directions, R is split into two parts R = R lat R ref : a reference rotation R ref that takes (0, 0, |p|) to p ref and a lattice rotation R lat from p ref to p. We take p ref = (0, 0, 1) and (1, 1, 0) for Dic 4 and Dic 2 , respectively. Finally, these continuum operators are subduced to obtain the lattice interpolators The subduction coefficients, Sη ,λ Λ,µ for Dic 2 and Dic 4 can be found in Table II of Ref. [15]. We construct the interpolators in all momentum polarizations. The correlation functions of different momentum polarizations for a given |p| are averaged before the GEVP analysis (discussed in section IV).
To illustrate the above procedure, consider the example of an interpolator in the (one-dimensional) B 2 irrep of Dic 4 (i.e. p = (0, 0, 1)), which receives contributions from λ = ±2 (see Table I

IV. LATTICE SETUP
The charmonium spectrum has been determined on an N f = 2 + 1 ensemble produced by the CLS consortium [16,17], labeled as U101. The configurations have been generated with a non-perturbatively O(a) improved Wilson fermion action and a Symanzik gauge action. The pion and the kaon masses are m π 280 MeV and m K 467 MeV respectively. The U101 ensemble lies along the Tr(m) = 2m u/d + m s = const line, meaning that the strange quark mass becomes heavier and the light quark mass lighter as the physical point is approached. The volume of the lattice is 24 3 × 128 and the lattice spacing is a = 0.08636(98)(40) fm [23] giving a physical spatial volume of (2.07 fm) 3 . The gauge and fermion fields fulfill open boundary conditions in the time direction and translational invariance is only expected in the bulk [24]. Hence we measure two-point correlation functions in the middle of our lattice at least 28 time slices away from the boundaries. No effects related to the finite temporal extent are seen in this time interval for the pion and charmed meson correlation functions. We utilize a total of 1638 source time slices from 255 configurations to compute the correlation functions.
The charm quark is quenched in our simulations and vacuum charm quark loops are absent in the theory. The charm quark mass is therefore set independently when measurements are performed. There are various possible observables that can be used for tuning the charm quark mass. Our long-term objective is the determination of the properties of exotic charmonium resonances, therefore the position of theD-D decay threshold needs to be taken into account. In particular, a pion mass heavier than its experimental value might result in a situation where a physical resonance will appear as a bound state. Another important issue is how the properties of the exotic charmonium states depend on the charm quark mass. We have therefore determined the charmonium spectrum for two different values of the charm quark mass, corresponding to κ c = 0.12315 and κ c = 0.12522, resulting in D-meson masses approximately 80 MeV above and below the physical value, respectively. In the present work, we focus on the results from κ c = 0.12522, and we will utilize the data from κ c = 0.12315 for our future meson scattering analysis.
We employ the distillation method [25] to compute matrices of correlation functions between a basis of meson interpolating fields. This method is powerful enough to implement local and non-local hadron interpolators as well as multi-hadron interpolators, and to perform definite non-zero momentum projections at the source and at the sink, which are crucial in a study of scattering amplitudes and resonances. This is achieved with a reasonable amount of computational resources and disk space requirements. Distillation is equivalent to standard Gaussian quark field smearing algorithms, written in terms of the eigenvectors of the gauge covariant lattice Laplacian in three dimensions. We construct our fermion sources from 90 eigenvectors and we employ full distillation, meaning that we compute the quark propagation between the full set of eigenvectors at the source and at the sink.
The computation of a given entry of the correlation matrix proceeds in four steps. First, the Arnoldi algorithm provides the smallest eigenvectors and eigenvalues of the Laplacian operator on all time slices t fulfilling the equation where spatial and color indices are suppressed and n labels the eigenvector index. Second, we solve the Dirac equation to obtain the quark propagators η k (t , α ) from each source ψ k (t, α). The source is zero everywhere except for time slice t and spin component α, where it is equal to the eigenvector Ψ k (t). The perambulators τ (t, t , α, α , k, k ) are constructed from the η k (t , α ) using In the third step, we compute the matrices φ The operator O αα (k, k ) has a spin structure and possibly also covariant derivatives acting on the Laplacian eigenvectors, therefore in general φ is not diagonal in the distillation space. The matrices φ carry the information about the quantum numbers of the operators that define the correlation function. Both φ and τ can be considered as squared matrices of dimension (90 × 4) × (90 × 4) defined on each time slice or connecting two time slices. From these matrices, we can finally compute the sum of all the Wick contractions, neglecting only the charm annihilation diagrams that are OZI suppressed. For instance, a simpleqq-meson correlation function reads explicitly The correlation matrices (Eq. (1)) are computed for all irreps in Table I usingcc operators constructed as in Section III. In Table II, we list the number of interpolators employed in each of the irreps in the three inertial frames we study. For multi-dimensional irreps, we restrict the computation of correlation matrices to one single row of the irrep. For brevity, we drop the subscript µ from the lattice operators used in Eqs. (10) and (12) throughout this section. Energies E n and overlaps Z n i = O i,Λ C |n of the lattice levels are extracted from solutions of the GEVP with t 0 = 2. E n are extracted by one or two exponential correlated fits to λ n (t), whereas Z n i s are extracted from constant fits to a plateau in the function, The quality of the fits to λ n (t) is illustrated for the example of the E − irrep of Dic 4 in Figure 16 of Appendix A. We illustrate the plateaus in Z n i for selected interpolators and the respective fits for the example of n = 1, 2 and 6 levels in the spectrum for the E − irrep of Dic 4 in Figure 17 of Appendix A.

V. RESULTS
In this section, we present the charmonium spectra for all lattice irreps in the three inertial frames considered. Using this energy spectrum and the corresponding operator state overlaps, we illustrate the spin identification procedure outlined in Section II and present the spin-identified charmonium spectrum in the three frames studied.

A. Charmonia at rest
In Fig. 1, we present the spin-identified charmonium spectrum in the rest frame for all lattice irreps, Λ P C . The color coding for the spectrum is similar to that used in Ref. [15], allowing for an immediate comparison of the qualitative features : The spin assignment of the charmonium spectrum in the rest frame is relatively straightforward and follows immediately from Table I and from the guidelines in Section II A.
Trivial assignments : Among J ≤ 3, the spectrum in the A 1 irreps can contain only J = 0 states, hence J = 0 assignment for these levels is immediate. Any effects of higher spins in these levels can only be due to J ≥ 4 states, which we assume to be negligible for two reasons: Firstly, no such low lying high spin charmonium states have been discovered. Secondly, we have not employed interpolators with such continuum quantum numbers. Analogously, levels in the spectrum of the E and A 2 irreps can only be assigned J = 2 and J = 3, respectively.
Energy degeneracies : Continuum states with J = 1 can appear only in T 1 . Hence initial J = 1 assignments can be made for the lowest two levels in T 1 irreps (one in T −+ 1 , which is an exotic quantum number), as they appear alone with no equivalent partner levels across other lattice irreps with the same P and C. The next higher spin that can occur in T 1 is J = 3, with near degenerate energy levels appearing in the T 2 and A 2 spectra. Thus among the third and the fourth levels in T −− 1 (aE n ∼ 1.57 in Fig. 1), the fourth level can be assigned J = 3, as is supported by the near degenerate energy levels in T −− 2 and A −− 2 . J = 1 is assigned to the third level with no remaining near degenerate partner levels in T −− 2 and A −− 2 . In general, the degeneracy of energy levels representing the same continuum state provide a handle on initial spin assignments of states with J ≥ 2.
Z factor enhancements : Lattice levels resembling J ≥ 2 often appear in dense energy bands along with states with different spins. In such cases, disentangling the quantum number information solely from the energies can be delicate. Initial spin assignments need to be further confirmed by investigating theZ factors (Eq. (2)). In Fig. 2 from the T −− 1 irrep to demonstrate the spin assignment for the third and fourth levels. The x-axis labels refer to the continuum interpolators O J P C , from which the respective lattice interpolators are built. Other indices are suppressed for brevity. Clearly, the third level (n = 3) has the largest and dominantZ factors for lattice interpolators subduced from J = 1 continuum interpolators, so this state is assigned spin J = 1. The fourth state (n = 4) has dominantZ factors for lattice interpolators subduced from J = 3 and is assigned spin J = 3.

, we show theZ factors for a subset of operators
Z factor degeneracies : By studying the expected degeneracies in Z factors (Eq. (4)), one can further confirm the J ≥ 2 assignments. In Fig. 3, we show such and E ++ . This confirms a J = 2 assignment for the levels at aE n ∼ 1.431 in these two irreps (c.f. Fig. 1).
By considering the pattern of energies and operator state overlaps of the lattice spectrum in the rest frame, we arrive at the spin-identified charmonium spectrum as depicted in Fig. 1. Within the single hadron approach, we assume the lattice energy levels represent the single hadron spectrum neglecting any effects of allowed strong decay thresholds. The interpretation of energy levels close to scattering thresholds, such as the first excitation in the T ++ 1 spectrum, are subtle. The presence of the nearby non-interacting levels DD * might have a strong influence on the determination of the energy of these excitations. We postpone these discussions to a follow up publication.

B. Charmonia in flight
|λ|−assignments of the lattice levels can be made following the guidelines discussed in Section II B. This is straightforward if the different helicity sectors are decoupled, i.e. if the breaking of U (1) symmetry is small. As the lattice interpolators for mesons in flight are subduced from continuum interpolators with good helicity, this is indeed observed in our correlation matrices. To demonstrate this, we present the normalized correlation matrix (C ij / C ii C jj ) on time slice 5 for the A + 1 irrep in Dic 2 in Fig. 4. A linear color gradient is used from deep blue representing the largest value (unity) to white representing the smallest value (zero). On the left, we present the correlation matrix with the operators ordered based on the J of the continuum rest frame interpolators used to build the continuum helicity interpolators. Strong crosscorrelations between different J sectors are evident from the figure. In the right figure, we order the interpolators such that λ = 0 (operators from 1 to 17) are collected before |λ| = 2 (operators from 18 to 28). Clearly the correlation matrix is almost completely block diagonal in helicity, indicating a clean signature for the helicity of the lattice levels. Hence in what follows, we first discuss the determination of the helicity-identified lattice spectrum. Then, we discuss our J P assignments based on the pattern of energies and overlaps.
the Dic 4 symmetry group, are relatively straightforward. According to Table I, λ = 0 is the only helicity appearing in A 1 or A 2 among |λ| < 4. Thus all levels in the A 1 and A 2 irreps in Dic 4 have helicity 0. Similarly |λ| = 2 is the only helicity appearing in the B 1 or B 2 irreps among |λ| < 4. Hence all the levels observed in B 1 and B 2 can trivially be assigned |λ| = 2. Linear combinations of λ = ±2 continuum states are expected to appear as pairs of energy levels, degenerate up to the effects of the reduced symmetry on the lattice (and finite volume effects), one each in the spectrum of B 1 and B 2 . Thus the spectra of these two irreps are expected to be similar, with all levels representing |λ| = 2 continuum states. The only non-trivial situation in Dic 4 arises in the E irreps, where |λ| = 1 as well as 3 can contribute. Disentangling the helicity 1 and 3 levels requires information from theZ i n factors as will be discussed shortly.
Energy degeneracies : In the charmonium spectrum with momentum p = (1, 1, 0), i.e. the Dic 2 symmetry group, each of the four lattice irreps can have contributions from two different helicities. Hence there are no trivial helicity assignments possible. According to the subduction patterns in Table I for the Dic 2 frame, a λ = 0 state can appear only in either A 1 or A 2 , depending on the underlying J P of the continuum state. Linear combinations of the λ = ±2 continuum states appear in the spectrum of both A 1 and A 2 , degenerate up to the effects of reduced symmetry on the lattice. Hence appearance of energy levels in the A 1 (A 2 ) spectrum with no near degenerate partner levels in the A 2 (A 1 ) spectrum indicates a λ = 0 assignment. However, appearance of a pair of near degenerate energy levels one each in both A 1 and A 2 could either indicate a possible |λ| = 2 assignment or an accidental degeneracy of two λ = 0 levels. A |λ| = 2 assignment for these levels implies a possible minimum spin assignment, i.e. J ≥ 2. In this case, one should also find another level corresponding to the λ = 0 component nearly degenerate with the |λ| = 2 candidates, in either A 1 or A 2 depending on the J P of the state in the continuum. Thus the levels around aE n ∼ 1.47 in the spectrum of A + 1 and A + 2 could be composed of a pair of |λ| = 2 levels one each in A + 1 and A + 2 and a λ = 0 level in A + 1 . Similarly one can make |λ| assignments for levels around aE n ∼ 1.61 in the spectrum of A + 1 and A + 2 . The helicity assignment for the levels around aE n ∼ 1.6 in the spectrum of A − 1 and A − 2 is more complicated and requires some guidance from the spectrum of B − 1 and B − 2 or inputs from a study of overlap factors. Similar to |λ| = 2, linear combinations of states with λ = ±1 and λ = ±3 appear as near degenerate levels in the spectrum of both B 1 and B 2 in Dic 2 . Consequently they are expected to mimic each other, similar to the case of B 1 and B 2 irreps in Dic 4 . However, they will be a mixture of |λ| = 1 and |λ| = 3 continuum states. A pair of near degenerate energy levels, one each in the spectrum of both B 1 and B 2 , with no other near degenerate levels suggests a |λ| = 1 assignment. Thus the lowest three levels in the spectrum of B 1 and B 2 for both C-parities can be assigned |λ| = 1, indicating also a possible minimum spin assignment J ≥ 1. A non-zero spin assignment implies the presence of another near degenerate λ = 0 level in either A 1 or A 2 , depending on the J P of the continuum state. A possible |λ| = 3 assignment to any pair of levels in the spectrum of B 1 and B 2 is indicated by appearance of additional near degenerate levels corresponding to the λ = ±1 components of the same continuum state at rest, similar to what was argued in the case of |λ| = 2. However, such an assignment is complicated as a λ = ±1 state appears in both B 1 and B 2 , unlike the case of a λ = 0 state, which appears only in either A 1 or A 2 . Beyond this point, a study of Z factors becomes crucial.
A flavor of J P assignments is already evident from these arguments on the expected degeneracies between different helicity components of the same continuum state across different lattice irreps. In the discussion above, they are intended to affirm the initial helicity assignments based on energy degeneracies.  Z factor enhancements : A study of Z factor enhancements on this partly helicity assigned spectrum not only verifies initial |λ| assignments based on energy degeneracies, but also disentangles the helicity composition of all energy levels in the spectrum in a transparent way. To demonstrate this, we showZ n i factors for the sixth lattice level for a subset of interpolators in the A + 2 and B − 1 irreps in Dic 2 in Fig. 5. The levels in the A + 2 and B − 1 irreps can clearly be identified to have dominantZ factors with lattice interpolators, O λ=0 and O |λ|=1 , respectively 4 . In this way one can make reliable helicity assignments for all the levels in the lattice spectrum.
Z factor degeneracies : Similar to the study in the rest frame, one can further investigate degeneracies in Z n i factors of pairs of near degenerate energy levels assigned with the same helicity across different lattice irreps (Eq. (6)). In Fig. 6 6. A comparison of the Z n i factors to confirm the quantum number assignments. Shown are the Z factors for two continuum interpolators that determine the levels close to aEn = 1.61 in all four irreps within Dic2.
to confirm the quantum number assignment for five levels close to aE n = 1.61 in the C = + spectrum of all four irreps within Dic 2 (c.f. Fig. 12). Plotted are Z factors of these states for lattice interpolators across all lattice irreps originating from two rest frame continuum interpolators, O 2 −+ i=1,2 . As summarized in Table I, λ = ±2 components of O 2 −+ i=1,2 are subduced to A + 1 and A + 2 , while its λ = ±1 components are distributed on to B + 1 and B + 2 . The λ = 0 component of this interpolator falls into the A + 2 irrep, as it has negative parity. Since lattice interpolators are built from moving frame continuum operators with good helicity, Z n i factors for the levels representing the same continuum state |λ| across different irreps are expected to be degenerate (Eq. (6)). It is evident from the figure that Z n i factors of the level in A + 1 and one level in A + 2 are degenerate, indicating a |λ| = 2 assignment. Similarly, the degenerate Z n i factors for the levels in B + 1 and B + 2 indicate a |λ| = 1 assignment. The only remaining level in A + 2 with no partners to associate can trivially be assigned λ = 0. Following these procedures, we arrive at the charmonium spectrum, as shown in Figs. 10 and 12, with reliably identified helicities.
Identifying J P from known |λ|η Once we have the helicity-identified charmonium spectrum, we can proceed to spin-parity assignments, which we perform considering the helicity components and Z factors as described below.
Helicity components : A meson with continuum J P can have helicities |λ| ≤ J. These helicity components are expected to be distributed across the lattice irreps as per Table I. The reduced symmetry of moving frames in the continuum (U (1)) does not constrain energies of different |λ| to be equal. However, constraints from Lorentz symmetry enforce different helicity components of a continuum J P to have the same energy [15]. Thus the easiest way to make J P assignments is by investigating the distribution of helicity-identified lattice energy levels across the little group irreps (Table I) and matching them with expected degeneracies.
First, we discuss the simple case of identifying the isolated 2 −+ state represented by five levels close to aE n = 1.61 in the C = + spectrum of all four irreps within Dic 2 (c.f. Fig. 12). Clear identification of all expected helicity components up to 2 with no remaining energy levels to consider and remarkable degeneracy in their energies, as evident from Fig. 6, suggest all those lattice levels to represent the same J = 2 continuum state. Furthermore, a J = 2 assignment plus the λ = 0 level appearing in A + 2 indicate a P = − assignment for the continuum state in the rest frame. Alternate spin assignments are excluded due to the presence of helicity 2 levels. The presence of any additional continuum state is excluded as all five levels under consideration are exhausted in describing helicity components of the same J = 2 continuum state. Thus one arrives at a rest frame quantum number J P C = 2 −+ for the state in the C = + Dic 2 spectrum around aE n ∼ 1.61. Now we elaborate on our procedure for a more complicated example of multiple continuum states appearing as a band of levels in a narrow energy region. To this end, we consider the energy levels in the region aE n ∈ (1.57, 1.63) in Dic 2 with C = − as highlighted in  The only remaining quantum number of these two states to be inferred is parity, which is determined bỹ η = P (−1) J as discussed earlier. However, this is not immediate as it is not evident which of the λ = 0 levels in A − 1 and A − 2 corresponds to the J = 2 and J = 3 continuum state. Nevertheless, an immediate inference from these two λ = 0 levels is that both the J = 2 and J = 3 states should have the same parity. This is independent of the final spin assignments of these λ = 0 levels. This means if the λ = 0 level in A − 1 is related to the J = 2 continuum state, then both J = 2 and J = 3 states have positive parity. If instead the λ = 0 level in A − 1 is related to the J = 3 continuum state, then both of them have negative parity. In other words, the continuum J P C composition of energy levels in this band can be either (2 +− , 3 +− ) or (2 −− , 3 −− ). However, due to the exotic nature of J P C = 2 +− , a P C = +− assignment is less plausible at this energy range. Hence, we make a (2 −− , 3 −− ) assignment for these levels. For the levels identified with λ = 0 and |λ| = 3 the spin assignment is complete, whereas for the levels identified with |λ| = 1 and |λ| = 2 the spin assignment (between J = 2 or J = 3) requires additional information contained in theZ factors, along the lines described in the next subsection.
for J P = J P . To illustrate this, we present theZ factors for the lowest six levels in A + 2 of Dic 2 and a subset of interpolators in Fig. 8 that can overlap with states of quantum numbers J P in the continuum at non-zero momentum. This information can be deduced from Tables X-XV in Ref. [15] and the corresponding overlaps are proportional to p or to higher powers. Black is used for Z factors of operators that cannot overlap with quantum numbers J P in the continuum either at rest or with non-zero momentum. Hence theZ factors in red are expected to be dominant in comparison with those in blue and black. This is evident from the figure 5 . The largẽ Z factors in red unambiguously confirm the J P assignments for all the lattice levels we have extracted. Note that the information in Table X-XV of Ref. [15], which helps to differentiate the operators into blue and black, was not essential for identifying J P in Fig. 11 and Fig.  13. In this way, we confirm all previously assigned J P and identify the J P for all the remaining ambiguous levels. Thus taking information from the expected patterns of the energy spectrum in moving frames and from thẽ Z factors into account, we arrive at the J P C assigned lattice spectrum in the moving frames, as shown in Figs. 11 and 13. The color coding is given in Eq. (21).
We remark that the Z factors for different helicities of a continuum state with non-zero momentum are not constrained to be degenerate by the U (1) symmetry group. We indeed find statistically different overlaps for different |λ| (c.f. Fig. 6). Hence a comparative study such as in Fig. 3 is not useful in making the continuum spin assignments in a moving frame.
Consistency check using the dispersion relation: Finally, as a consistency check, we verify that the energies of the spin identified spectra from different inertial frames roughly follow a dispersion relation E(p). The charmed hadrons at our lattice spacing a are not expected to satisfy the continuum dispersion relation exactly. So we test the extracted energies based on a lattice dispersion relation cosh(bE(p)) = cosh(bM ) + i=x,y,z where M = E(0) is the energy extracted at rest and the parameter b can be viewed as an effective lattice spacing relevant for this momentum dependence. This dispersion relation interpolates between the continuum relation E(p) = M 2 + p 2 for b → 0 and the free-boson dispersion relation for b = a. Figure 9 shows an example, where the parameter b 0.46 a is fit from J P C = 1 +− and the same parameter is employed for other 1P charmonia, indicating that (24) results in a reasonable momentum dependence for those.

C.
Charmonium summary for various J P The J P C -identified charmonium spectra with various momenta within the single hadron approach are summarized in Figs. 14 and 15. Only the states with reliably identified spin and parity are shown and the energy values are based on Figs. 1, 11 and 13 6 .
5 TheZ is not significantly enhanced for the n = 4 state in A + 2 of Dic 2 . But its spin-parity J P = 0 − is clear, given the absence of near degenerate levels in other irreps (c.f. Fig. 12). 6 Eigenstates with certain J P C are present in several lattice irreps and we choose the average of energies from all the irreps for FIG. 9. The lattice energies at various p 2 for the 1P charmonium levels (red) together with the dispersion relation (24) at integer p 2 joined by lines (blue). The parameter b in the dispersion relation is fitted from the 1 +− charmonium and the same parameter is employed for the other states. Figure 14 shows the charmonium mass spectrum determined at rest in terms of mass splittings with respect to the mass of the η c meson. These mass splittings are compared to those from experiments [28]. Except for 1 −+ , the extracted levels are likely candidates for conventional quark model charmoniacc. The 1 −+ state is a candidate for the hybrid charmoniumcGc with exotic J P C , i.e. quantum numbers that cannot result from a quarkeach jackknife sample. Note that the levels in different irreps receive different shifts from the finite volume and the discretization used. Averaging those values provides a spectrum suitable for a qualitative comparison to experiment rather than a quantitative prediction.
antiquark configuration. The features of the lattice and experimental spectrum agree qualitatively, while a quantitative agreement is not expected due to several shortcomings in our calculation. First of all, we utilize a lattice QCD ensemble with heavier than physical pion mass (m π 280 MeV) and our charm quark mass corresponds to a D meson mass about 80 MeV lower than its physical value, as discussed in Section IV. This work ignores the resonant nature of charmonia above the thresholds (indicated for DD and D s D s by dashed lines) and the resulting mass estimates are naively expected to be correct to the order of the respective decay widths. We also ignore the effects of thresholds on the near-threshold states. This renders the first excited state 1 ++ state too high with respect to experimental X(3872), which is common to all lattice results within the single hadron approach. The features of the extracted lattice spectra in Fig. 14 also roughly agree with other lattice results based on the single hadron approach (all in the rest frame) aimed at excitedcc [4][5][6][7][8].
Now we turn to the spectra with non-zero momenta, as the identification of their spin and parity was the main purpose of our work. The J P assigned energy spectra of charmonium in different inertial frames are compared in Fig. 15. The energies increase with the momentum roughly as expected by the dispersion relation (Eq. (24)); an example of a more detailed comparison was shown in Fig. 9. One does not expect the dispersion relation E(p) to be respected exactly for states near and above thresholds, since energies at different momenta receive different shifts, depending on the positions of the discrete noninteracting two-meson states. This plot presents only the spectral levels with reliably identified J P and reflects the fact that this identification gets more challenging as momentum increases due to the reduction of the discrete symmetries on the lattice (Table I) as well as due to the degradation of signal-to-noise at larger energies. Therefore, some of the higher-lying charmonia could be reliably assigned J P only for the lower momenta.  10. |λ|-identified charmonium spectrum in the moving frame with p = (0, 0, 1). All irreps Λ C of the corresponding little group Dic4 are presented ( Table I) FIG. 11. J P -identified charmonium spectrum in the moving frame with p = (0, 0, 1). Irreps Λ C of group Dic4 are presented. The colors indicate J P of states according to the color-coding (21).

VI. CONCLUSIONS
This paper studies the spectra and spin-parities of charmonia at rest and in flight. The analysis is performed using the distillation method on a single N f = 2 + 1 CLS ensemble with m π 280 MeV and a somewhat smaller than physical charm quark mass.
In the continuum, the spin J and parity P of a hadron are good quantum numbers in its rest frame, while the helicity λ is a good quantum number for a particle in flight. In experiment, the J P of a hadron in flight is determined by making a Lorentz transformation of its decay products to its rest frame. In lattice QCD, there is no direct analogue of this procedure due to the reduced Lorentz symmetry.
We show a reliable procedure to determine the underlying rest-frame J P of charmonia in flight within lattice QCD. To identify J P of a certain eigenstate, it is not enough to consider this state in isolation. One should consider most of the eigenstates up to the desired energy in all lattice irreducible representations. The challenge is that eigenstates with many different J P contribute to a given irreducible representation. The strategy is to determine the energies of all these eigenstates and their overlaps to lattice versions of carefully-constructed operators O [J P ,|λ|] (p) cc [14,15], which in the continuum couple only to charmonia with spin-parity J P at p = 0 and to states with helicity λ at p = 0. The J P are determined by considering degeneracies of energies or overlaps across different irreducible representations, and by considering the sizes of appropriately normalized overlaps. In this way, we reliably identify the corresponding spin and parities for all ground and excited charmonia with masses m ≤ 4.0 GeV, J ≤ 3 and p 2 ≤ 2 (2π/L) 2 .
It is straightforward to apply the procedure illustrated in this work to study the excited hadron spectrum on any other lattice QCD ensemble. The relevance of different interpolators could change depending on the system (light mesons, heavy mesons or heavy quarkonium) under investigation. Naively one may have to repeat the whole analysis for each system being studied and on each ensemble. Even if this is the case, for heavy quarkonium spectroscopy qualitative information on the opera-tor state overlaps from this investigation can alleviate the efforts to extract the spin identified spectrum on a different ensemble 7 . Note that the energy ordering of near degenerate states can vary depending on the parameters of the ensemble. Hence one cannot simply carry over the J P assignments from one ensemble to the next based on energy ordering alone.
The results in this paper will be valuable for our ongoing study of charmonium resonances in (coupled-channel) meson-meson scattering, where we combine the quarkantiquark basis used in this study with meson-meson interpolators to go beyond the single hadron approach.

ACKNOWLEDGMENTS
We would like to thank G. Bali  We are grateful to the Mainz Institute for Theoretical Physics (MITP) for its hospitality and its partial support during the completion of this work. We thank our colleagues in CLS for the joint effort in the generation of the gauge field ensembles which form a basis for the here described computation. The simulations were performed on the Regensburg iDataCool and Athene2 clusters and the SFB/TRR 55 QPACE2 [29] and QPACE3 machines. Part of the simulations were preformed at the local cluster at the Department of Theoretical physics at Jozef Stefan Institute. The CHROMA [30] software package was used extensively along with the multigrid solver implementation of Refs. [31][32][33][34] (see also Ref. [35]).    12. |λ|-identified charmonium spectrum in the moving frame with p = (1, 1, 0). All irreps Λ C of the corresponding little group Dic2 are presented (Table I)  FIG. 14. Summary of the charmonium masses extracted from our lattice simulation for various J P C within the single-hadron approach (using only quark-antiquark interpolators) (magenta). Mass splittings from the mass of the ηc are shown together with the experimental values [26] (black). We emphasize that these results are a qualitative illustration from a single lattice ensemble at unphysical pion mass, and likely do not capture mesons close to decay thresholds or broad resonances correctly. The candidate for the first excited scalar charmonia from experiment is not yet settled and we show the mass of the recently discovered X(3860) [27]. . Shown are all states with reliably-identified spin J and parity P , which are found at rest and at least for one non-zero momentum frame.