{'fraction_non_alphanumeric': 0.03737760749255002, 'fraction_numerical': 0.01772669220945083, 'mean_word_length': 4.280640230195126, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 8, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 19, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Exotic and excited-state radiative transitions in charmonium from lattice QCD

We compute, for the first time using lattice QCD methods, radiative transition rates involving excited charmonium states, states of high spin and exotics. Utilizing a large basis of interpolating fields we are able to project out various excited state contributions to three-point correlators computed on quenched anisotropic lattices. In the first lattice QCD calculation of the exotic 1-+ eta_c1 radiative decay, we find a large partial width Gamma(eta_c1 ->J/psi gamma) ~ 100 keV. We find clear signals for electric dipole and magnetic quadrupole transition form factors in chi_c2 ->J/psi gamma, calculated for the first time in this framework, and study transitions involving excited psi and chi_c1,2 states. We calculate hindered magnetic dipole transition widths without the sensitivity to assumptions made in model studies and find statistically significant signals, including a non-exotic vector hybrid candidate Y_hyb? ->eta_c gamma. As well as comparison to experimental data, we discuss in some detail the phenomenology suggested by our results and the extent to which it mirrors that of quark potential models and make suggestions for the interpretation of our results involving exotic quantum numbered states.


I. INTRODUCTION
The charmonium system is often described as the "hydrogen atom" of meson spectroscopy. Being reasonably non-relativistic, it is explained fairly well by potential models, at least below the open-charm (DD) threshold. More recently it has also been studied using effective field theory approaches (such as pNRQCD) and QCD sum rules. There has lately been a resurgence of interest in the charmonium system with the B-factories, CLEO-c and BES finding missing states, making more accurate measurements of properties of these states and discovering a number of new resonances that are not easily explained by the quark model. This has spurred renewed theoretical interest with much speculation as to whether these states are hybrids or multiquark/molecular mesons. To date there are no charmonium states having manifestly exotic J P C such as 1 −+ , 0 +− , 2 +− , that would directly signal physics not present in potential models.
The states below open-charm (DD) threshold can not decay via an OZI allowed strong decay and so have reasonably narrow widths. Their radiative transitions can therefore have significant branching ratios and are experimentally accessible. The transitions from and production of the, as yet unobserved, exotic 1 −+ are particularly interesting. A lattice calculation of transition form factors of excited charmonia is therefore timely and this is the first such study. The corresponding excited charmonium spectrum was calculated in lattice QCD in Ref. [1]. Transition form factors of the lightest few charmonia, those ground states accessible with interpolating fieldsψΓψ, were calculated in Ref. [2]; this work brought to the attention of CLEO-c experimentalists the discrepancy between the lattice calculated value of Γ(J/ψ → η c γ) (and indeed the values predicted in most model calculations) and the single experimental measurement of this from Crystal Ball [3]. In a tour-de-force analysis [4], a much more reliable value was extracted from CLEO-c data that is in much better agreement with theoretical estimates.
The calculation we will present is performed in the quenched approximation, neglecting altogether the effect of light quark degrees of freedom. As such it is rather directly related to the simplest quark potential models in which charm quarks move in a static potential of assumed gluonic origin. Attempts have been made to add in the effects of light-quark loops to these models [5,6,7], in some cases finding that these effects can be large [8]. We will address this possibility in light of our results.
A strong motivation for developing the lattice QCD techniques required to extract excited and exotic state radiative transition matrix elements is the versatility of the method. We can use these methods, tested here in charmonium, at any computationally feasible quark mass. In particular this opens up the possibility of computing in a framework close to QCD the meson photocouplings that appear in the meson photo-production process to be utilized in the JLab 12 GeV GlueX experiment [9]. In this paper we will also compare results with the flux-tube model of gluonic excitations which to date is the only theoretical guide to the size of the hybrid couplings and hence production rates [10,11].
The paper is structured as follows: We begin in section II with a description of the technology used to construct three-point correlators and to project on to the contribution due to various excited states. In section III we present our results for the transition formfactors between various meson states. A discussion of the phenomenology of these results in terms of experiment and models follows in section IV before we conclude in section V. Two appendices which consider some technical details complete the manuscript.

II. TECHNOLOGY
In this paper we will explore radiative transitions with charm-mass quarks using the quenched anisotropic lattices and clover fermion action described in [1]. Radiative transition matrix-elements follow from vector current three-point functions, whose construction using sequential source technology is described in [2]. As in [2], only connected diagrams are considered with the assumption that disconnected diagrams are negligible in comparison. The new item of technology in this study is the use of a large basis of meson interpolating fields as explored in two-point functions in [1] with further interpretation in [12].
These operators are lattice-discretised versions of the gauge-invariant fermion bilinears, D is the gauge-covariant derivative operator. In order to improve the overlap on to lower lying mesons, the quark fields may also be smeared over space as e 1 4 σ 2 D i D i ψ(x). Linear combinations of these basic operators have been constructed which transform irreducibly under the lattice cubic symmetry and which have rather simple interpretations in the limit of zero lattice spacing, for details see [1,12]. The operator set is large enough to have a considerable redundancy within a given quantum number sector, which can be utilized to extract excited states.

A. "Ideal" operators and eigenvector projection
In [1] the spectrum of charmonium was extracted by solution of the generalised eigenvalue problem 1 which constitutes a best solution for the spectrum in the variational sense. The matrix of two-point correlators is constructed as C ij (t) = 0| x O i ( x, t) O j ( 0, 0)|0 where the sum over lattice sites forces the three-momentum of single particle states to be zero 2 . This method relies upon a redundancy of operators O i in any given quantum number sector. Fitting the time-dependence of the principal eigenvalues, λ n (t), gives us the masses (or energies at finite three-momentum) of the states with the quantum numbers of the hermitian operators O i,j . Details of the two-point analysis, including the importance of the choice of timeslice t 0 , are given in [1]. The interpretation of the eigenvector v n in the solution of Eqn. 1 can be expressed as follows: weighting the operators by this vector gives the optimal operator within the limited operator space for the state labeled by n. It is convenient to normalise these "ideal" operators in a manner which accounts for the value of t 0 : Ω n = √ 2E n e −En t 0 /2 v n i O i . Note that these eigenvectors are trivially related to the vacuum-operator-state matrix elements or "overlaps", Z n i = n|O i |0 , that appear in a bound-state spectral decomposition, The procedure we shall follow in this paper is to compute three-point correlators having at the source (located on a fixed timeslice, t i ) a (smeared) local operator of the form ψ( 0, t i )Γψ( 0, t i ) -all possible gamma matrices can be considered for the computing cost of a single "forward" propagator. In this work we consider Γ = γ 5 , γ i , 1 giving access to pseudoscalar, vector and scalar states at the source.
At the sink (located on a fixed timeslice, t f ) we use sequential-source technology using a broad selection of local and derivative-based operators as described in [1] -for each operator at a given momentum (usually p f = 000) we have the computing cost of a single "backward" propagator. At the vector current insertion (inserted on all timeslices t i < t < t f ) we insert the local vector current with all possible lattice three-momenta q up to | q| 2 = 4. Translational invariance ensures momentum conservation, selecting the correct value of p i out of the sum over all momenta produced by a local operator.
1 repeated indices are summed 2 In some cases we also considered non-zero three-momentum correlators -this will be discussed later Three-point correlators for "ideal" operators are constructed from the "raw" correlators, eqn 3 by projecting with the appropriate eigenvector 3 , The net effect is to provide us with a three-point correlator that, while it has multiple states contributing through the local source operatorψ( 0, t i )Γψ( 0, t i ), should have only a single state, n, contributing from the sink operator Ω n . In principle we could improve this further by using the full basis of operators at the source as well as the sink, but using sequential source technology this quickly becomes rather expensive 4 .
In fact we can be a little more precise about the statement that only a single state contributes at the sink -in considering the analysis of the two-point correlators we determined that t 0 can be considered to be the timeslice on which the correlator matrix (of dimension N ) is saturated by N states -closer to the operator we require more states to describe the correlator matrix [1,14]. It follows that if we are t 0 timeslices away from the sink operator position t f we can be quite confident that we have contribution only from the single state labeled by n. When we come to fit the time-dependence of the three-point correlator we will not use timeslices any closer to the sink than t = t f − t 0 .
In figure 1 we show the effect of the eigenvector projection on to "ideal" states in the case of a pseudoscalar source operatorψγ 5 ψ and a selection of pseudoscalar sink operators. We clearly see that the projection on to the ground state "ideal" operator produces a considerable "flattening" of the correlator moving toward the sink. We do see some curvature beginning within t 0 = 8 timeslices away from the sink as expected 5 . More importantly for this study we see that projection on to the "ideal" operator for the first excited state yields a clear non-zero signal from which one can extract a transition matrix element.
There is in fact a slight technical subtlety to be dealt with in the eigenvector projection on to "ideal" operators that arises from the method of solution of the two-point problem. Solving the eigenvalue problem gives eigenvectors on each timeslice, v n j (t 2pt ) -if the solution is to be a true spectral representation there should be no time dependence in the eigenvectors, at least for t 2pt > t 0 where the correlator matrix is supposed to be saturated. In [1] we saw that over a considerable range of t 2pt , the Z values (obtained trivially from the v, see equation 2) were flat. We explore the effect of any such possible time-dependence on the three-point functions by performing the projection for eigenvectors belonging to a range of timeslices, v n j (t 2pt ). In figure 2 we show for a ground state and first excited state projection the t 2pt dependence on the projected correlator. It is clear that there is no strong dependence on t 2pt and we choose to average over the projections for various t 2pt to reduce configuration-by-configuration fluctuations that may be caused by the generalized eigensystem solver.

B. Transition form-factors
The "ideal" operator-projected three-point functions are related to the transition matrix elements between states as follows: While all possible eigenstates with the quantum numbers of O i will contribute at the source (the sum over m), the eigenvector projection ensures that we need only consider the single state n at the sink. Our normalisation of Ω n ensures that if f is spin-zero, The Minkowski-space transition matrix element f n ( p f ) j µ ( 0, 0) i m ( p i ) can be decomposed in terms of multipole form-factors multiplied by Lorentz covariant combinations of the momenta and (if appropriate) polarization tensors of the particles labeled by i, f : The general technique for obtaining these covariant multipole decompositions is given in the appendix of [2]. Inserting this decomposition into eqn. 4, and performing the implicit sums over helicity (in the case of spin ≥ 1), the three-point correlator on each timeslice, t, can be expressed as a linear sum of known "kinematic" and "propagation" factors times the unknown multipole form-factors 6 , In practice we opt to solve this linear system 7 on each timeslice, t, for a set of "effective", t-dependent form-factors,F kn (Q 2 , t), using the propagation and kinematic factors for the ground state at the source(m = 0): The propagation factor P 0,n ( p i , p f , t i , t, t f ) contains the overlap Z 0 Γ ( p i ) and the energy E 0 ( p i ). These are determined in the analysis of two-point function at finite momentum. Since the state masses at zero momentum are rather precisely determined by the variational solution, we use the continuum dispersion relation to obtain E( p) = m 2 + | p| 2 . Z values then follow from a linear fit to the finite-momentum two-point correlators supplying the known exp (−E n ( p)t) factors.
The extracted "effective" form factors have the property thatF k,n (Q 2 , t) ≈ F k,n (Q 2 ) + O e −(E m=1 −E m=0 )(t−t i ) so that away from the source the excited state contributions die away and the form-factor plateaus to the i m=0 → f n value. We fit the time-dependence of these extracted "effective" form-factors to a sum of exponentials where the energy dependence is that extracted from spectrum studies 8 . A typical example is shown in figure 3. We retain only the plateau value, discarding the relatively unreliable source excited state information. 9 The expansion of transition matrix elements as performed above relies upon the continuum Lorentz symmetry which is broken on our cubic lattice. A rigorous study would consider the expansion in terms of irreducible representations of the cubic group at finite momentum and would typically involve a larger number of form-factors, some of which would tend to the multipole form-factors as a → 0 and some of which would have to vanish. We have not performed such a decomposition but suspect that it is not important at this lattice spacing since when we solve the correlator-form-factor linear system assuming continuumlike decompositions we typically obtain χ 2 per d.o.f. very close to 1, suggesting that there is little need to enlarge the basis space.
The fact that charmonium states are eigenstates of charge conjugation (C), coupled with the photon having C = −1, means that all charmonium radiation transitions are between states of opposite C. The origin of this is that the photon couples equally to the charm quark and the charm anti-quark. In our lattice calculation we need not impose this symmetry; by coupling only to the quark we can obtain transitions between states of equal C. One way of viewing this is that it is like having a u, d pair but as heavy as the charm quark -this can useful for comparison with models and we will show results for a number of C-violating transitions. 9 Strictly speaking the kinematic factors for excited m states need not be proportional to the ground state values such that excited state multipoles may "leak" into the wrong ground state multipole -clearly this effect will fall off for t ≫ ti. An alternative approach is to include the excited state terms directly into the decomposition fit, expanding the space to also be over timeslices, the time-dependence being the discriminator for the various states contributing at the source -this guarantees the right kinematic factors for excited states. In all cases to be presented the analysis was cross-checked using this fitting method (the excited state energies obtained from applying the dispersion relation to the excited state masses reliably extracted using the variational method) and agreement within statistical fluctuations found.
level mass / MeV suggested state model assignment 0 spectrum extracted from two-point functions, suggested experimental state analogues and quark model bound-state assignments made in [12].

A. Vector current renormalization
Since we use the local vector current in the construction of our three-point functions we need to determine the renormalization constant, Z V , to relate the extracted matrix elements to physical matrix elements. We do this by insisting that the pseudoscalar form-factor at zero Q 2 takes the value 1. On an anisotropic lattice we should allow there to be different renormalization constants for temporally and spatially directed currents, indeed we find Z V (s) = 1.23 (2) and Z V (t) = 1.118 (6). In all results presented below only the spatially directed current is used. See appendix B for a discussion of the effect of improvement of the vector current.

B. Scalar -Vector transitions
The first results we will present concern transitions between scalar (0 ++ ) and vector (1 −− ) states of charmonium. Using a quark-smeared operatorψψ at the source and 13 vector operators 10 at the sink we extracted transitions between the ground state scalar χ c0 and the lowest six vector states. In table I we show the spectrum reported in [12], possible comparable experimental states and a model-dependent state assignment. This transition is characterized by two multipole amplitudes, a transverse electric dipole E 1 (Q 2 ) and a longitudinal C 1 (Q 2 ), the matrix element decomposition being where Ω( The "vector" operators used are, in fact, all in a particular irreducible representation of the cubic group at rest, namely T 1 which in the continuum limit contains spins 1, 3, 4 . . . . By (12) 1.14(5) 199 (6) 131(14) 1 considering the degeneracy pattern and the values of the overlaps Z (n) j , the two-point function analysis strongly suggests that the 2nd excited state in the T −− 1 sector is in fact a spin-3 state and as such we do not present results involving this state here [1,12]. Additionally we do not show the amplitudes for the 4th excited state 11 as the signals are statistically consistent with zero for all Q 2 . Shown in figure 4 are the electric dipole transition form-factors for the four relevant states. See appendix A for a discussion of the inclusion of p f = (100) correlators in this analysis.
In each case the Q 2 dependence of the lattice data has been fitted with the form, whose motivation is described in [2]. The results of the fits are given in Table II. One point to make here is that in our lattice calculations, the photon only couples to the quark and not the antiquark and we do not explicitly include the electric charge of the quark. Therefore we actually compute and presentF k (Q 2 ) which are related to the physical multipole amplitudes by F k (Q 2 ) = 2 × 2 3 e ×F k (Q 2 ). Note also that the multipole amplitudes have mass dimension 1 and hence we plot them in temporal lattice units, where a −1 t = 6.05 GeV is determined from the static potential. Partial decay widths follow by averaging over initial helicities and summing over final helicities which gives, for a general A → Bγ decay: Where we have experimental masses for states we use these to compute the phase space, otherwise we use the value extracted from the lattice calculation [1,12]. The ground state transition χ c0 → J/ψγ form-factor shows behavior rather similar to that found in [2] which used the same lattices but Domain Wall fermions rather than Clover and which did not make use of multiple sink-operators and operator projection. We will consider this in a little more detail in Appendix B where we will consider the effect of O(ma) improvement of the local vector current. We will discuss the results in comparison with experiment and with quark potential models in section IV B.
FIG. 4: Electric dipole transition form-factors χ c0 ↔ ψ. Plotted is the form-factor in temporal lattice units against the photon virtuality in GeV 2 . Fits to the lattice Q 2 dependence as described in the text. Experimental points at Q 2 = 0 are extracted from experimental decay widths taken from [15,16]. Results of fit to lattice data using equation 6. Partial decay width computed using fitted value ofV (0) and physical phase space (where known). All errors are purely lattice statistical. Experimental partial decay widths from [4,15].

C. Vector -Pseudoscalar transitions
Using the same T −− 1 operator set at the sink and the quark-smearedψγ 5 ψ operator at the source we obtained results for the single magnetic dipole form-factor in the vectorpseudoscalar transition. The decomposition used 12 is  table III. We refer the reader to the paper [2] for a discussion of the systematic error introduced into the phase space by an inaccurate lattice estimate of the hyperfine splitting.

D. Axial -Vector transitions
In transitions between axial (1 ++ ) and vector (1 −− ) states there are two transverse formfactors, electric dipole (E 1 ) and magnetic quadrupole (M 2 ), and one longitudinal form-factor (C 1 ). The matrix element decomposition takes the form  In this case 13 we computed with the vector at the source, produced using quark-smeared ψγ j ψ, and a set of ten axial (actually T ++ 1 ) operators at the sink 14 . In figure 6(a) we present the two transverse form-factors for the ground state axial meson and in figure 6(b) their ratio. The amplitudes were each fitted with a form like eqn. 6 which unfortunately is not as well constrained as in the scalar case owing to the kinematic factors preventing a slightly time-like Q 2 point corresponding to p i = p f = (000). The results of the fits are a tÊ1 (0) = 0.23(3), β E1 = 440(40) MeV, λ E1 = 0.71(30) GeV −2 and a tM2 (0) = −0.020 (17), β M 2 = 450(50) MeV, λ M 2 = 5(6) GeV −2 . This corresponds to a partial decay width of Γ(χ c1 → J/ψγ) = 270(70) keV which is in reasonable agreement with the PDG's average of 320(25) keV.
The ratio M 2 E 1 (Q 2 ), shown in figure 6(b), was fitted with various functional forms shown 13 We could have usedψγ 5 γ j ψ at the source to produce the axial meson, this operator at finite momentum also has overlap with the pseudoscalar state which severely limits its usefulness. 14 γiγ5, ρ × ∇T 1, ρ (2) × ∇T 1, a1 × DT 1, b1 × BT 1 in both quark-smeared and unsmeared versions. by the shaded bands yielding M 2 E 1 (0) = −0.20 (6) where the error includes a crude estimate of the systematic error due to the uncertainty in fitting form. The ratio of the extrapolated values from the separate fits to (7). Clearly without data points at smaller Q 2 or some certainty about the expected Q 2 dependence, we cannot constrain this any further and hence cannot make a particularly meaningful comparison with the PDG average .008 −0.017 . Form-factors for the transition from the first excited axial state, χ ′ c1 down to the J/ψ are shown in figure 7 where multiple fit forms were used, all returning a χ 2 /N dof close to one. The estimates for the physical photon point thus obtained are a tÊ1 (0) = 0.050 (15) and a tM2 (0) = −0.004 (14) where again we include a crude systematic error estimate for the fit-form variation. The E 1 transition corresponds, for a χ ′ c1 at 4.1 GeV, to a partial decay width Γ(χ ′ c1 → J/ψγ) = 21 (12) keV.

E. Tensor -Vector transition
The transition between the lightest 2 ++ state and the lightest vector state was not considered in [2] as that study used only local fermion bilinears (ψ(x)Γψ(x)) to produce states -a spin-2 particle cannot be produced by any such operator. Here we use a set of six operators projected into T ++ Here In order to have a constrained linear system for extraction of the five form-factors we are required to consider all five spin-2 helicities, which with respect to the lattice cubic symmetry are distributed in T 2 (3) and E(2) irreps. Strictly speaking these are independent irreps on the lattice and we should be careful about combining them, but any lattice symmetry breaking must vanish as a → 0 and in [1] we found that there were strong signals that the continuum rotational symmetry was restored to a good approximation already at this value of a. For example, for the ground states in the T 2 , E channels we have a high degree of mass degeneracy, and the values of the overlaps (Z T 2 /Z E ) are compatible at the 1% level. Details of this analysis are in [1].
The extracted form-factors and fits to the Q 2 dependence are shown in figure 8 where the expected hierarchy |E 1 (0)| > |M 2 (0)| ≫ |E 3 (0)|, to be discussed in section IV B, is observed. The description of the lattice three-point correlators by the continuum decomposition, eqn 7 (which does not take any account of T 2 /E discretization differences) is excellent, with typically χ 2 /N dof ∼ 1. As a simple systematic test of the degree to which T 2 /E discretization differences could affect the determination of the form-factors we tried deliberately introducing such a difference by arbitrarily multiplying all E correlators by a factor. Even with a 10% increase in the magnitude of all E correlators, the form-factor values changed by less than the statistical error, while the χ 2 /N dof increased to around 6. If we where to see such large χ 2 in our raw data (which we never do) we would not trust the results -we conclude that our form-factor extractions are not being strongly affected by T 2 /E discretization differences.
The amplitudes at Q 2 = 0 correspond to a partial decay width of Γ(χ c2 → J/ψγ) = 380(50) keV to be compared with the PDG average of 406 (31) (7) and (11) which in the case of the first ratio is considerably larger than the PDG average of −0.13(5) but which does appear to be of the correct sign.
We also extracted the χ ′ c2 → J/ψγ transition form-factors 16 as presented in figure 9(a). The large value of E 3 might be surprising when compared to the ground state result -a number of simple systematic tests were performed to investigate if this might come about through various lattice effects with the result that we were unable to change the values outside statistical error bars by any reasonable adjustment. The amplitudes at Q 2 = 0 correspond to a partial decay width of 20(13) keV. In figure 9(b) we present the transition form-factors for the next excited tensor state 17 , χ ′′ c2 → J/ψγ where the hierarchy |E 1 (0)| ≫ |M 2 (0)|, |E 3 (0)| appears to be restored and where we predict a partial decay width of 88(13) keV. In section IV B we will propose a simple explanation for all these observations in the framework of a non-relativistic quark model.

F. Exotic transitions
Our principal focus here is the exotic 1 −+ state, η c1 , here 18 found to be at 4300(50) MeV. On a quenched lattice, provided we can eliminate the possibility that the lightest state in T −+ 1 is part of a non-exotic 4 −+ (see [12] for support of this elimination), we can be fairly certain 16 We found a mass of 4115(28) MeV for the χ ′ c2 17 We found a mass of 4165(30) MeV for the χ ′′ that this state is a hybrid, having an excited gluonic field in addition to a charm-anticharm quark pair. This is in contrast to having a higher quark number Fock state which, since we lack light-quarks altogether in this calculation, could only arise for states having mass near 4m c 19 . A charge-conjugation allowed decay of this meson would be η c1 → J/ψγ having transverse magnetic dipole and electric quadrupole multipole contributions. The decomposition of the vector current matrix element between two non-identical vector mesons takes the following form: . We used a set of eight operators at the sink to produce the T −+ 1 state 20 . Our results projected on to the ground state in T −+ 1 are shown in figure 10 -the lattice data are not fitted as a function of Q 2 as very little extrapolation is required to associate the point at Q 2 = 0.06 GeV with the real photon point. The value of M 1 (0) corresponds to a partial decay width Γ(η c1 → J/ψγ) = 115(16) keV. Note that this is no different in scale to many measured conventional charmonium transitions. If we also consider charge-conjugation symmetry violating decays (by coupling only to one quark) we have access to 1 −(+) → 0 −(+) γ, 1 −(+) → 0 +(+) γ. These results can be compared to models of hybrid bound-state structure as we will do in the following section. Results are shown in figure 11.
Another exotic J P C to which we have easy access is the 0 +− state which appears in the A +− 1 channel at 4465(65) MeV. We have computed a charge-conjugation symmetry violating decay of the lightest state with these quantum numbers to the conventional ground state vector 1 −(−) . This transition matrix element has the same decomposition as the vectorscalar in equation 5. We used a set of three operators 21 to create the 0 +(−) state; in figure  12 we show the electric dipole form-factor for 0 +(−) → 1 −(−) γ.

G. Other assorted transitions
We also computed some C-violating transitions which have allowed analogues in the light quark sector 22 , 1 +(+) → 0 −(+) γ and 2 +(+) → 0 −(+) γ. The first of these has a multipole decomposition identical to eqn 5 while the second takes the form Results are shown in figure 13. A discussion in terms of non-relativistic quark models follows in the next section.

IV. PHENOMENOLOGY
The results presented in the previous section can be compared with experimental measurements and additionally with the predictions of various models describing heavy quark bound states. A rather successful approach to describing charmonium has come by considering heavy charm quarks to be moving non-relativistically (or nearly so) in a static potential motivated by QCD. These quark potential models have been used to compute spectra, radiative transition rates, and, in many-body extensions of the theory, hadronic decay rates (e.g. [6,18,19]). In their simplest formulation, the gluonic field plays no dynamical role and one does not have exotic quantum numbered states having excited glue (hybrids). In order to consider such states one can construct models making specific assumptions about the nature of the gluonic field and its excitations, examples being the flux-tube model [20], "constituent" gluon models and models based upon a many-body treatment of QCD in a physical gauge [21]. In what follows we will address to what extent our lattice QCD results inform these models.
Firstly we should discuss the degree to which the approximations we have made in our computation introduce systematic error to our results. 22 namely a ± 1 → π ± γ and a ± 2 → π ± γ.

A. Lattice Systematics
This calculation has been performed within the quenched truncation -as such one way to view it is as a calculation of a version of QCD having just one heavy flavor of quark with the neglect of heavy flavor quark loops being justified by their large mass. In this regard this calculation is rather directly comparable to quark potential models which in their simplest form also neglect the effect of dynamical light quarks. We set the lattice scale in our calculation using the Sommer parameter which is itself related to the static quark potential and the charmonium spectrum, so here too the comparison is fairly direct. In figure 14 we show the static potential extracted on this lattice along with the phenomenological forms used in a selection of potential models.
Clearly then our calculation does not include the effect of charmonium states coupling to multi-hadron states containing light quarks, e.g. DD. There have been suggestions within extensions of quark-potential models that including such physics can have a considerable effect on the spectrum [5,6] and radiative transition rates [8] of charmonium. We will attempt to address this in the discussion to follow.
In terms of direct comparison with reality, in which there are three flavors of quark lighter than charm, another failing introduced by the quenched approximation is the incorrect run-FIG. 14: Lattice static potential along with phenomenological potentials used in [7,18,22]. ning of the strong coupling owing to the beta function containing the wrong number of quark degrees of freedom. We expect this to show up in terms of inconsistent scale settingsomething that has been observed previously [23].
Our calculation has been performed at only one value of the lattice spacing; as is well known one expects there to be systematic shifts in quantities owing to the lattice discretization that go away as one approaches the continuum. These are suggested to be particularly serious for the discretization of heavy quarks. At a single lattice spacing we cannot accurately estimate the size of such effects, although our use of the improved Clover action and an anisotropic lattice 23 should help a great deal in the reduction. We present in Appendix B a limited study of improving the vector current insertion in the manner described in [24]. In general the improvement effects are not large, suggesting that scaling to the continuum should not overwrite our results.
In the spectrum calculations using the same action on the same lattices [1,12], the possibility that the spatial volume is too small to comfortably house highly excited states was raised.
One way to reduce the systematic error introduced by extrapolating from finite Q 2 to Q 2 = 0 would be to consider utilizing twisted boundary conditions on the propagator inversions to get momenta and hence Q 2 values rather close to Q 2 = 0 -this has proven successful in form-factor computations [25].
Modulo the caveats that we have raised, we remain convinced that the results we present likely represent a faithful description of the pattern of physics of radiative transitions in charmonium, including the properties of gluonic hybrids. Most importantly we have clearly demonstrated that the technology of projecting three-point functions using "ideal" operator eigenvectors obtained in two-point function calculations works well in giving us access to excited state transitions. All the possible sources of systematic error can be addressed in future calculations using a set of sufficiently large dynamical lattices of various lattice spacings.

B. Conventional state transitions and the quark-potential model
Ref. [12] used information extracted from two-point functions to identify excited charmonium states found in a lattice calculation with the corresponding states expected in the quark model. The transition form-factors can also be compared to calculations in quark-potential models. Our quenched lattice calculations offer a rather direct comparison with these potential models: in both cases loops of mesons containing light quarks have been ignored. As discussed in Ref. [2,26], the results of quark models are sensitive to approximations, such as the choice of frame, which are not an issue in lattice calculations -we will not discuss these further here. In the following, for simplicity, we use harmonic oscillator wavefunctions and consider the rest frame of the initial meson.
Within the simplest non-relativistic model of a meson emitting a photon of momentum q, we find the following transition form factors for non-radially excited states undergoing a change of orbital angular momentum of one unit: where β is the harmonic oscillator wavefunction parameter. Here the r-values follow from the coupling of quark spin and orbital angular momentum as 24 Transforming from | q| 2 to the invariant virtuality of the photon (Q 2 ) in the rest frame of the decaying meson this corresponds to Note that this is of the general form given in eqn. 6. In the same model, we have for transitions involving both a change in orbital angular momentum and a quark spin-flip, the expression This applies to the charge-conjugation forbidden decays 1 In addition, this same form is predicted for the M 2 formfactor in the transitions χ c1,2 ( 3 P 1,2 ) → J/ψ( 3 S 1 )γ and the charge-conjugation forbidden In terms of Q 2 this gives which we note is also of the form of eqn 6. 24 note that a typographical error in [2] is corrected here Within this model we find that E 3 (Q 2 ) = 0 for χ c2 ( 3 P 2 ) → J/ψ( 3 S 1 )γ. In fact this result follows from the general assumption of a single quark transition [27] and does not depend critically upon the details of the model. For a transition J → J ′ (where J and J ′ are respectively the initial and final meson spins) the allowed k (for E k or M k depending on k and the parity change) are in general |J − J ′ | ≤ k ≤ J + J ′ . However, if only a single quark is involved in the transition there is a further restriction on the values allowed. If the interacting quark in the initial (final) meson has total angular momentum j (j ′ ) then |j − j ′ | ≤ k ≤ j + j ′ . Hence in general for an n 3 P 2 → n ′ 3 S 1 transition the allowed k are k = 1, 2, 3. However in the single quark transition assumption we have j = 3/2 → j ′ = 1/2 and k is restricted to k = 1, 2. Note that it is possible to have E 3 = 0 if the tensor meson is a 3 F 2 (j ′ = 5/2) state, in which case k = 2, 3 and E 1 = 0. It is also possible to have E 3 = 0 in transitions involving hybrid, multiquark or molecular mesons where there are other degrees-of-freedom able to carry some angular momentum.

ψ − η c transitions
Our result for Γ(J/ψ → η c γ) is in reasonable agreement with the new experimental value from CLEO-c [4], the potential model calculations of [18,28], the pNRQCD result of [29] and the QCD sum-rule estimate of [30]. Similarly our result for the "hindered" M 1 transition Γ(ψ ′ → η c γ) is consistent with the experimental values [4,15] and Eichten et. al. [28], but significantly smaller than the values in Ref. [18] which depend sensitively on the overlap of orthogonal wavefunctions with a small recoil factor breaking the orthogonality.
Li and Zhao [8] study the possibility of hadronic loop contributions in J/ψ and ψ ′ radiative decays to η c and η ′ c . They find the addition of loop contributions involving pairs of D and D * mesons can bring the model values from Ref. [18] in to line with experiment. Most significantly, they find large loop contributions, of order 10 keV, to the ψ ′ → η c γ transition. A cancellation between the quark model contribution of order 10 keV (from Ref. [18]) and these loop contributions results in a small transition width of ∼ 1 keV. This appears to present a problem with respect to our result, which being quenched does not include loop contributions, since there is no room for a large loop contribution when our Γ(ψ ′ → η c γ) is compared to the experimental value. This might indicate that the loop contributions of Ref. [8] are overestimated and we have identified one possible source for this in the coupling g D * Dγ used. The partial decay width Γ(D * 0 → D 0 γ) ∼ 800 keV assumed in Ref. [8] appears to be rather large for an M 1 radiative transition and, with the branching ratio from the PDG [15], implies a D * 0 total width of 2.1 MeV which is at the upper limit allowed by the PDG. For comparison, Close and Swanson [31], using a non-relativistic quark model predict Γ(D * 0 → D 0 γ) = 32 keV. The corresponding decrease in the loop amplitudes would remove the problem in comparison with our result. In addition, there are considerable uncertainties arising from the estimated hadronic coupling constants in the loop contributions which are not discussed in Ref. [8]. These uncertainties may be particularly important because of the delicate nature of the cancellation.
We have also been able to extract a signal for the ψ ′′ → η c γ which in the standard interpretation of the ψ ′′ would be a 1 3 D 1 → 1 S 0 transition. In a harmonic oscillator wavefunction basis the form-factor has the same leading | q| 3 behaviour as 2 3 S 1 → 1 S 0 γ indicating that it should suffer "hindered" suppression just like the ψ ′ → η c γ relative to J/ψ → η c γ. The lattice data for the amplitude appears to be in line with this, within a large statistical uncertainty. Discussion of the next excited vector state will follow in section IV C.

χ-ψ transitions
The χ cJ → J/ψγ partial widths we obtain are dominated by the electric dipole component and compare favorably with the experimental values. The largest discrepancy is in the statistically precise χ c0 transition. The lattice systematic errors discussed above may be to blame for this, without further calculation we cannot determine this definitively, but we note that [7] found that approximately "unquenching" the quark model for this transition did not induce a large change in the rate. In Appendix B we consider the effect of using an improved vector current and we find that the χ c0 → J/ψγ transition is one case where a statistically significant change with respect to the local current does occur. This might indicate that this transition is particular sensitive to scaling in the lattice spacing.
The fitted values of β in eqn 6 are in rough agreement for all of χ c0,1,2 as one might expect for states that in the potential models differ only by small spin-orbit effects. The fitted values of λ deviate somewhat from the quark-model expected 1 : 1 2 : − 1 2 ratio for χ c0,1,2 , being 1 : 0.62(26) : −0.34 (2). Given the approximations and frame-dependence inherent in the quark model even such rough agreement is surprising.
In the χ c1,2 → J/ψγ transitions, the hierarchy of multipoles expected in the quark model (|E 1 (0)| > |M 2 (0)| ≫ |E 3 (0)|) is observed. The precise degree of suppression of M 2 compared with E 1 is not easy to compute in non-relativistic quark potential models. However, the single quark transition assumption predicts that E 3 (Q 2 ) = 0 independent of the frame and details of the potential and the lattice data do agree with this. Our extracted values of M 2 (0) or the ratio M 2 (0)/E 1 (0) all depend upon theoretically undetermined extrapolation in Q 2 and so the apparent disagreement with the very small values found in experiment is not yet overly concerning.
The electric dipole couplings for ψ ′ → χ c0 γ and ψ ′′ → χ c0 γ are, within reasonably large statistical errors, in agreement with experimental values and with the "GI" values tabulated in [18] 25 . The estimates of D-meson loop effects from [7] are comparable to the level of our statistical uncertainty.
For the transition χ ′ c1 → J/ψγ we extract a partial width of 21 (12) keV which is in reasonable agreement with the quark-model estimates of [18], being 14 or 71 keV depending upon model details. Since in our study this state has a likely interpretation as the 2 3 P 1 state of the quark model, without any DD * effects included, this can act as a benchmark for models that consider the experimental X(3872) state as being 1 ++ .
The large E 3 (0) in the lattice calculation of the χ ′ c2 → J/ψγ transition is in stark disagreement with the general result that E 3 = 0 for a 3 P 2 → 3 S 1 transition. As discussed above, this does not appear to be a lattice artifact and so another interpretation must be found. Possible explanations include that the χ ′ c2 state is a hybrid where the gluonic field carries spin, but here we concentrate on the most conservative interpretation, that χ ′ c2 is a conventional 3 F 2 state.
An F -wave tensor state ( 3 F 2 ) is expected to have E 1 (Q 2 ) = 0 in general from the single quark transition assumptions outlined earlier. In addition, a quark model calculation shows that M 2 (0) is highly suppressed: whereas E 3 is less suppressed: with the same leading | q| 2 dependence as the M 2 form factor in a 3 P 2 → 3 S 1 transition (Eqn. 11). This pattern that |E 3 (0)| > |M 2 (0)| ≫ |E 1 (0)| is consistent with the lattice results for this χ ′ c2 → J/ψγ transition. These lattice results support the interpretation of this χ ′ c2 as the lightest 3 F 2 state having a partial decay width to J/ψγ of 20 (13) keV. We are unaware of any model calculations of this transition rate.
The χ ′′ c2 → J/ψγ transition appears to have reverted to the |E 1 (0)| > |M 2 (0)| ≫ |E 3 (0)| hierarchy expected for a 3 P 2 state and we propose that this state is the 2 3 P 2 . There is of course the possibility that the P and F states are mixed, but this mixing cannot be anything like maximal given a large E 3 appears in only one case and a large E 1 only in the other case. Within a non-relativistic quark model, the 3 P 2 and 3 F 2 states can mix via any tensor potential in the Hamiltonian. However in most models this term is suppressed by the inverse charm mass squared and so is small. The tensor term, like the hyperfine interaction, is a short distance effect and so may not be reliable in these quenched lattice calculations. Because this term is suppressed anyway, the errors introduced from this should be small. Another possible source of mixing which we lack in the quenched theory comes through DD meson loopssuch effects are discussed in [5,32].
The partial decay width Γ(χ ′′ c2 → J/ψγ) extracted from the lattice data is 88(13) keV which is comparable with the 53 or 81 keV computed by [18]. We note that the interpretation of excited states in this channel as F and P waves was hinted at in the two-point function analysis of [12] where states compatible with being 3 F 2 and 2 3 P 2 were found to be nearly degenerate. Within the quark potential models used in [18], the F wave state is expected 60 − 100 MeV heavier than the P wave. The degeneracy we found may be an artifact of "squeezing" these spatially larger states into a 1.2 fm box.
A discussion of the excited χ c2 states is timely given the observation of a candidate state at 3929(5) MeV by Belle in γγ → DD which has been associated with the 2 3 P 2 state [33]. We lack reliable theoretical estimates of the γγ widths of 3 F 2 and 2 3 P 2 states -this calls for an extension of the work done in [34] to consider excited states which could use very similar operator projection technology to that discussed in this paper. We note here that an optimist viewing figure 2 of [33] might hope that the statistically insignificant excess at 4080 MeV could, with increased statistics, become a signal for the other state in the 2 3 P 2 / 1 3 F 2 pair.

C. Exotic & crypto-exotic transitions
The only charge-conjugation allowed transition involving an exotic we compute, η c1 → J/ψγ, has a rather large partial decay width 115 (16) MeV which is dominantly through a magnetic dipole transition. Even accounting for the large phase space this is very large on the usual scale of magnetic dipole transitions. Conventional cc states can only have magnetic dipole transitions if there is a quark spin-flip and this is suppressed by the large charm-quark mass. In a hybrid meson the extra gluonic degree-of-freedom allows an M 1 transition without spin-flip. This can be seen explicitly within the flux-tube model, where in an M 1 transition between a conventional L = 0 meson and a hybrid meson the tube absorbs the angular momentum 26 .
The two-point function analysis of [12] suggested that the excited vector state we have called Y hyb. (having mass around 4.4 GeV) is a crypto-exotic hybrid meson, having conventional J P C = 1 −− but internally an excited gluonic field. That analysis preferred the quarks to be in a spin-singlet in this state. In that case we would expect a large non spin-flip M 1 transition Y → η c γ which is precisely what is seen in our lattice data. Here Γ(Y → η c γ) = 42(18) keV which is considerably larger than any other vector to pseudoscalar transition. Within the flux-tube model, the non-exotic Y (1 −− hyb. ) and the exotic η c1 (1 −+ hyb ) (which we find at a mass ∼ 4.3 GeV) differ only in being quark spin singlet and triplet respectively and Γ(Y → η c γ) ≈ Γ(η c1 → J/ψγ). The lattice data does not strongly disagree with this. It would be very interesting to see analogous calculations in other models of excited glue to see if this is a general result or one peculiar to the flux-tube model.
Further evidence for the quark spin triplet nature of the η c1 (1 −+ hyb ) and the quark spin singlet nature of the Y (1 −− hyb. ) comes from the fact that E 1 (0) for 1 −(+) → 0 +(+) γ is 0.06(1) while for 1 −(−) hyb → 0 +(+) γ it is consistent with zero. The first of these then is S qq = 1 → S qq = 1 while the second requires a spin-flip, expected to be suppressed by the heavy-quark mass. Note that the magnetic dipole transition 1 −(+) → 0 −(+) γ in figure 11 is at the scale of other "hindered" M 1 transitions, consistent with requiring both gluonic excitation and quark spin-flip. We note that the spin singlet nature of our Y state does not make it a good candidate for the experimental Y (4260) state whose potentially large decay rate into ππJ/ψ would tend to suggest dominance of spin-triplet.
Here we briefly comment that if this heavy-quark hybrid physics is any guide to the behaviour of light-quark hybrid systems, we should expect the GlueX experiment to copiously photoproduce hybrid mesons off the meson cloud around a baryonic target.

V. CONCLUSIONS
We have demonstrated that the "ideal" excited state operators within a basis of operators can be used to successfully extract excited state transitions from three-point correlators. Using this technique we have carried out an extensive survey of radiative transitions in charmonium with detailed consideration of the phenomenology suggested by the results.
Notably we have performed the first lattice QCD calculation of the exotic η c1 radiative decay and found a large Γ(η c1 → J/ψγ) = 115(16) keV. We found statistically significant electric dipole and magnetic quadrupole form factors in χ c2 → J/ψγ, calculated for the first time in this framework, and have studied excited χ c1,2 transitions. Our results for the excited χ c2 states suggest that there could be a radially excited 2 3 P 2 state and a 1 3 F 2 state rather close in mass. This signal followed from a clear observation of dominance of E 1 over E 3 in one case and E 3 over E 1 in the other and matches the general expectations of potential models. This is relevant given the recent observation of a candidate 2 ++ excited state in γγ → DD [33].
Our calculation of magnetic dipole transition widths, such as ψ → η c γ and hindered excited state transitions reflected the expected suppression of the excited state transitions. We note that, modulo lattice systematic effects which can be reduced, our method is not troubled by the uncontrolled approximations and model dependence inherent in model calculations of such suppressed transitions. In light of this we were able to make statements regarding the possible influence of closed channel D-meson loops on the hindered transitions, determining that they cannot be as large as suggested in certain studies. We identified a putative nonexotic hybrid state having a large magnetic dipole transition amplitude to η c -the possibility that this reflects the non-spin-flip magnetic dipole excitation allowed within the flux-tube model (and likely within other models having more than cc content) was raised.
While the current numerical results may be affected by lattice systematic errors, future calculations using the now proven methodology can address these shortcomings by using dynamical lattices of sufficient size having a number of different lattice spacings. At some point after the introduction of dynamical lattices, the difficulty of dealing with resonant states embedded in a multi-meson continuum will have to be addressed.
We have demonstrated that our results are in agreement with general predictions of quark potential models but we are able to go beyond this to make statements about states in which there is an excited gluonic field. These results are appropriate for comparison with models proposing particular forms for the gluonic excitation.
The great advantage of this development within lattice QCD as compared to models relying upon the non-relativistic motion of quarks, is its immediate applicability to the light quark systems. Future efforts will consider photocouplings of light quark mesons, and in particular exotics, as these are of central important in the production rates for the GlueX project. If the large couplings we find for the 1 −+ state with heavy quarks persist into the light quark sector this will confirm the intuition and model results used to motivate the GlueX production process.

APPENDIX A: DERIVATIVE OPERATORS AT ZERO AND NON-ZERO MOMENTUM
In this appendix we give a few details of the spin-structure we deal with when considering two-and three-point correlators. We take the example of the T 1 −− "vector" operators considered at p = (000) and p = (100). Using the decompositions in the appendix of [1], we find that at zero momentum all the operators we have used have the following behaviour (at least in the continuum limit), and there are no non-zero overlaps onto particle of any other J P C (except those lattice artifacts that appear suppressed by powers of a). The three-point correlator then involves a sum λ ǫ * i ( 0, λ)ǫ j ( 0, λ) = −g ij . At non-zero momentum things are not so simple, with one unit of lattice momentum (e.g. p = (100)) we have the following non-zero overlaps in the continuum limit, These forms are derived following the decompositions given in Appendix A of [1] 27 . So here we have additional contributions from (0, 1, 2) +− particles -the (0, 2) +− contributions we can neglect as these states are very heavy, but the 1 +− we must worry about. We can see the entry of such particles in the two-point function spectral analysis; in figure  15 we show the spectrum at p = (100) with either i = x, where 1 +− does not contribute, or with i = y, z where 1 +− can contribute. Also shown are the T −− 1 and T +− 1 spectra extracted at p = (000) extrapolated to p = (100) using the continuum dispersion relation E = m 2 + | p| 2 . We see the entry of the lightest 1 +− state into the T −− 1 ( p = (100)) spectrum in the case i = y, z but not in the case i = x as expected. We can explicitly exclude this state from our three-point analysis by not projecting onto the eigenvector belonging to this level. In the results presented in figures 4,5 we separately consider the cases i = x and i = y, z shown by the triangles and the circles. Note that we did not simply rely upon energy matching to determine the levels -a more precise mapping between p = (000) levels and p = (100) comes from considering the overlaps, Z. As a concrete example consider the matrix elements 1 −− ( p, λ)|ψγ i ψ|0 = Zǫ * i ( p, λ) where Z is proportional to the vector decay constant and where, for unsmeared quark fields in the continuum limit, Z should be independent of momentum p. As can be clearly seen in Table IV, one can identify the J/ψ state with the lightest in each case, while the ψ ′ state is the next lightest in two of the cases, but is the third state in the case where the 1 +− can contribute. We considered the Z values for the entire set of operators when making the state assignments at finite momentum. This approach becomes increasingly more cumbersome as the momentum increases so we have not considered it any further. . State proposed to be J/ψ in bold, state proposed to be ψ ′ in italics.

APPENDIX B: IMPROVEMENT OF THE VECTOR CURRENT
In [24] the improvement of the vector current to O(ma) in a manner compatible with the improvement in the anisotropic Clover action was presented. Therein the analysis was for heavy-light currents and the renormalization constants were determined perturbatively. Here we are interested in heavy-heavy currents and we shall determine Z V non-perturbatively. Nevertheless we can consider modifying our local current to include extra terms suggested by the improvement scheme to investigate any change in the form-factor values. This can be considered to give a crude estimate of how much we might anticipate scaling to the continuum to affect our results.
Our Clover action uses r = 1 and as such our local vector current is not automatically improved -at tree level the improved current is given byΨγ µ Ψ where Ψ(x) ∝ 1 + a s d 1 γ j − → D j ψ(x) is the "rotated" field used in the construction of the improved action.
The improvement parameter d 1 = 1 4 (1 − ξr) + O(m 0 a t ) ≈ −0.5 for renormalised anisotropy ξ = 3.0. The quark mass parameter that appears in our Clover action, m 0 a t has the value 0.0401 which is clearly small, while the mass in spatial lattice units is three times as large, but still might be argued to be small.
We have attempted to construct the improved current by using the tree-level equation of motion a t γ 4 ∇ 4 + a s ν ξ 0 γ j ∇ j + m 0 a t ψ = 0 (B1) to eliminate the derivative in the expansion to O(a s ) of the improved currentΨγ µ Ψ. This yields the resultΨ where a common proportionality factor is ignored given that we determine Z V nonperturbatively using the pseudoscalar form-factor at Q 2 = 0. At q = (000) the second term does not contribute and we would expect ν d 1 ≈ 1.11 for the parameters in our action. Note that in section III A we found 1.11 (1) for this ratio using a non-perturbative extraction.
We computed form-factors using the improved currents for a few of the transitions considered in this paper. The only effect of a considerable size was found in the χ c0 → J/ψγ transition as shown in figure 16(a). Note that the addition of the improvement brings the Clover data into better agreement with the DWF data on the same lattices as we might expect given the automatic O(a) improvement one has with DWF. We note that for the same correlators projected on to the excited state ψ ′ such a large difference with respect to the local current was not seen, see figure 16(b). We were unable to find any other large effects due to improvement, e.g. consider J/ψ → η c γ shown in figure 17. These observations may suggest that except in certain notable cases (the scalar), the discretization errors on our results are relatively small. FIG. 17: J/ψ → η c γ transition form-factor using local and O(ma) improved current with anisotropic Clover quarks and the fit from a study using the local current with DWF fermions on the same anisotropic lattices.