Radiative Transitions in Charmonium from Lattice QCD

The coupling of various charmonium mesons to a photon is studied using lattice QCD, giving access to radiative form factors and transitions, and probing the mesons' structure. Methods are developed which allow the robust determination of amplitudes, including those involving an excited state as well as multiple form factors, for a range of kinematics. These are applied in dynamical lattice QCD calculations using the distillation technique to compute the underlying three-point correlation functions. Form factors and amplitudes involving the low-lying charmonia, $\eta_c$, $J/\psi$, $\chi_{c0}$ and $\eta_c'$, are calculated, demonstrating the methods and providing the first dynamical lattice QCD results for some quantities.


Introduction
The study of charmonium and charmonium-like mesons provides a useful window into the non-perturbative regime of QCD.The mass scale of the charm quark, falling between that of light quarks and the bottom quark, means that these systems are reasonably nonrelativistic and hence cleaner to study than light hadrons theoretically, but they are light enough to be accessible in many different experiments.A plethora of charmonia have been observed [1], including a number of states, known as exotic hadrons, that do not fit into the quark model picture of a meson as a bound state of a quark and an antiquark (see e.g.[2][3][4] for some reviews).
Spectroscopy, the study of the masses and widths of stable hadrons and hadronic resonances, is one tool that can be used to probe QCD.However, more detailed information on the structure of hadrons can be obtained from other quantities; for example, from radiative transitions and electromagnetic form factors, related to a hadron's coupling to a photon.
Radiative partial widths such as Γ(J/ψ → η c γ) and Γ(χ c0 → J/ψ γ) have been measured by, and are continuing to be studied by, a range of experiments.To complement experimental measurements, it is important to perform theoretical calculations to compare against and aid in the interpretation of experimental results, e.g. to determine whether a state is a conventional charmonium meson, a hybrid meson, a molecular state, a tetraquark, etc.There is a long history of computing radiative transition amplitudes involving charmonia in phenomenological models, but first-principles QCD-based calculations are also necessary.Lattice QCD, where the theory is discretised on a lattice in a finite spacetime volume and quantities are computed numerically, is a method that allows for systematically-improvable computations within QCD.Lattice calculations of a few radiative transition amplitudes involving low-lying charmonia, such as J/ψ → η c γ, ψ(2S) → η c γ, χ c0 → J/ψ γ, are now mature [5][6][7][8][9][10][11], but there have been fewer calculations involving higher-lying and excited charmonia [12][13][14].
The aim of this paper is to establish and test a robust methodology that can be applied to radiative transitions involving exotic and excited charmonia, building on work on transitions in the light-quark sector [15], studies of the charmonium spectrum [16,17], and earlier quenched calculations of radiative transitions involving charmonia [5,12].The work in the charmonium sector showed how lattice calculations and their phenomenological interpretation can be used to probe the structure of charmonia, e.g. to identify a hybrid meson candidate.We will discuss how amplitudes and form factors can be robustly extract from computations of three-point correlation functions involving carefully constructed interpolating operators.To demonstrate the techniques, transitions involving some low-lying charmonia, η c (J P C = 0 −+ ), J/ψ (1 −− ), χ c0 (0 ++ ) and η ′ c (0 −+ ), are studied in dynamical lattice QCD calculations.
The rest of this paper is organised as follows.In Section 2, radiative transitions, form factors and the Lorentz structure of matrix elements are discussed.Section 3 presents the techniques employed to extract the relevant information from lattice QCD calculations and Section 4 provides specific details of the calculations used in this study.Results from applying the methodology are presented in Section 5 and a summary is given in Section 6.

Radiative Transitions and Form Factors
The coupling of hadrons to photons induces radiative transitions between hadrons.At leading order in QED, the amplitude for h → h ′ γ is proportional to the transition matrix element involving the electromagnetic current, where the hadrons h, h ′ have angular momentum J, J ′ , helicity λ, λ ′ and 4-momentum p µ = (E h , ⃗ p), p ′µ = (E ′ h , ⃗ p ′ ).The electromagnetic current is j µ = i q i ψi γ µ ψ i with quark charge q i and the sum taken over quark flavours i.This matrix element can be decomposed as, where K µ i [h ′ J ′ (λ ′ , ⃗ p ′ ); h J (λ, ⃗ p )] are Lorentz vectors referred to as kinematic factors and F i (Q 2 ) are Lorentz scalars referred to as transition form factors (or just form factors if h ′ = h).The form factors only depend on , the virtuality of the photon.
The decomposition and kinematic factors for a particular matrix element can be constructed by considering all possible contractions of the polarization tensors, ϵ, and 4momenta associated with each hadron, provided that the contracted object carries a single Lorentz index, i.e. it has the correct Lorentz structure, and is linear in the polarization tensors.Constraints from Ward identities and discrete symmetries such as parity must then be taken into account, removing form factors that are redundant or not allowed.
As an example, the matrix element of the pseudoscalar meson, h ′ = h = η c , has a decomposition involving a single form factor, (2.3) For a vector meson, h ′ = h = J/ψ, there are three independent form factors There is no unique decomposition of a matrix element in terms of form factors: as well as the freedom to divide out a form factor by a Lorentz scalar and put it in the corresponding kinematic factor, when more one form factor is present any linearly independent set of form factors can be used.In this paper we generally use the multipole basis [18].
From the form factors we can infer qualitative and quantitative information on the hadrons.For example, by looking at the hierarchy of form factors partial radiative widths and charge distributions which are useful for comparison to experiment can be determined.The partial radiative width for h → h ′ γ can be related to the multipole form factors by [12], where ⃗ q = ⃗ p ′ − ⃗ p is the spatial component of the photon momentum, m and J are the mass and angular momentum of hadron h respectively, and the sum is over all the relevant transverse multipole form factors, {F k (Q 2 )}.In the next section we describe how form factors are determined in lattice QCD calculations.

Methodology
Lattice QCD calculations are performed in a finite volume and so the spectrum is discrete.If h and h ′ are stable with respect to QCD (i.e. they do not decay strongly), and n i and n f are the corresponding finite-volume energy eigenstates then matrix elements describing a transition h → h ′ γ are determined by considering three-point correlators of the form, where Ω † n i (0) and Ω n f (∆t) are operators located on timeslices 0 and ∆t which create n i and annihilate n f respectively.The reduced symmetry of the finite-volume lattice compared to an infinite-volume continuum introduces some complications.Spin, J, is no longer a good quantum number and instead states are classified by the irreducible representations (irreps) of the relevant finite symmetry group.All the operators used in this work are projected onto definite momentum and constructed to transform in the appropriate irrep as described in [19,20]; more details of the consequences of the reduced symmetry for extracting matrix elements can be found in Section D of [15].Before discussing how to extract matrix elements from the correlators, we will first describe the construction of operators that efficiently interpolate the states of interest and this requires an analysis of two-point correlators.

Lattice Spectroscopy and Optimized Operators
We form an 'optimized operator' Ω † n as a linear combination of a large basis of operators , where coefficients {v i } are found such that there is a large overlap ⟨n| Ω † n |0⟩ onto the desired state n, but very small overlap with other states.Following the method described in [19,20], in each channel a large basis of fermion-bilinear operators is used, built from Dirac gamma matrices and gauge-covariant derivatives, projected onto the desired quantum numbers and momentum and transforming in the appropriate finitevolume irrep.With N operators an N ×N matrix of two-point correlators can be computed, Optimal coefficients are then determined following the methods of [21,22]: a generalized eigenvalue problem for the correlator matrix is solved, where t 0 is an appropriate reference timeslice.The generalized eigenvalues, λ n (t, t 0 ), are related to the energies of the states as λ n (t) ∼ e −En(t−t 0 ) , where E n is the energy of the n'th state in this channel.The generalized eigenvectors v give the variationally-optimal coefficients and optimised operators are then constructed as, The normalisation of the operator, i.e. the overlap ⟨n| Ω † n |0⟩, depends on the details of the implementation of the variational method and the choice of t 0 -this is discussed in Appendix A.
As an example to demonstrate the effectiveness of these optimised operators, Figure 1 shows two-point correlators in the η c channel at zero momentum with the leading time Illustration of the effectiveness of optimised operators in the η c channel at zero momentum.Points show correlators with the leading time dependence divided out as described in the text.Shown in blue and orange are correlators with optimised operators Ω 0 and Ω 1 designed to overlap onto the η c and η ′ c respectively at both source and sink.In red is the correlator obtained with a simple fermion bilinear operator O 0 ∼ ψγ 5 ψ at both source and sink, normalised so that the plateau value is 1 at large time.In green is the "off-diagonal" correlator with Ω 1 at the source and Ω 0 at the sink -these points have been shifted up by 1 to display them on the same plot.dependence divided out, i.e.C(t) e En(t−t 0 ) , computed using the lattice setup described below in Section 4. The correlator with an optimised ground-state η c operator Ω 0 at the source and sink (blue points) plateaus at a much earlier time than the correlator from a simple fermion-bilinear operator O 0 ∼ ψγ 5 ψ (red points) with the plateau starting at t/a t ≈ 10 compared to t/a t ≈ 18.This means that signals can be extracted at earlier times where the correlator is more statistically precise and for three-point correlators this enables the matrix element to be determined more precisely.In addition, a good signal is found with the optimised operator Ω 1 for the first excited state η ′ c (orange points), something which would not be possible with a simple fermion-bilinear operator.This plateau value can be slightly less than 1 because of the normalisation of the operators -for more details see the discussion around Eq. (A.6) in Appendix A. The green points show the correlator with Ω † 1 at the source and Ω 0 at the sink -as well as having the leading time dependence divided out, these have been shifted up by 1 for display on the same plot.This early plateau to 1, i.e. plateau of unshifted correlator to 0, demonstrates the good orthogonality between these two operators.

Correlator Analysis and Extracting Form Factors
Having established how to construct operators that efficiently interpolate states of interest, we now turn focus back to three-point correlators of the form in eq.(3.1).Exposing the overlap of the operator onto the desired state and a single additional state, . ., where ε n ′ i is small and the ellipsis denotes other suppressed terms from overlaps onto higher states, we can expand the correlator as The desired matrix element can then be extracted from . ., where terms involving excited state matrix elements are combined in the f n ′ f , f n ′ i terms.It can be seen that at appropriate intermediate times 0 < t < ∆t, the second and third terms should be highly suppressed and the correlator should plateau to a value corresponding to the matrix element.
The general kinematic decomposition outlined in eq.(2.2) shows that multiple form factors can contribute to a single matrix element.If only a single kinematic factor appears, then dividing the matrix element by the kinematic factor, K µ , for K µ ̸ = 0, yields the desired form factor, where there is no summation over µ.In this case, we consider a time-dependent form factor, F (Q 2 ; t) = Γ(Q 2 ; t)/K µ , and fits of the time dependence to a constant, F (Q 2 ; t) = F (Q 2 ), i.e. assuming contamination is negligible, and to The final fit form and time range used are chosen to achieve a good fit while permitting as large a time range as possible.When defining the χ 2 to be minimised in the fitting procedure, the covariance matrices used to account for correlations between different timeslices are estimated using a "shrinkage" approach presented in [23][24][25][26].In particular, we use an implementation [27] of the "Oracle Approximating Shrinkage Estimator" for the covariance matrix.
0, t) obtained from a three-point correlator featuring an optimised η c operator with ⃗ p = 2π L (0, 0, 1) at both the source and the sink.The black line shows the result of a fit to a constant for 10 ≤ t/a t ≤ 30 with the band indicating the ±1σ statistical uncertainty (darker points are included in fit whilst translucent points are not).
As an example, Figure 2 shows F (Q 2 = 0, t) obtained from a three-point correlator featuring an optimised η c operator with ⃗ p = 2π L (0, 0, 1) at both the source and the sink, computed using the lattice setup described below in Section 4, where L is the spatial extent of the lattice.This can be seen to plateau rapidly at intermediate times, 0 < t < ∆t, allowing a fit over a large time range and so enabling the form factor to be extracted robustly and with high statistical precision.
To allow propagation of correlations into subsequent analyses, such as the Q 2 dependence of the form factor, the fitting procedure is repeated on jackknife samples of the data.This allows an additional verification of the robustness of the fit under small perturbations.Fits to subsets of data consisting of only correlators containing a temporal current or only those containing a spatial current are also considered.
When multiple form factors enter, there is in general an under-specified linear system at each Q 2 because a single matrix element cannot determine multiple form factors.In this case we bin together sets of matrix elements with the same or similar values of Q 2 .There is a compromise between the number of data points in a given bin and the range of Q 2 that a bin spans.Bigger bins provide more constraints on the form factors in a bin, but at the expense of assuming the form factors do not vary appreciably over a larger Q 2 interval (i.e. over the bin) and having fewer form factor estimates across the range of Q 2 .
From the relativistic dispersion relation, it can be seen that if two matrix elements have the same initial momentum and the same final momentum they will have the same Q 2 .Such a redundancy can be generated by varying the Lorentz index of the current and the spin/helicity component of any hadrons with non-zero spin.A more coarse-grained approach is to group matrix elements that have the same magnitude of current momentum , ⃗ p ′ 2 then these matrix elements will have similar Q 2 .Enumerating the matrix elements within a Q 2 bin with an index a, the linear system can be written as, where i enumerates the form factors in the decomposition.Expressing the matrix of kinematic factors as K and the vectors of form factors and matrix elements as F and Γ, these form factors can be determined by linear regression [15] using the normal equation, We then fit the time dependence of the form factor to obtain F i (Q 2 ) as in the case described above where only a single form factor appears.To check the robustness of the form factor determinations, we repeat the process a number of times systematically removing different matrix elements in the bin from consideration.Equation (3.5) is solved on each time slice and each gauge configuration (via jackknife) to yield N cfgs samples of F i (Q 2 ; t).This binning procedure can also be used when only a single form factor appears, but in this case it is generally preferable to consider each matrix element separately because this allows more straightforward independent cross-checks, and the isolation of contributions from different lattice irreps and from the spatial current v.s. the temporal current.
With multiple form factors an alternative approach would be to fit the time dependence of Γ a (Q 2 ; t) to get Γ a (Q 2 ), and from this obtain F i (Q 2 ).A potential advantage of this method is that the expected time dependence of contamination from unwanted states in Γ a (Q 2 ; t) is clearer than in F i (Q 2 ; t).However, because we use optimised operators, we find that matrix elements plateau a relatively short distance away from the source/sink (see Fig. 1) and a reliable plateau in F i (Q 2 ; t) is consistently found.We find that our procedure leads to more stable determinations of the form factors -the alternative method requires very many fits to the time dependence and when a form factor is associated with a small kinematic factor a small change to one fit can have a significant impact on the result.Where both methods give robust determinations of the form factors, the results from each are consistent.

Exp1
F 0 e Table 1.Forms considered when fitting the Q 2 dependence of form factors in this study.
Equipped with a set of form factors {F (Q 2 )}, the Q 2 dependence can be parameterised to enable extrapolation/interpolation to Q 2 = 0. To avoid model dependence from assuming a particular Q 2 form, we attempt fits with the variety of phenomenologically-motivated fit forms listed in Table 1.As discussed later, we discard any fits that are significantly worse at describing the data.Correlations are again taken into account through a covariance matrix that is estimated using the shrinkage approach as described above.

Calculation Details
This study uses a single ensemble of gauge fields, generated with a tree-level Symanzikimproved gauge action and a tadpole-improved Sheikholeslami-Wohlert (clover) fermion action with stout smeared spatial gauge fields and 2+1 flavours of dynamical quarks (two degenerate light up and down quarks, and a heavier strange quark) [28,29].Correlation functions were computed on 603 gauge configurations.The discretisation is anisotropic with a spatial lattice spacing, a s ∼ 0.12 fm, and a finer temporal lattice spacing, a t , with anisotropy ξ = a s /a t ∼ 3.5.The lattice volume is (L/a s ) 3 × (T /a t ) = 20 3 × 128 -this relatively large volume with m π L ∼ 5 means that finite-volume effects should not be significant for charmonia below D D threshold as can be seen in [16].The cubic spatial volume with periodic boundary conditions quantizes momentum as ⃗ p = 2π L (n x , n y , n z ), n i ∈ Z, and in this work we will consider momenta with n 2 = n 2 x + n 2 y + n 2 z ≤ 4. Strange quarks are tuned to approximate their physical mass while the unphysicallyheavy light quarks yield m π ≈ 391 MeV. 2 The tuning of the parameters in the charm-quark action and the impact of any imperfect tuning are described in Ref. [16].To present results in physical units, the lattice spacing is determined using the mass of the Ω baryon which (L/a s )3 × (T /a t ) m π / MeV N cfgs N vecs a s / fm m π L 20 3   Ω atm lat.Ω = 5.666 GeV [30].A summary of the lattice parameters can be found in Table 2.
We compute correlators using the distillation framework [31] and further details on the calculation of three-point correlators in this approach can be found in [15].Quark fields at the source and the sink are smeared using the distillation operator which projects onto the low energy modes.The quark fields in the vector current insertion are unsmeared because otherwise a complicated energy-dependent renormalisation would be required to connect with physical amplitudes.The operators at the source and sink, and the current are each projected onto definite momentum, and the source-sink separation is ∆t/a t = 40.
Diagrams representing the relevant Wick contractions for the coupling of a photon to a meson are shown in Fig 3.In this work we only compute the contribution where the photon couples to the charm quark (a) and not where the photon couples to the charm antiquark (b), and we do not explicitly include the electric charge of the quark -this means that we actually compute scaled versions of the multipole form factors, Fk , defined by . It also implies that charge-conjugation (C) violating processes can be investigated, including the form factors of charmonia, and used to probe the structure of states. 3In addition, disconnected diagrams like (c) where the charm quark and antiquark annihilate are not computed -these are OZI suppressed and are expected to only give small contributions in charmonium [32].We also neglect diagrams where the photon couples to a disconnected quark loop (d) which are again expected to only give small contributions -the contribution from a charm quark loop is suppressed by the large charm-quark mass and in the SU(3) flavour symmetry limit the contributions from light and strange quark loops cancel.
Working with a finite lattice spacing leads to discretisation effects and with only one lattice spacing we can not fully quantify them.The hyperfine splitting, ∆m = m J/ψ − m ηc , is particularly sensitive to these and on this lattice we find ∆m lat.= 80.36 (16)MeV compared with the experimental value ∆m exp.= 113.0(4)MeV-Ref.[16] discusses this tension and a rough estimate of the size of the discretisation uncertainties.
We use an anisotropic lattice with different temporal and spatial lattice spacings and, while parameters in the action are tuned to restore hypercubic symmetry as far as possible, this leads to an additional type of discretisation effect.The anisotropy experienced by a stable hadron can be determined from the relativistic dispersion relation, where m is the mass of the hadron and E is its energy at momentum magnitude The mass and anisotropy determined from fits to eq. ( 4.1) with n 2 ≤ 3 for relevant charmonia are presented in Table 3 and, as an example, the momentum dependence of the η c energy is shown in Figure 4.The fits are generally found to describe the data reasonably well without the need to add higher-order discretization terms such as (a s p) 4 .
The measured anisotropies are all similar and in line with what has been found in other flavour sectors on this lattice [16,[33][34][35][36], suggesting that these discretisation effects are small.The small discrepancy between the values of ξ measured for the different helicity magnitudes of the J/ψ (λ = 0 and |λ| = 1) is due to the reduced symmetry of the lattice and is consistent with what has been observed in other studies [20,36,37].To estimate the systematic effect from the small discrepancies between different determinations of the anisotropy, tests were performed with ξ varied over an envelope of ξ n i ± σ, ξ n f ± σ, where σ is the corresponding statistical uncertainty.In all tests, the resulting effect is much smaller than the statistical uncertainty associated to the form factors. 4 In addition, we avoid using higher momenta where discretisation effects are expected to become more important.In the next section we discuss the renormalisation of the vector current and how we use an improved current to reduce discretisation uncertainties in calculations of the three-point correlation functions.

Renormalizing and Improving the Vector Current
The local vector current is not conserved in our lattice formulation and must be renormalised in order to relate the extracted matrix elements to physical results.We perform a non-perturbative renormalization by imposing that F ηc (Q 2 = 0) be unity.The matrix elements we extract must therefore be multiplied by, On an anisotropic lattice, the renormalisation constants for the temporal and spatial currents, Z t V and Z s V , are distinct and must be determined seperately.Following the methods of Section 3.2, F lat. ηc (0) is extracted from ⟨0| Ω ηc (∆t, ⃗ p) j µ (t, ⃗ q = 0) Ω † ηc (0, ⃗ p) |0⟩ for a range of momenta ⃗ p and the resulting Z V are shown in Figure 5. Fitting the data for to a constant yields, where χ 2 /N dof = 14/(4 − 1) = 4.7 and χ 2 /N dof = 24/(7 − 1) = 4.0 for the temporal and spatial fits respectively, and the uncertainties are statistical.Results presented in subsequent sections have all been appropriately renormalized.
To reduce discretisation effects we add an O(a) improvement term to the local vector current which, for the fermionic action used here, leads to improved renormalised temporal and spatial vector currents [15], where ν s = 1.078 is a parameter in the fermion action [28] related to the spatial clover coefficient.Note that the improvement term vanishes for Q 2 = 0 and so it does not modify the determination of Z V .The effect of the improvement term was investigated and it was found to improve the agreement between the determination of physical observables from the spatial and temporal currents.As a demonstration, the Q 2 dependence of the η c radiative form factor is shown using both the unimproved and the improved vector currents in Figure 6.The most pronounced effect of the improvement is to reduce the discrepancy between results obtained from the temporal and spatial current insertions.In addition, the spread for a given temporal/spatial current with different source and sink momentum combinations is reduced.This effect appears to be more pronounced at larger Q 2 where the coefficient of the improvement term and discretization effects increase.All the results we present in subsequent sections were obtained using the improved currents.

Improved
Figure 6.The η c radiative form factor, F ηc (Q 2 ), defined in eq.(2.3), plotted against Q 2 calculated using the unimproved current (left plot) and the improved current (right plot).Each point corresponds to a single correlator involving a temporal (red) or spatial (blue) current with the methods of Section 3.2 used to extract the form factor and the appropriate renormalisation factor (Z t V or Z s V ) applied.The effectiveness of the improvement term can clearly be seen in the reduced spread of the form factors for each Q 2 in the right plot compared to the left plot.

Results
In this section results for radiative and transition form factors involving some S and Pwave charmonia are presented.In all cases the improved vector current was used and the renormalisation factor has been applied as discussed in Section 4.1.Comparisons to previous lattice calculations and experimental results are discussed where appropriate.Our aim is to develop and test techniques that can be used to study transitions between various charmonia and to gain insight into their structure, and not to quantify all the systematic uncertainties.The uncertainties quoted in final results are obtained by taking a range over the allowed fit forms, variations over the masses and anisotropies, and the statistical uncertainties.

η c Form Factor
The simplest application of the methodology outlined above is to cases with a single form factor, and perhaps the most straightforward to consider is the form factor of the lowestlying charmonium meson, the η c .The relevant kinematic decomposition is given in eq. ( 2.3).This transition can not occur in nature since it violates charge-conjugation parity conservation, or, looked at another way, coupling of the current to the charm quark would cancel with the contribution from coupling to the anti-charm quark.Nevertheless, it can be determined in a lattice calculation by only coupling the current to the charm quark as discussed in Section 4, serves as a useful test of lattice methods and can provide information about the meson.
To determine the form factor, three-point correlators were computed for a range of Q 2 values by varying the source and sink momenta, ⃗ p and ⃗ p ′ .Following the approach outlined above, each correlator was divided by the leading exponential time dependence and the associated kinematic factor, and the residual time-dependence was fit to yield F ηc (Q 2 ) for each correlator.The Q 2 dependence of the form factor was fit using the parameterisations in Table 1 -summaries of the results are given in Table 4 in Appendix B. The extracted F ηc (Q 2 ) are shown as points in Figure 7 and the bands indicate the results of the fitsthose which have a considerably worse goodness of fit, shown in italics in the table, are excluded from the plot and subsequent analysis.In all the parameterisations, F ηc (0) is a free parameter and its fitted value is unity within statistical uncertainties as expected from the renormalisation of the vector current, giving confidence in the results and the determination of Z V .
While the fits appear to reasonably capture the Q 2 dependence of the form factor, even the best have χ 2 /N dof ∼ 4. We note that the data points are statistically precise but with an unquantified systematic uncertainty and this combination may explain the poor goodness of fit.The spread of data points at each Q 2 (resulting from different combinations of source and sink momenta) is an indication of lattice artifacts.When values of the form factor determined at larger Q 2 , corresponding to correlators with |⃗ p| ≥ 2 2π L at the source, sink and/or insertion, are included in fits, the χ 2 /N dof deteriorates further lending support to the argument that lattice artifacts are a significant factor causing this poor goodness of fit.In agreement with Ref. [5], for parameterisations with two parameters we find that the L are translucent and not included in the fits.Bands show results of fits to the Q 2 dependence, as described in the text, giving 1σ statistical uncertainties around the mean.The dashed grey line indicates Q 2 = 0.
VMD-inspired fit form does a poor job of describing the data in contrast to the Gaussian form.
The charge radius, defined by can be determined from the Q 2 dependence of the form factor at Q 2 = 0 and compared to previous lattice results and other approaches.To investigate if there is any systematic difference between the charge radius extracted from a spatial or temporal vector current, as well as fits to all the data (Table 4), we also consider fits involving only data from correlators featuring a temporal insertion and only data from those featuring a spatial insertion -summaries of the results are given in, respectively, Tables 5 and 6 in Appendix B. It can be seen that the spread of results from fitting to these subsets of the data is similar in size to the spread from different fit forms.Taking an envelope over the fits, excluding those shown in italics in the tables, yields, where the quoted value and uncertainty have been chosen to encompass the spread over statistical uncertainties from the different fits -the spread is much larger than the statistical uncertainty from any given fit and the uncertainty from varying m ηc and ξ ηc within their uncertainties.
These results include quenched calculations [5] as well as those with dynamical light quarks (N f = 2) [6, 10], compared to this work which has dynamical light and strange quarks (N f = 2 + 1), used different lattice actions and methods, and do not quantify systematic uncertainties, and so a detailed comparison is not appropriate.However, we note that our result is within the range of those values.Now that we have demonstrated the approach in a calculation of the η c form factor, we will test it further by studying the form factor of an excited pseudoscalar, the η ′ c .

η ′ c
Form Factor A natural extension of the calculation of the η c form factor is to consider the form factor of the first excited pseudoscalar charmonium meson, the η ′ c .Correlators involving the η ′ c are expected to be statistically noisier than those involving the η c because it is an excited state.In addition, the reduced symmetry at non-zero momentum means that the η ′ c can be up to the third excited state in a given lattice irrep 6 .Therefore, a reliable determination of the form factor relies on the use of appropriate operators in the three-point correlators.These additional challenges mean that it provides a good test of the methodology.As seen in the example case in Figure 1, the optimised operators constructed following the methods of Section 3.1 perform well in both extracting a signal for the η ′ c and minimising contamination from the subleading η c ground state.
The kinematic decomposition for the η ′ c is identical to that for the η c in eq.(2.3).The analysis proceeds as in Section 5.1 and the resulting form factor is shown in Figure 8 along with results of fits to the Q 2 dependence.Summaries of the fits are given in Table 7 in Appendix B. Additional fits (not shown on the plot) to only data featuring a temporal current or a spatial current are shown in Tables 8 and 9 respectively.As expected, the data are considerably noisier than for the η c and, probably due to this, the goodness of fits are much better, χ 2 /N dof ∼ 1.5.In all cases F (0) is a free parameter and its fitted value is consistent with unity.The determination of Z t V and Z s V did not use the η ′ c and so, in contrast to the η c form factor, this provides a non-trivial check.
We determine the η ′ c charge radius in the same way as for the η c .Taking an envelope over fits yields, (18) fm. 6where the χc1, χc2 and ηc can all appear with lower energy in the same channel This is larger than the η c charge radius, as one would expect for a radial excitation.It also has a much larger statistical uncertainty, as expected for an excited state.In the next section, instead of a radially excited charmonium meson we study an orbital excitation, the scalar χ c0 .0.0 0.2 0.4 0.6 0.8 1.0 7 but for the η ′ c form factor.

χ c0 Form Factor
As for the η c , the radiative form factor of the χ c0 is determined by computing correlators with a χ c0 at both the source and sink.The kinematic decomposition is the same as for the η c in equation (2.3) and the analysis proceeds as before.The form factor and results of fits to the Q 2 dependence are shown in Figure 9. Summaries of the fits are given in Table 10 in Appendix B -additional fits (not shown on the plot) to only data featuring a temporal current or a spatial current are shown in Tables 11 and 12 respectively.The data here are noisier than the ground state η c but less so than for the η ′ c .In all fits F (0) is a free parameter and is found to be one within statistical uncertainty, giving further confidence in the results and the vector current renormalisation.
As previously discussed, results in the literature used different fermion discretisations, different numbers of dynamical flavours and have unquantified systematic uncertainties and so should not be compared in detail.Nevertheless, we find agreement with [5] and the very broad estimate of [10], though there is some tension with the result in [6].
The hierarchy of charge radii, r 2 ηc , is consistent with expectations from quark-potential models where the η c , χ c0 , η ′ c are nL = 1S, 1P , 2S states respectively, with n the principal quantum number and L the orbital angular momentum [5].Now that we have demonstrated the methods in some simple cases, we move to experimentally-accessible transitions from which partial widths can be determined and then more complicated cases involving more than one form factor.

J/ψ → η c Transition Form Factor
The simplest physical process to study is the J/ψ → η c γ radiative transition for which the kinematic decomposition is There is still only a single form factor present but the Lorentz structure is different from that for the η c and involves a polarization vector, ϵ σ (λ, ⃗ p ), where λ is the helicity of the J/ψ.
Following the same analysis procedure as above, the resulting form factor is shown in Figure 10 along with the results of fits to the Q 2 dependence.Summaries of the fits are given in Table 13 in Appendix B. Additional fits (not shown on the plot) to only data featuring a temporal current or a spatial current are shown in Tables 14 and 15 respectively.
The points with small Q 2 , which correspond to the source and sink having identical momenta, and one point with Q 2 ∼ 1 (GeV) 2 have much larger error bars than the other points.For these points the kinematic factor is proportional to a t ∆E = a t (E ηc − E J/ψ ) which is small (≪ a t E ηc ).These kinematic factors are much smaller than the others and so the signal in the relevant correlators is also much smaller leading to a larger statistical uncertainty on the extracted form factor.Excluding these points from the fits does not change the results significantly because of their large uncertainties.
The form factor at Q 2 = 0 can be related to the radiative partial width by where ⃗ q is the momentum of the photon in the rest frame of the J/ψ, which depends strongly on the hyperfine splitting ∆m = m J/ψ −m ηc .There is a discrepancy between the hyperfine splitting measured on this lattice [16] and the experimental value, likely due to it being very sensitive to lattice artifacts.Therefore, the value of |⃗ q | is dependent on whether lattice or experimental masses are used.This issue also affects the kinematic factors at small Q 2 with large uncertainty discussed above. 7 One possible approach is to compare |F J/ψ ηc (0)| rather than Γ(J/ψ → η c γ), avoiding the ambiguity in the choice of ∆m.In addition, any points with a kinematic factor proportional to the hyperfine splitting can be excluded from the fits, though we find that excluding such points does not change the results significantly because they have large uncertainties.Following this approach and taking a spread over fit forms we find, 7 There is a negligible difference between ∆m and ∆E when the current has zero spatial momentum.
This range of |F (0)| lies at the upper end of the other lattice estimates [5-8, 12, 38] but within 1σ of most results presented.Lattice results are consistently larger than the derived PDG value of F J/ψ ηc (0) = 1.57(18) [1].However the most recent CLEO and KEDR results [39,40] have central values |F (0)| = 1.62(15), 2.06 (11), in better agreement with this study. 8ll the matrix elements considered thus far have only a single form factor in their decompositions.In general this is not the case and to test the methods further we will now consider two channels where multiple form factors occur.factors of the J/ψ -these are unphysical because they violate charge-conjugation parity conservation, but they provide information on the J/ψ and serve as a further test of the methodology.-25 -

J/ψ Form Factors
A vector meson has three radiative form factors and the relevant matrix element can be written as in eq.(2.4).The form factors G 1 , G 2 and G 3 therein are related to the charge E 0 , magnetic-dipole M 1 and electric-quadrupole E 2 form factors by, (5.8) ) and E 2 (Q 2 ), determined using the methods outlined above, are shown in Figures 13,14  For the charge form factor, the Q 2 dependence clearly follows a pattern similar to that seen in the η c , η ′ c and χ c0 charge form factors (Figures 7, 8 and 9).We find E 0 (0) = 1.005(4), roughly consistent with unity as expected from the vector current renormalisation.Following the same procedure as for the η c and taking an envelope over the fits gives a charge radius, which is comparable to r 2 ηc 1 2 = 0.243 (7) fm, as would be expected in a quark model as the η c and J/Ψ are both 1S states and so have similar spatial wavefunctions.This result is within the range of previous values in the literature, r 2 J/ψ 1 2 = 0.257(4) fm and 0.2570(38) fm, from [5] (quenched) and [10] (N f = 2) respectively.
In contrast to E 0 (Q 2 ) and M 1 (Q 2 ), the E 2 (Q 2 ) data is much noisier and the Q 2 dependence is not clear.The kinematic factor associated to the quadrupole form factor, K E 2 , is much smaller than those for the other form factors leading to the magnitude of the contribution to the overall matrix element, K E 2 E 2 (Q 2 ), being much smaller than the statistical error on the matrix element.This makes extraction of E 2 (Q 2 ) challenging.The kinematic factor is larger at larger Q 2 improving the quality of the determination -this can also be seen in [5] where the points at larger Q 2 have smaller statistical uncertainties.Given these large statistical uncertainties, we note that fitting to functions with many parameters is not appropriate.As such, in Table 19 we only consider fits to a constant, E 2 (Q 2 ) = E 2 (0), and two-parameter fit forms.Taking a spread over them yields E 2 (0) = −0.172(15).This is comparable with previous lattice determinations which found E 2 (0) = −0.109(56)[9] and −0.23(2) [5].Despite the challenges involved, it can be seen that the methods developed in this work can be used to robustly and precisely determine multiple form factors.

Conclusions
In this study, building on previous work by the Hadron Spectrum Collaboration [15], a lattice QCD methodology to compute radiative form factors and transition amplitudes involving charmonia has been developed.The methods have been tested in computations of amplitudes involving some of the lower-lying charmonia, η c , J/ψ, χ c0 and η ′ c .While we should reiterate that we have used a single lattice spacing and not quantified the various systematic uncertainties, the results demonstrate that a high level of statistical precision can be achieved, even for excited states, and multiple form factors can be robustly determined.Whilst the results at Q 2 = 0 have been highlighted, transitions and form factors have been determined at a range of momenta combinations offering insight into the Q 2 dependence.These are the first studies of a number of these amplitudes in a lattice QCD calculation with dynamical light and strange quarks.Now that the techniques have been demonstrated, in future work they will be applied to study radiative transitions involving higher-lying exotic charmonia.As has been shown in quenched lattice QCD calculations [12], studying radiative transitions can provide information on the structure of charmonia and aid in their identification.For example, it will be interesting to study transitions involving hybrid charmonia with J P C = 1 −− and J P C = 1 −+ -there are suggestions of experimental candidates for the former, but not currently for the latter [1,2].To study charmonia unstable with respect to QCD, the methods can be combined with those developed and demonstrated in Refs.[41][42][43][44].These calculations are timely given ongoing investigations at BESIII and the LHC, and the planned PANDA experiment.
Advanced Scientific Computing Research and Office of Nuclear Physics, Scientific Discovery through Advanced Computing (SciDAC) program; also acknowledged is support from the Exascale Computing Project (17-SC-20-SC), a collaborative effort of the U.S. Department of Energy Office of Science and the National Nuclear Security Administration.
The software codes Chroma [48] and QUDA [49,50] were used to compute the propagators required for this project on clusters at Jefferson Laboratory under the USQCD Initiative and the LQCD ARRA project.Gauge configurations were generated using an award of computer time provided by the U.S.

Data Access Statement
Results of the fits to the Q 2 dependence of the form factors are given in the supplemental material.Reasonable requests for further data such as matrix elements or correlation functions can be directed to the authors and will be considered in accordance with the Hadron Spectrum Collaboration's policy on sharing data.for large t, and so, comparing eq.(A.Thus, for each three-point correlator there are two independent normalization factors to be corrected for, one depending on n and one on m, when determining the matrix elements ⟨n| j µ (0) |m⟩.So for all three-point correlators, the matrix elements are given by, ⟨n| j µ (0) |m⟩ = e En(∆t−t) e Emt (1 where the quantites E n , E m and A n , A m are determined on each configuration via jackknife.

B Form Factor Parameterizations
In this appendix we present information about the fits to the Q 2 dependence of the form factors using the parameterisations in Table 1.Further details of the fits can be found in the supplemental material.Tables 4, 5 and 6 show the fits for the η c form factor including all data, only data extracted from correlators featuring a temporal current insertion and only data extracted from correlators featuring a spatial current insertion respectively.Tables 7, 8 and 9 give analogous results for the η ′ c form factor, Tables 10, 11 and 12 results for the χ c0 form factor, and Tables 13, 14 and 15 results for the J/ψ → η c γ form factor. Results of fits for the χ c0 → J/ψ γ E 1 form factor are shown in  4 but for fits to the η c form factor involving only data extracted from correlators featuring a spatial current insertion.The fit shown in italics has a considerably worse goodness of fit compared to the others and is excluded from subsequent analysis.

Figure 3 .
Figure 3. Relevant Wick contraction topologies.Diagrams (a) and (b) show fully connected diagrams, with the photon coupling to the charm quark and charm antiquark respectively.Diagrams (c) and (d) are disconnected.There is also another disconnected diagram like (c) where the photon couples to the sink rather than the source and another like (d) where the charm quark and antiquark annihilate as in (c).As discussed in the text, these disconnected diagram are expected to give only very small contributions and in this study only diagram (a) is computed.

Figure 4 .
Figure 4. Momentum dependence of the η c energy on the lattice used in this study.The points show the computed energies and the line shows the result of a fit to the relativistic dispersion relation, eq.(4.1), which has χ 2 /N dof = 1.64/(4 − 2) = 0.82.The point in grey with n 2 = 4 was not included in the fit.The inset shows the residuals, a t (E − E fit ) × 10 4 .

Figure 5 .
Figure 5. Renormalization factors as defined in eq.(4.2) for the temporal current Z t V (red points) and the spatial current Z s V (blue points), for a range of equal source and sink momenta ⃗ p = 2π L ⃗ n, determined using the methods discussed in Section 3.2.Where multiple points occur at a given n 2 they are slightly horizontally displaced to aid visibility.The red and blue bands are from fits to a constant and show the mean and ±1σ statistical uncertainty.Points corresponding to correlators where the source and sink have |⃗ p | 2 ≥ 4 2π L 2 are translucent and are not used in the fits.

2 )Figure 7 .
Figure 7.The η c form factor, F ηc (Q 2 ).Red and blue points are the form factor extracted from correlation functions featuring a temporal or spatial current insertion respectively.Points corresponding to correlators where at least one of the source, sink and current have momentum magnitude ≥ 2 2πL are translucent and not included in the fits.Bands show results of fits to the Q 2 dependence, as described in the text, giving 1σ statistical uncertainties around the mean.The dashed grey line indicates Q 2 = 0.
and 15 respectively.Results of fits to the Q 2 dependences are also shown in the figures and summaries of the results are given inTables 17,18 and 19 in Appendix B. Due to the structure of the kinematic factors, direct calculation at Q 2 ≈ 0 is only possible for E 0 , while the other form factors rely on an extrapolation in Q 2 .

Table 2 .
Lattice parameters for the lattice used in this study.N vecs is the number of distillation vectors.

Table 3 .
Relevant hadron masses and anisotropies from fits to eq. (4.1) for n 2 ≤ 3.5 Department of Energy INCITE program and supported in part under an ALCC award, and resources at: the Oak Ridge Leadership Computing Facility, which is a DOE Office of Science User Facility supported under Contract DE-AC05-00OR22725; the National Energy Research Scientific Computing Center (NERSC), a U.S. Department of Energy Office of Science User Facility located at Lawrence Berkeley National Laboratory, operated under Contract No. DE-AC02-05CH11231; the Texas Advanced Computing Center (TACC) at The University of Texas at Austin; the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation Grant No. ACI-1548562; and Jefferson Lab.

Table 16 ,
and results for the J/ψ E 0 , M 1 and E 2 form factors are shown in Tables 17, 18 and 19 respectively.Any fits that are rejected are shown in italics and are not included in subsequent analyses.The reason(s) for rejecting the fits are given in the table caption.

Table 4 .
Fits to the Q 2 dependence of the η c form factor with the fits including data extracted from correlators featuring a temporal or spatial current insertion.For each fit, the form factor at Q 2 = 0, F 0 , charge radius, r 2 ηc 1 2 , number of parameters, N p , and goodness of fit are given, where N dof is the number of degrees of freedom.

Table 5 .
As Table4but for fits to the η c form factor involving only data extracted from correlators featuring a temporal current insertion.The fit shown in italics has a considerably worse goodness of fit compared to the others and is excluded from subsequent analysis.

Table 6 .
As Table