A Rotations and helicity operators

# Excited meson radiative transitions from lattice QCD using variationally optimized operators

## Abstract

We explore the use of ‘optimized’ operators, designed to interpolate only a single meson eigenstate, in three-point correlation functions with a vector-current insertion. These operators are constructed as linear combinations in a large basis of meson interpolating fields using a variational analysis of matrices of two-point correlation functions. After performing such a determination at both zero and non-zero momentum, we compute three-point functions and are able to study radiative transition matrix elements featuring excited state mesons. The required two- and three-point correlation functions are efficiently computed using the distillation framework in which there is a factorization between quark propagation and operator construction, allowing for a large number of meson operators of definite momentum to be considered. We illustrate the method with a calculation using anisotopic lattices having three flavors of dynamical quark all tuned to the physical strange quark mass, considering form-factors and transitions of pseudoscalar and vector meson excitations. The dependence on photon virtuality for a number of form-factors and transitions is extracted and some discussion of excited-state phenomenology is presented.

###### pacs:
12.38.Gc, 13.40.Gp, 13.40.Hq, 14.40.Be
1

## I Introduction

The coupling of mesons and baryons to photons at leading order in is given by matrix elements of the quark-field vector current where . In the case that the initial and final state hadron are the same we express the matrix element in terms of Lorentz-invariant form-factors, whose dependence on the virtuality of the photon, , can be related to quark charge and current distributions within the hadron. The current can also induce a transition from one hadron eigenstate to another , in which case we speak of transition form-factors. For hadrons with non-zero spin, there are multiple possible amplitudes which can be labeled by the helicity of the hadrons, or we may expand the current in terms of multipoles to provide another convenient physically motivated basis.

In the meson sector, these matrix elements appear in photo– and electro–production of mesons from nucleon and nuclear targets, where the coupling of the photon to a -channel meson exchange is described by transition form-factors. Measurements of these processes with unprecedented statistics will be made in the GlueX and CLAS12 detectors at the 12 GeV upgraded CEBAF Dudek et al. (2012a). In particular, photoproduction has been proposed as a means to produce large numbers of exotic hybrids, those mesons which contain an excitation of the gluonic field as well as the usual quark-antiquark pair Horn and Mandula (1978); ?; ?; ?, Dudek (2011).

Production of excited mesons at very forward angles using pion beams is also driven by the vector-current transition form-factor. In the Primakoff process the pion absorbs a nearly on-shell photon from a nucleon or nuclear target, with an amplitude to transition to another meson species given by transition form-factors at . Recent such measurements of the couplings and have been made at COMPASS Adolph et al. (2014).

In charmonium and bottomonium, the relatively small total widths of the low-lying states means that radiative transition rates between them constitute significant branching fractions, and can be measured directly, as can rates of decay to a photon plus light-quark mesons where the heavy quark-antiquark pair annihilates Beringer et al. (2012); Ablikim et al. (2011a); ?; ?; ?.

In the baryon sector, the dependence of transition form-factors can be measured quite directly in electro-production of excited nucleons off proton and neutron targets. The relative magnitudes of the various multipole or helicity amplitudes and the variation with photon virtuality have been discussed as a means to study the internal quark-gluon structure of excited and states Aznauryan et al. (2013).

The properties of hadrons constructed from strongly interacting quarks and gluons should be calculable within the relevant gauge field theory, Quantum Chromodynamics (QCD). At the energy scale of hadrons, QCD does not have a small coupling constant and must be treated non-perturbatively. The tool we will use to achieve this is lattice QCD, in which the field theory is discretised on a finite grid of Euclidean space-time points, and where we can compute correlation functions as an average over a finite but large number of possible gauge-field configurations.

The vector current matrix elements that we are interested in, , appear in three-point correlation functions of generic form, . Here , are operators constructed from quark and gluon fields, capable of interpolating mesons from the vacuum. In general, such hadron operators with definite quantum numbers interpolate not just one QCD eigenstate, but rather a linear superposition of all states with those quantum numbers. Each state propagates through Euclidean time with a factor such that at large times only the state of lowest energy survives, suggesting that if we separate the three operators in the correlation function by large time intervals we will obtain the transition between the lightest states with the quantum numbers of , . However in many cases we actually wish to study states which are not the lightest with a given set of quantum numbers. If we use a generic operator to interpolate them from the vacuum their contribution will be as subleading exponential dependences in the correlation function, which become dominated by lighter states as we propagate in Euclidean time. Determining the amplitude of subleading exponential contributions through fitting the time-dependence proves to be unreliable.

Our solution to this problem is to form ‘optimal’ operators for interpolation of each state that we wish to study, that is operators which have a dominant amplitude to produce a particular state, and significantly reduced amplitude to produce all other states, particularly those lighter than the state in question. In this case the three-point correlation functions are dominated by the desired initial and final states, even if they are not the lightest in a given quantum number sector. Optimal operators can be constructed as a linear superposition of operators in a large basis, where there is a particular linear superposition for each state in the tower of excited states. The different superpositions are orthogonal in a suitable sense. We are able to find appropriate superpositions through a variational analysis of the matrix of two-point correlation functions, , for a set of operators .

The technology we will utilize, constructing ‘optimized’ operators as linear superpositions of a large basis of interpolators, is useful not only for extracting transitions featuring excited states, but also to improve ground-state signals by significantly reducing the unwanted contributions of excited states to correlation functions. The presence of such contributions remains a problem for calculations attempting precision extraction of ground-state matrix elements Aznauryan et al. (2013); Alexandrou et al. (2011); Horsley et al. (2014); Edwards et al. (2006); Lin (2012); Bali et al. (2014); Dinter et al. (2011); Green et al. (2014). Increase in statistical noise precludes separating operators by large time separations, so instead the use of ‘optimized’ operators seeks to deal with the problem directly by suppressing creation of the unwanted excited states.

The use of optimized operators as a tool to extract excited-state radiative transitions was previously considered in Dudek et al. (2009a), where the meson source operators, were optimized, while the meson sink operators, , were simple local fermion bilinears. Quark propagation from the meson sink proceeded using the ‘sequential source’ method, which proved to be a significant limitation on what could be achieved. In this study we will make use of the distillation framework for correlator construction Peardon et al. (2009) within which we will be able to use optimized operators of definite momentum at both source and sink as well as to insert a vector current operator of definite momentum. This is the first use of distillation to compute three-point functions.

In order to demonstrate the technology we perform a calculation on dynamical lattices having three flavors of quark all tuned to approximately the physical strange quark mass - the spectrum of isovector mesons on these lattices was previously presented in Dudek et al. (2009b, 2010). We consider pseudoscalar and vector mesons in the SU(3) flavor octet, extracting form-factors and transition form-factors for both ground and excited states in these channels. This should be expected to be a challenging undertaking – transitions between excited-state and ground-state vector and pseudoscalar mesons are ‘hindered’ and have relatively small magnitude for small photon virtualities – as such we will be extracting small signals from correlation functions built using ‘optimal’ excited state operators.

In section II we introduce the decomposition of vector current matrix elements in terms of Lorentz covariant kinematic factors multiplying the unknown form-factors we seek to extract. We proceed to introduce the gauge configurations, operator basis, and background information relevant to two-point calculations in Sections III and IV. Section V concerns three-point functions; we describe their calculation within the distillation framework and their relation to the matrix elements of interest while also demonstrating the efficacy of optimized operators. We present the results of our calculation in Section VI, comparing with relevant previous calculations, and then conclude with a summary and outlook in Section VII. Appendices describing helicity operators, momentum conservation in a finite-volume, the improvement of the vector current and the flavor structure of the current follow.

## Ii Form-factors and Transitions

The photon couples to the electric charges of quarks via the vector current , up to a factor of the magnitude of the electron charge, . In general a transition induced by this current between a hadron, , of spin- and a hadron of spin- is described by the matrix element,

 Missing or unrecognized delimiter for \big

where the spin-state of is specified in terms of its helicity, , the projection of along the direction of momentum . These matrix elements are simply related to the helicity amplitude for the transition by including the initial-state photon’s polarization vector,

 Missing or unrecognized delimiter for \big (1)

where and where the photon has a virtuality .

There are relations between these amplitudes which follow from the constraints of Lorentz invariance, current conservation and invariance under parity transformations. These can be accounted for if we write a matrix-element decomposition in terms of a number of Lorentz invariant form-factors, ,

 Missing or unrecognized delimiter for \big (2)

The Lorentz covariant ‘kinematic factors’, , are constructed from the meson four-momenta, , , and initial and final state polarization tensors relevant to the spin of the mesons, , . For any given pair of mesons , of definite spin and parity, there are only a limited number of possible constructions consistent with parity invariance and with the additional constraint of current conservation, we can write explicit decompositions in terms of a few independent form-factors.

For example, a pseudoscalar particle like the pion has only a single form-factor appearing in its decomposition,

 ⟨π+(→p′)∣∣jμ∣∣π+(→p)⟩=(p+p′)μFπ(Q2). (3)

At , the vector current measures the charge of the pion in units of , so exactly.

In a transition between two different pseudoscalar particles, there is again only one form-factor, but the kinematic factor differs owing to the differing masses () of the pseudoscalar particles (),

 ⟨π′+ (→p′)∣∣jμ∣∣π+(→p)⟩ =[(p+p′)μ+m′2−m2Q2(p′−p)μ]Fπ′π(Q2). (4)

The transition matrix-element between a vector particle and a pseudoscalar can be expressed as

 ⟨π+(→p′)∣∣jμ∣∣ρ+(λ,→p)⟩ =ϵμνρσp′νpρϵσ(λ,→p)2mπ+mρFρπ(Q2), (5)

and for a vector meson stable under the strong interactions, the transition form-factor at can be related to the radiative decay width , where is the momentum of the final-state photon in the rest-frame of the decaying meson.

The vector current matrix element for a stable vector hadron of mass has three possible covariant structures having the right parity transformation properties once current conservation is demanded Arnold et al. (1980),

 ⟨ρ +(λ′,→p′)∣∣jμ∣∣ρ+(λ,→p)⟩ = −[(p+p′)μϵ∗(λ′,→p′)⋅ϵ(λ,→p)]G1(Q2) +[ϵμ(λ,→p)ϵ∗(λ′,→p′)⋅p+ϵμ∗(λ′,→p′)ϵ(λ,→p)⋅p′]G2(Q2) −[(p+p′)μϵ∗(λ′,→p′)⋅pϵ(λ,→p)⋅p′12m2]G3(Q2), (6)

with a corresponding set of three independent dimensionless form-factors , , . A convenient basis having a clearer physical motivation is provided by the expansion of the vector current in terms of multipoles Durand et al. (1962), which in this case leads to a set of form-factors,

 GC =(1+Q26m2)G1−Q26m2G2+Q26m2(1+Q24m2)G3 GM =G2 GQ =G1−G2+(1+Q24m2)G3, (7)

which are proportional to the charge (), magnetic dipole (), and quadrupole () multipoles respectively2. At they are related to the charge, magnetic moment and quadrupole moment of the vector meson: , , .

The other form-factors we considered above may also be identified with a particular multipolarity – in the transition case the single form-factor is of magnetic dipole () type, while for the cases it is a charge form-factor . For stable meson states, time-reversal invariance indicates that the form-factors are real functions of .

## Iii Calculation Details

In this first investigation of the extraction of excited state form-factors using distillation, we restrict ourselves to a single ensemble of gauge-field configurations, having three degenerate flavors of dynamical quarks tuned to approximately the physical strange quark mass. This set of anisotropic Clover lattices Edwards et al. (2008); Lin et al. (2009) has been used previously in studies of the meson spectrum Dudek et al. (2009b, 2010, 2011a); Liu et al. (2012); Dudek et al. (2013a), meson decay constants Mastropas and Richards (2014), baryon spectrum Edwards et al. (2011); Dudek and Edwards (2012); Edwards et al. (2013); Padmanath et al. (2014) and meson-meson scattering Dudek et al. (2011b, 2012b, 2013b, 2014). For the calculations reported on in this paper, we used 535 configurations of lattice volume , with a spatial grid spacing of and a temporal spacing roughly 3.5 times smaller.

In this calculation we have an exact flavor symmetry such that all the octet mesons (, , ) are degenerate with a mass close to MeV. Where results are expressed in dimensionful units, they are determined from the dimensionless quantities using the scale-setting procedure,

 E=atEatmΩ⋅mphys.Ω.

where is the baryon mass calculated on this lattice and is the experimental value Beringer et al. (2012).

Our use of a Clover-improved anisotropic quark action introduces an improvement term into the vector current which appears at tree-level. Discussion of the effect of improvement and renormalisation of the vector current will appear in Section V.5.

## Iv Two-Point Function Analysis

In order to determine radiative transition amplitudes between meson states within QCD we must first obtain the spectrum of states and find operators, constructed from quark and gluon fields, that reliably interpolate the states of interest from the vacuum. In general a color-singlet operator having definite can produce all QCD eigenstates having those quantum numbers,

 O†i|0⟩=∑n|n⟩⟨n|O†i|0⟩.

We seek to determine optimized interpolators, , which when acting on the vacuum strongly interpolate only a single state with much reduced contributions from other states,

 Ω†n|0⟩ =|n⟩⟨n|Ω†n|0⟩+∑m≠n|m⟩⟨m|Ω†n|0⟩ =|n⟩⟨n|Ω†n|0⟩+∑m≠n|m⟩εm.

In essence we seek a procedure by which we can minimize the () relative to the strength with which our operator creates the ’th state, .

We will proceed by using a basis of interpolators, , to construct two-point correlation functions of the form,

 Cij(t)=⟨0|Oi(t)O†j(0)|0⟩,

where operators are color-singlet constructions built from the basic quark and gluon fields of QCD, having the quantum numbers of the desired hadrons. Such correlation functions can be expressed as,

 Cij(t)=∑n12En⟨0|Oi(0)|n⟩⟨n|O†j(0)|0⟩e−Ent,

where the spectrum of eigenstates is seen to control the Euclidean time dependence3.

### iv.1 Variational Analysis

We propose that within any basis of operators there is a particular linear combination that is most suited to interpolate the lightest state of the spectrum, another linear combination that optimally interpolates the first excited state, a third combination for the second excited state and so on. Thus optimized interpolators take the form , where one can show that the best estimate for the weights , in a variational sense, comes from solving the generalized eigenvalue problem Michael (1985); Luscher and Wolff (1990); Dudek et al. (2008); Blossier et al. (2009),

 C(t)v(n)=λn(t)C(t0)v(n). (8)

Here is the matrix whose elements are the correlation functions constructed from the basis of operators, , and is a generalized eigenvector. The generalized eigenvalues, or principal correlators, behave like at large times, and can be used to determine the spectrum of energy eigenstates. The vectors are orthogonal on a metric, , where is a reference timeslice. Examination of the orthogonality condition suggests that should be chosen to be sufficiently large such that the correlation functions are dominated by the lowest-lying states, with heavier states having decayed exponentially to a negligible level. Further considerations on the choice of are presented in Dudek et al. (2008); Blossier et al. (2009).

In practice we solve Eqn 8 independently on each timeslice, , so that for each state, , we obtain a time series of generalized eigenvectors , which we observe to be essentially time-independent with a suitably large choice of . In practice we use the mean values (over the ensemble of gauge configurations) of the elements chosen on a single timeslice to construct the optimized operators as

 Ω†n=√2Ene−Ent0/2∑iv(n)iO†i, (9)

where the coefficients multiplying the sum are chosen to give the normalization , which will prove to be convenient when considering three-point functions.

The procedure of solving the generalized eigenvalue problem in a basis of interpolating fields is carried out independently for each quantum number channel at each possible allowed momentum value.

### iv.2 Meson operator construction

A straightforward approach to constructing a basis of operators capable of interpolating meson states is to use fermion bilinears containing some number of spatially directed gauge-covariant derivatives, that is operators of generic structure,

 O∼¯ψΓ↔D⋯↔Dψ, (10)

where . By expressing the vector-like derivatives and Gamma matrices in a circular basis, we can easily construct operators of definite spin using the standard Clebsch-Gordon coefficients Dudek et al. (2009b, 2010). Operators of this type with non-zero momentum can be constructed to have definite helicity, as described in Thomas et al. (2012).

In this calculation, QCD is discretised on a grid of points whose spatial structure has a cubic symmetry and as such we have not the complete continuous rotational symmetry of the continuum, but rather a reduced symmetry: the symmetry of the cubic group at rest and the relevant little group in a moving frame. A consequence is that instead of having an infinite number of irreducible representations labeled by integer spin (at rest) we only have access to the finite number of irreducible representations, or irreps, of the cube, labelled , , , , and . A corresponding argument applies in-flight, where the continuum helicity labeling is broken down to a finite number of irreps of the little group, see Thomas et al. (2012) for details.

The operators of definite (or helicity) constructed above can be projected into irreps of the relevant symmetry group using a procedure called subduction,

 O[J]Λμ=J∑M=−JSJMΛμOJM (11)

where labels the cubic irrep and is the ‘row’ of the irrep (. The subduction coefficients, are tabulated in Thomas et al. (2012); Dudek et al. (2010). These operators have been used extensively to study the excited spectrum of mesons Dudek et al. (2013a); Liu et al. (2012); Dudek et al. (2011a, 2009b, 2010, 2013b, 2014).

### iv.3 Correlator construction through distillation

Operators which interpolate hadrons from the vacuum have long been known to do so more effectively if the quark fields are suitably smeared over space Gusken et al. (1989); Allton et al. (1993). An extremely convenient method to do this is provided by distillation Peardon et al. (2009), where the smearing operator is constructed on time-slice as an outer product of vectors in color and space,

 □→x→y(t)=ND∑n=1ξ(n)→x(t)ξ(n)†→y(t), (12)

where the vectors should be constructed to have strong overlap onto the low-energy quark modes most relevant to low-lying hadron states. A suitable choice for the vectors are the eigenvectors of the gauge-covariant three-dimensional Laplacian on a time-slice ordered by their eigenvalue.

Smearing each quark field, a meson creation operator of fermion bilinear form with momentum is

 O†(→p)=¯ψ→x□→x→ye−i→p⋅→yΓ→y→z□→z→wψ→w

where time, color, and spin indices have been suppressed for brevity and where repeated position indices are summed. The object can be non-local in -space, and may for example feature gauge-covariant derivatives as discussed in the previous section.

An advantage of the distillation framework is that it leads to a factorization of two-point functions into matrices describing quark propagation, called perambulators, and matrices describing operator construction. A generic connected two-point function using fermion bilinear constructions, in which we explicitly show the smearing operators, can be decomposed as

 ⟨0∣∣ ¯ψ□Γf□ψ(t)¯ψ□Γi□ψ(0)∣∣0⟩ =−τnm(0,t)Φfmp(t)τpq(t,0)Φiqn(0)

where the perambulators, , can be obtained by inverting the Dirac matrix on sources at timeslice 0, and where encodes the operator construction (here we include the momentum projection in ).

### iv.4 Meson spectra & optimized operators

We computed correlation matrices for irreps corresponding to at rest, and to magnitude of helicity with nonzero momentum. We made use of , where the momentum is expressed in units of , . The quark flavor constructions were chosen to give access to the members of the octet – in this case the two-point function Wick contraction contains only a single connected diagram.

In the rest frame we used operator constructions of the type in Eq. 10 with up to three derivatives while for nonzero momentum we used up to two-derivative constructions. Further details about the construction of derivative operators and the basis used in this calculation can be found in Thomas et al. (2012); Dudek et al. (2009b, 2010, 2013a).

The resulting correlation matrices were analyzed using Eq. 8, with each principal correlator, , being fit with the form,

 λn(t)=(1−An)e−En(t−t0)+Ane−E′n(t−t0), (13)

where the fit parameters are , , and – the second exponential is present to absorb the effect of any states other than remaining in the principal correlator. In practice, for suitably large values of , we find that the energy scale of is typically at or above the largest energy extracted, .

A typical example is presented in Figure 1, where the principal correlators for the lightest three states in the , irreps, which contain mesons, are shown. We make use of a basis of 26 operators in the channel and 12 operators for .

Another example is shown in Figure 2(a) for the case , , which contains the helicity zero component of mesons of ‘unnatural parity’ (), where we present the effective mass of the principal correlator, with . As will be discussed in the next subsection, these states appear to correspond (in order of increasing energy) to the ground-state , the , the first-excited , and the .

We note that in both figures we observe that there is negligible curvature present for , indicating that a single state is dominating the correlation function. This observation only holds for sufficiently large choices of , and is one guide in the selection of a suitable value. The time beyond which excited state contributions are negligible, which we can call , plays an important role in the construction of our three-point functions in terms of selecting the time separation of the source and sink meson operators. We desire our three-point functions to feature a time region in which the correlation function is dominated by the transition of interest, with contributions from other states having decayed away. In order to achieve this we should separate the source and sink projected operators (Eq. 9) by at least . It remains possible that the vector current insertion may act to suppress or amplify the contribution of unwanted excited states, but since we do not have this information in advance of the calculation, the above time separation serves as a reasonable estimation of the minimum.

The source and sink cannot be separated arbitrarily far as the statistical noise on the entire three-point correlation function grows exponentially with increasing separation. We can obtain an estimate of the maximum practical time separation for three-point functions by examining the growth of noise in two-point function principal correlators. As an example in Figure 1 we see that the ground-state signal remains of high statistical quality out to 25 timeslices (and beyond), while the first and second excited states begin to show significant fluctuations above . Later we will find that while optimized ground-state three-point correlation functions are still statistically precise for time separations as high as 36 timeslices, excited-state correlation functions are not well determined for separations larger than around 20 timeslices.

In practice we solve the matrix problem, Eq. 8, independently on each timeslice (as described in Dudek et al. (2010)). For sufficiently large we find that the elements of the eigenvectors so obtained, , are essentially flat for , and in practice we construct our projected operators, Eq. 9, using ensemble mean values taken from a single timeslice. In Figure 2(b) we show the correlation functions , observing that they behave in the manner we expect for optimized operators. The corresponding off-diagonal correlation functions for are statistically compatible with zero for .

We conclude this section by demonstrating that an optimized ground-state operator constructed as described above, does indeed significantly reduce the contribution of excited states to a correlation function, relaxing to the ground-state more quickly. This can be seen in Figure 3 where the optimized ground-state operator in is compared with the (distillation smeared) operator.

### iv.5 Meson dispersion relations

As previously mentioned, we independently compute the energy spectrum of states at each allowed value of for each relevant irrep, but we expect to see the same mesons appearing at each momentum with an energy determined by their rest mass and the relevant dispersion relation. After identifying the mesons at each momentum (by their overlap with characteristic operators Thomas et al. (2012)), we may examine their dispersion relation, . This is presented in Figure 4 for mesons, and mesons (not used in this analysis) where they are all observed to be compatible with the relativistic dispersion relation, , or in temporal lattice units,

 (atE)2=(atm)2+(2πξ(L/as))2|→n→p|2,

for a meson of mass , with momentum , with the anisotropy taking the value . Making use of optimized operators with momenta up to allows us to sample many values of in the form-factor extraction.

## V Three-point functions

We now turn to the three-point correlation functions used in this analysis which contain the vector-current matrix elements of interest. Their basic form is

 Missing or unrecognized delimiter for \big

where the operators are capable of interpolating meson states of definite momentum from the vacuum – a suitable basis was discussed in the previous section. The relation between the correlation function and the desired matrix-element is exposed by inserting complete sets of eigenstates with the quantum numbers of and , and evolving all operators back to the origin of Euclidean time,

 Cfμi(Δt,t)=∑ni,nf 12Enf12Enie−Enf(Δt−t)e−Enit ×⟨0∣∣Of(0)∣∣nf⟩⟨nf∣∣jμ(0)∣∣ni⟩⟨ni∣∣O†i(0)∣∣0⟩.

The summation runs over all states, but clearly if the separations between the operators are large, , only the lightest states in the and channels will contribute and we can extract the vector-current matrix element between them. However, at modest time separations there will remain subleading exponential contributions from excited states, and these ‘polluting’ terms can be a source of systematic error in the extraction of ground state matrix elements Aznauryan et al. (2013); Alexandrou et al. (2011); Horsley et al. (2014); Edwards et al. (2006); Yamazaki et al. (2008); Lin (2012); Bali et al. (2014); Dinter et al. (2011); Green et al. (2014). Reducing excited state pollution by simply separating operators by longer Euclidean times is not always practical, due to the increase in statistical noise with increasing separation.

One of our major aims here is to extract excited-state matrix elements, which we may access using the optimized operators described in the previous section. Using the optimized operator for an excited state should lead to a three-point correlation function whose leading behavior at large times is given not by the ground state, but rather by the relevant excited state. If we are interested in the ground-state there is also an advantage to using the appropriate optimized operator in that it will have much reduced overlap onto low-lying excitations (relative to any single operator in the original basis, for example , c.f. Figure 3), leading to a corresponding reduction in the excited state pollution in the three-point correlator.

Three-point correlation functions using optimized operators take the form,

 Cnfμni(Δt,t) =⟨0∣∣Ωnf(Δt)jμ(t)Ω†ni(0)∣∣0⟩ Extra open brace or missing close brace (14)

where the leading time-dependence is that of the states , selected by the choice of optimized operators. The ellipsis represents the residual contributions of other states, which should be significantly suppressed when using optimized operators – we will explore the degree to which this is manifested in explicit calculation.

### v.1 Correlator construction & distillation

Three-point correlation functions featuring a current insertion require a slight extension of the distillation framework presented in Section IV.3 since the quark fields in the current should be those which appear in the action, which are not smeared. Exposing the smearing operators, the general form required is,

 ⟨0∣∣¯ψ□Γf□ψ(Δt)⋅¯ψΓψ(t)⋅¯ψ□Γi□ψ(0)∣∣0⟩,

which, considering for now only the completely connected Wick contraction, can be decomposed as

 ⟨ 0∣∣¯ψ□Γf□ψ(Δt)⋅¯ψΓψ(t)⋅¯ψ□Γi□ψ(0)∣∣0⟩ =−τnm(0,Δt)Φfmp(Δt)[ξ(p)†(Δt)M−1Δt,tΓM−1t,0ξ(q)(0)]Φiqn(0) =−τnm(0,Δt)Φfmp(Δt)GΓpq(Δt,t,0)Φiqn(0),

where the object in square brackets, , is a generalized perambulator. It can be obtained through inversion from sources at timeslice 0 to obtain , inversion from sources at timeslice to obtain which can be related to using hermiticity, and contraction with the operator insertion, , at each value of between and . In this calculation we do not average over multiple time-sources separated by the same value of , although this could be done to increase statistics.

In our application, exposing the Dirac spin, space and color indices, the current insertion is . In the case of an improved vector current, which we consider later, we require also , yielding another set of generalized perambulators.

Within this construction we are able to project each operator into definite momentum, and as such we only compute correlation functions in which the momentum is conserved, . Some discussion of momentum conservation in a finite volume appears in Appendix B.

### v.2 Correlation functions using optimized operators

Turning first to the case of three-point functions with pion-like operators at the source and sink, we plot in Figure 6 the form-factor (as defined in Eqn. 3) extracted from the three-point function,

 ⟨0|Oπ(Δt,→pf)jμ(t,→q)O†π(0,→pi)|0⟩

where represents either (in red) or the optimized operator (in blue). The sink operator, located at , is in the irrep of momentum , while the source operator, located at , is at rest in the irrep. We clearly observe that the optimized operators give rise to a signal which is flat over a number of timeslices away from the source and sink, corresponding to the contribution of just the ground-state pion, while the simpler operators over this time range always retain a non-negligible pollution from excited states. Such behavior is expected from our two-point function analysis: for example, at rest we find , so the distillation smeared operator , acting on the vacuum, creates both the ground and first excited state with comparable strength.

Our principal motivation for using optimized operators is to get access to transitions involving excited states. In Figure 7 we show matrix elements extracted from three-point correlation functions computed using either the ground-state or first-excited state optimized operator at the source () and either the ground-state or first-excited state operator at the sink (). We observe that there are clear statistically significant signals for excited-state transitions when using the appropriate optimized operators.

In general, even for optimized operators, there may still be some residual contamination coming from states that lie beyond the reach of our variational basis, and indeed curvature away from flat behavior as we approach the source or sink timeslice is observed in Figures 6 and 7.

In order to make maximal use of the time-series data, in particular in those regions where there remains some unwanted excited-state contribution, we opt to perform a correlated fit over a time range with the form,

 F(Q2;t)=F(Q2)+ffe−δEf(Δt−t)+fie−δEit (15)

where and are real fit parameters. We make further use only of the constant term, which corresponds to the desired form-factor. Fitting the data to this form also exposes the energy scale of the pollution terms, and . Generically, when present, we find that these energies lie at or above the scale of the largest energies we reliably extract in our two-point function variational analysis. In cases where there is a clear extended plateau region, we may exclude the exponential terms and perform a fit to a constant value.

The dependence upon source-sink separation, , for the ground-state pion form-factor can be seen in Figure 8. The lower panel shows the illustrative case , , where we observe that there is only a visible plateau region for , while for the shorter separations, , we make use of a fit using Eqn. 15. The resulting values of are observed to be compatible – other time separations were also explored with similar results. The upper panel shows that this procedure is generally applicable and leads to form-factors from each time-separation that are in agreement across a range of .

In practice, while we extract a very large number of form-factor determinations at many -values, we choose to make use of only those where application of Eq. 15 to shows modest excited-state contributions. Any cases where a clear trend toward a constant value is not visible are discarded.

### v.3 Extracting multiple form-factors

Eqn. 2 presents the general form of the decomposition of a vector-current matrix element into independent form-factors, , and the corresponding kinematical factors, , which depend upon momenta and helicities. Moving to a more complete notation including momentum and helicity labels, our three-point correlation functions may be written,

 ⟨0∣∣Ωnf,→pf,λf(Δt)jμ→q(t)Ω†ni,→pi,λi(0)∣∣0⟩ =e−Enf(Δt−t)e−Enit⟨nf,→pf,λf∣∣jμ(0)∣∣ni,→pi,λi⟩+… =e−Enf(Δt−t)e−Enit∑iKμi(nf,→pf,λf;ni,→pi,λi)Fi(Q2)+…,

where as previously the ellipsis represents possible pollution from states other than , which will have residual time-dependence, but which as shown in the previous section are suppressed when using optimized operators.

Any one correlator provides, in general, an underdetermined linear system for the multiple form-factors we wish to extract. By using different combinations of initial and final state helicities and momenta all at the same , we can build a linear system which is constrained or over-constrained, from which we can determine the set . Removing the known Euclidean time-dependence, , from the correlation functions, we should be left with objects which are time-independent up to pollution from other states, which should be modest for optimized operators. For fixed state choices, , using an indexing , we can write the linear system,

 Γa(Δt,t) ≡eEnf(Δt−t)eEnit⟨0∣∣Ωnf,→pf,λf(Δt)jμ→q(t)Ω†ni,→pi,λi(0)∣∣0⟩ =∑iKi(a)Fi(Q2)+…, (16)

which is of the form with a vector over , a rectangular matrix with indices , and a vector of form-factors, indexed by . This may be converted into a system featuring a square matrix: , which can be inverted using SVD. In the case where only a single form-factor contributes this procedure can still be followed as a way to average over rotationally equivalent momentum combinations.

In practice we solve this system independently for each value of between and – if our optimized operators were perfect, leading to no pollution from other states, we would obtain the same form-factors on each timeslice. In fact we obtain , where the time-dependence is fitted as described in the previous section to account for the presence of pollution from other states.

This approach was previously used in Dudek et al. (2006) and Dudek et al. (2009a) in the extraction of charmonium form-factors.

### v.4 Cubic symmetry

A consequence of discretizing QCD on a hypercubic grid is that the theory does not possess the full three dimensional rotational symmetry of the continuum. Instead, we are restricted to a subset of rotations that leave the cube invariant. This smaller symmetry group has only a finite number of irreducible representations into which the infinite set of continuum representations labelled by integer spin, , must be subduced. A simple example is , where the five equivalent ‘rows’ () get distributed into a three-dimensional irrep called and a two-dimensional irrep called . Because there are only a finite number of these irreps, they must accommodate multiple values of , such that also contains parts of . For systems with non-zero momentum, the symmetry group is called the ‘little group’, and the corresponding subduction is from helicity, . Tables of the spin/helicity content of cubic irreps can be found in Thomas et al. (2012).

To correctly reflect the symmetry of our theory then, we should label our correlation functions according to irreducible representations of the cubic symmetry. In practice this is what we do by computing using the subduced operators introduced in Section IV.2. Using these operators, the three-point functions take the form,

 ⟨0∣∣ΩΛf,μfnf,→pf(Δt)jΛγ,μγ→q(t)ΩΛi,μi†ni,→pi(0)∣∣0⟩, (17)

where the indices , label the cubic group irrep and the ‘row’ () of the irrep.

Considering only the cubic symmetry of the lattice, and not any underlying continuum-like symmetry, we would not expect there to be any relationship between different irreps. Furthermore matrix element decompositions should be defined in terms of the irreps of the cube, not in terms of hadrons of definite spin. For example a correlation function with a operator at the source should take values that need not be related to one with an operator at the source.

However, were there really to be no relation, we could hardly claim to be approximating QCD in a realistic manner. In practical calculations it should be the case that through a combination of sufficiently fine lattice spacing, reduction of discretization artifacts through improvement of the action Symanzik (1983), and interpolation of hadrons using operators smoothed over many lattice sites Davoudi and Savage (2012), that the continuum symmetry is manifested to a good approximation with only small deviations. For example we might expect to see a relation between and correlation functions corresponding to them originating from the same meson. In previous two-point function calculations we have observed that the rotational symmetry of the continuum theory is clearly visible in relations amongst the irreps both for eigenstate masses and the values of matrix elements Dudek et al. (2009b, 2010).

Since we expect to see a comparable restoration of the rotational symmetry in this calculation, we do not attempt to build decompositions according to the symmetries of the cube, rather making use of the continuum-like helicity decompositions presented earlier, subduced into irreducible representations of the cube.

A slight additional complication in this analysis arises from our use of anisotropic gauge configurations in which the space and time directions are discretized with different spacings. Spatially directed currents will need to be renormalized separately from temporal currents and the discretization effects along the two directions are expected to be different – in explicit calculation we will not mix spatially directed currents with their temporal counterparts. Had we used isotropic lattices the temporal component of the vector current would be related to the spatial components, however here we will keep them separate with the temporal component of the current subducing differently from the spatial components. For spatial components4, the subduced current is