The $\rho$-meson light-cone distribution amplitudes from lattice QCD

We present the results of a lattice study of the normalization constants and second moments of the light-cone distribution amplitudes of longitudinally and transversely polarized $\rho$ mesons. The calculation is performed using two flavors of dynamical clover fermions at lattice spacings between $0.060\,\text{fm}$ and $0.081\,\text{fm}$, different lattice volumes up to $m_\pi L = 6.7$ and pion masses down to $m_\pi=150\,\text{MeV}$. Bare lattice results are renormalized non-perturbatively using a variant of the RI'-MOM scheme and converted to the $\overline{\text{MS}}$ scheme. The necessary conversion coefficients, which are not available in the literature, are calculated. The chiral extrapolation for the relevant decay constants is worked out in detail. We obtain for the ratio of the tensor and vector coupling constants $f_\rho^T/f_\rho^{\vphantom{T}} = 0.629(8)$ and the values of the second Gegenbauer moments $a_2^\parallel = 0.132(27)$ and $a_2^\perp = 0.101(22)$ at the scale $\mu = 2\,\text{GeV}$ for the longitudinally and transversely polarized $\rho$ mesons, respectively. The errors include the statistical uncertainty and estimates of the systematics arising from renormalization. Discretization errors cannot be estimated reliably and are not included. In this calculation the possibility of $\rho\to\pi\pi$ decay at the smaller pion masses is not taken into account.


Introduction
In recent years exclusive reactions with a large momentum transfer to a light vector meson V = ρ, K * , φ in the final state are attracting increasing attention. Prominent examples are provided by B-meson weak decays, B → V π, B → V ν , B → V γ, B → V µ + µ − , etc. Their study constitutes a considerable part of the experimental program of the LHCb collaboration at CERN [1] and the future Belle II experiment at the upgraded KEK facility [2]. Among these processes, the decays B → K * µ + µ − and B s → φµ + µ − are of particular relevance as the angular distributions of the decay products give access to a host of observables that are sensitive to new physics, see, e.g., ref. [3] for a recent review. Another example is deeply-virtual exclusive ρ-meson production (DVMP) in electron-nucleon collisions at high energy, eN → eρN , that, besides deeply-virtual Compton scattering (DVCS), allows one to resolve the transverse distribution of partons inside the nucleon. The corresponding cross sections were measured by the HERA collider experiments H1 and ZEUS and the fixed target experiments HERMES (DESY), CLAS (JLAB), and Hall A (JLAB) at small and moderate values of the Bjorken momentum fraction x Bj , respectively. In the future, exclusive vector meson production will be studied with unprecedented precision at the electron-ion collider (EIC) [4].
The standard framework for the theoretical description of such processes is based on collinear factorization. In this approach the vector mesons are described in terms of lightcone distribution amplitudes (DAs) that specify the distribution of the longitudinal momentum amongst the quark and antiquark in the valence component of the wave function; the transverse degrees of freedom are integrated out. In general, meson and baryon DAs are scale-dependent non-perturbative functions and their moments (weighted integrals over the momentum fractions) are given by matrix elements of local operators. From the phenomenological point of view the normalization (representing the value of the wave function at the origin) and the first non-trivial Gegenbauer moment that characterizes the width of the DA are the most relevant quantities. For example, knowledge of the second moment of the DA of the longitudinally polarized ρ meson is crucial for global fits of generalized parton distributions from the DVMP and DVCS data [5].
The ρ-meson coupling to the vector current is known experimentally and the other parameters were estimated in the past using QCD sum rules [6], see also ref. [7] for an update. Lattice calculations of the tensor coupling have been reported in refs. [8][9][10][11][12] and the second moments in ref. [13].
In this work we present new results using two flavors of dynamical clover fermions at lattice spacings between 0.060 fm and 0.081 fm, different lattice volumes and pion masses down to m π = 150 MeV. Our approach is similar to the strategy used in our paper on the pion DA [14]. In addition to a much larger set of lattices as compared to the previous studies, a new element of our analysis is a consistent use of non-perturbative renormalization including mixing with the operators containing total derivatives. As the coefficients for the conversion between our non-perturbative renormalization scheme on the lattice and the MS scheme are not available in the literature for tensor operators, we have performed the necessary calculations in continuum perturbation theory to two loop accuracy. The chiral extrapolation for the relevant quantities is worked out in detail.
Although our calculation presents a considerable improvement as compared to earlier studies, there are still several issues that we do not address in this work. First and foremost, we only consider ρ mesons and leave the effects of the SU (3) flavor breaking for a future study. Likewise, we do not consider ω mesons that would require the calculation of disconnected diagrams and different techniques. We also do not attempt to take into account effects due to the ρ → ππ decay that becomes possible at the smaller pion masses used in our simulations, although for the short-distance observables considered in this work it seems unlikely that such effects are of principal importance. Last but not least, discretization errors cannot be estimated reliably using the set of lattices at our disposal, which may be an important problem in such calculations. We expect to be able to improve on some of these issues using new N f = 2 + 1 lattice configurations that are being generated in the framework of the CLS initiative [15]. This work, aiming in the long run at smaller lattice spacings with the help of open boundary conditions, is in progress.
The presentation is organized as follows. Section 2 is introductory, we collect the necessary definitions and specify the quantities that will be considered in this work. Section 3 contains a list of the correlation functions that we compute on the lattice. The lattice ensembles at our disposal and the procedure used to extract the signal are described in section 4. A non-perturbative calculation of the necessary renormalization factors is described in section 5, supplemented by appendix A, where we consider the renormalization of the same operators in the continuum and sketch a two loop calculation of the corresponding conversion factors. Complete results needed for the evaluation of the matching coefficients between our RI ′ -SMOM scheme (defined as in ref. [14]) and the MS scheme are presented in the auxiliary file attached to the electronic version of this paper. Section 6 is devoted to the data analysis and the extrapolation to the physical pion mass using, where available, chiral effective field theory expressions derived in appendix B. The final section 7 contains a summary of our results and a discussion.

Continuum formulation
The ρ meson has two independent leading twist (twist two) DAs, φ ρ and φ ρ [16], corresponding to longitudinal and transverse polarization, respectively. Neglecting isospin breaking and electromagnetic effects, the DAs of charged ρ ± and neutral ρ 0 mesons are related so that it is sufficient to consider one of them, for example, ρ + . The DAs are defined as meson-to-vacuum matrix elements of renormalized non-local quark-antiquark light-ray operators, ⟨0 d (z 1 n) n[z 1 n, z 2 n]u(z 2 n) ρ + (p, λ)⟩ = m ρ f ρ (e (λ) · n) 1 0 dx e −ip · n(z 1 (1−x)+z 2 x) φ ρ (x, µ) , ,µ n ν ⟨0 d (z 1 n)σ µν [z 1 n,z 2 n]u(z 2 n) ρ + (p, λ)⟩ = if T ρ e (λ ′ ) · e (λ) (p · n) 1 0 dx e −ip · n(z 1 (1−x)+z 2 x) φ ρ (x, µ) , (2.1b) where σ µν = i 2 [γ µ , γ ν ], z 1,2 are real numbers, n µ is an auxiliary light-like vector (n 2 = 0), and ρ + (p, λ)⟩ is the state of the ρ + meson with on-shell momentum p 2 = m 2 ρ and polarization λ. The straight-line path-ordered Wilson line connecting the quark fields, [z 1 n, z 2 n], is inserted to ensure gauge invariance. The ρ-meson polarization vector e (λ) µ has the following properties: 2) and we use the notation 3) The variable x has the meaning of the fraction of the ρ meson's light-cone momentum p · n which is carried by the u-quark, whereas 1 − x is the momentum fraction carried by the antiquarkd, and µ is the renormalization scale (we assume the MS scheme). The scale dependence will often be suppressed in what follows. The couplings f ρ and f T ρ appearing in (2.1) are defined as matrix elements of local operators: In the following, we will refer to them as vector and tensor couplings, respectively. The vector coupling f ρ is scale independent and can be extracted from experiment, see appendix C in ref. [17] for a detailed discussion. One obtains [17] where for the neutral ρ meson we quote separate values for theūu anddd currents. The difference in the given three values is due to isospin breaking and electromagnetic corrections, which will be neglected throughout this study.
The tensor coupling f T ρ is scale dependent and is not directly accessible from experiment. To leading order one obtains is the number of colors and N f the number of active flavors.
The DAs are normalized to unity, and, neglecting isospin breaking effects, are symmetric under the interchange of the momentum fractions of the quark and the antiquark, For convenience we introduce a generic notation ⟨⋯⟩ , for the moments of the DAs defined as weighted integrals of the type The symmetry property (2.8) implies and in addition we have the (momentum conservation) constraint Hence the set of moments (2.9) for positive integers k, l is overcomplete. We introduce the variable corresponding to the difference of the momentum fraction between the quark and the antiquark and consider ξ moments ⟨ξ n ⟩ , = ⟨(2x − 1) n ⟩ , , n = 2, 4, 6, . . . , (2.13) or, alternatively, Gegenbauer moments as independent non-perturbative parameters. The two sets are related by simple algebraic relations, e.g., The rationale for using Gegenbauer moments is that they have autonomous scale dependence at the one loop level a n (µ) = a n (µ 0 ) . (2.16b) The anomalous dimensions are given by As Gegenbauer polynomials form a complete set of functions, the DAs can be written as an expansion Typical integrals that one encounters in applications can also be expressed in terms of the Gegenbauer coefficients, e.g., Since the anomalous dimensions increase with n, the higher-order contributions in the Gegenbauer expansion are suppressed at large scales so that asymptotically only the leading term survives, usually referred to as the asymptotic DA: Beyond the leading order, higher Gegenbauer coefficients a n mix with the lower ones, a k , k < n [18,19]. This implies, in particular, that Gegenbauer coefficients with higher values of n are generated by the evolution even if they vanish at a low reference scale. This effect is numerically small, however, so that it is usually reasonable to employ the Gegenbauer expansion to some fixed order.

Lattice formulation
From now on we work in Euclidean space, using the same conventions as in ref. [14]. The renormalized light-ray operators entering the definition of the DAs are defined as the generating functions for the corresponding renormalized local operators, cf. ref. [20]. This means that moments of the DAs, by construction, are given by matrix elements of local operators and can be evaluated on the lattice using the Euclidean version of QCD. Our aim in this work is to calculate the couplings f ρ , f T ρ and the second DA moments. To this end we define bare operators On the lattice the covariant derivatives will be replaced by their discretized versions. Projection onto the leading twist corresponds to symmetrization over the maximal possible set of Lorentz indices and subtraction of traces. The operation of symmetrization and trace subtraction will be indicated by enclosing the involved Lorentz indices in parentheses, for instance, Note that for the operators involving the σ µν -matrix also those traces have to be subtracted which correspond to index pairs where one of the indices equals µ or ν.
Using the shorthand ⃗ and its matrix element between the vacuum and the ρ state is proportional to the bare value of the second moment ⟨ξ 2 ⟩ : 24) where N (µνρ) is a kinematical prefactor. The operator V + (µνρ) (x) in the continuum reduces to the second derivative of the vector current, so that with the same prefactor. While in the continuum ⟨1 2 ⟩ bare = 1 by construction, this is no longer true on the lattice because the Leibniz rule holds for discretized derivatives only up to lattice artefacts and hence (2.25) is violated. As we will see below, the deviation from unity for the renormalized ⟨1 2 ⟩ is small. Nevertheless, it still has to be taken into account and affects the relation between ⟨ξ 2 ⟩ and the Gegenbauer moment at finite lattice spacing [14]: The situation with the tensor operators T ± µ(νρσ) and the corresponding matrix elements ⟨⋯⟩ is similar.
The operators V − (µνρ) and V + (µνρ) mix under renormalization even in the continuum, as do T − µ(νρσ) and T + µ(νρσ) . Additional mixing could result from the fact that the continuous O(4) symmetry of Euclidean space is reduced to the discrete H(4) symmetry of the hypercubic lattice. This is particularly worrisome if operators of lower dimension are involved. Fortunately, in the case at hand it is possible to avoid additional mixing by using suitably chosen operators, which will be detailed below.

Lattice correlation functions
The basic objects from which moments of the ρ DAs can be extracted on the lattice are 2-point correlation functions. In order to "create" the ρ meson we use the interpolating current V ν (x), which is defined as V ν (x) with smeared quark fields. For details of our smearing algorithm see section 4. Let O be a local (unsmeared) operator, e.g., one of the operators defined in eq. (2.22) above. One then obtains for the 2-point function in the region where the ground state dominates Here T is the time extent of the lattice, a is the lattice spacing, and E denotes the energy of the ρ state. The sign factors σ are determined by the Dirac matrices in the creation operator (which is in our case always γ ν ), while σ O are the analogous factors for O (see  table 1), and n t is the number of time derivatives in O.
For the decay constants and the second DA moments of the ρ meson we have to evaluate the following set of correlation functions:

Decay constants
In order to determine the leading twist ρ-meson couplings we use the correlation functions with j = 1, 2, 3, assuming the dominance of the lowest one-particle state.
In the actual fits we average over the forward and backward running states. As in our simulations the signal disappears in the noise well before the middle point t = T 2 in the time direction is reached (see figure 1 for an example), the "mixing" of these two contributions is completely negligible. Therefore we work with simple exponential fits, where the averaging operatort ± is defined aŝ The decay constants f ρ and f T ρ can be obtained by simultaneously fitting the correlation functions (3.5a)-(3.5c). The result for the mass is then dominated by the two correlation functions (3.5b) and (3.5c) that contain an unsmeared operator at the sink, because they have much smaller statistical errors. However, they exhibit larger contributions from excited states so that the isolation of the ground state is less reliable. Therefore we first fit the correlator with a smeared operator at the sink, (3.5a), to extract Z ρ and m ρ . These values are then inserted in eqs. (3.5b) and (3.5c) in order to obtain f ρ and f T ρ as well as f T ρ f ρ from a second fit. This procedure is repeated on every bootstrap sample allowing an estimation of the statistical error.

Second moments -the longitudinal case
Multiplets of twist-2 operators suitable for the evaluation of the second longitudinal moments consist of the operators In order to be able to write these and some of the following formulae in a compact form we have introduced the notation R ± , where R + (R − ) is the bare value of ⟨1 2 ⟩ (⟨ξ 2 ⟩ ). We will try to increase the signal-to-noise ratio by considering only correlation functions with the smallest non-zero momentum in one spatial direction, which is equal to 2π L on a lattice of spatial extent L. Therefore we exclude O ± 4 from our calculation. After averaging over all suitable combinations as well as over forward and backward running states, the second longitudinal moments can be obtained from the ratio where momentum averaging is accounted for by the operatorp ± :

Second moments -the transverse case
In the transverse case we consider the following multiplets: As in the longitudinal case, we only consider correlation functions with the smallest non-zero momentum in one spatial direction and perform averages similar to those in eq. (3.9). This leads to the following ratio for the second transverse moments:  Table 2. Ensembles used for this work. For each ensemble we give the inverse coupling β, the hopping parameter κ, the pion mass m π , the finite volume corrected pion mass m ∞ π determined in ref. [22], the lattice size, the value of m π L, where L is the spatial lattice extent, the number of configurations N conf and the number of sources N src used on each configuration. Note that the pion masses have been slightly updated compared to the numbers in ref. [14]. The ensembles marked with † were generated on the QPACE systems of the SFB/TRR 55, while the others were generated earlier within the QCDSF collaboration.

Details of the lattice simulations
For this work we used gauge configurations which have been generated using the Wilson gauge action with N f = 2 flavors of non-perturbatively order a improved Wilson (clover) fermions. A list of the ensembles used is shown in table 2. We used lattices with three different inverse couplings β = 5.20, 5.29, 5.40, which correspond to lattice spacings between 0.06 fm and 0.081 fm. The pion masses vary between 150 MeV and 500 MeV, with spatial volumes between (1.9 fm) 3 and (4.5 fm) 3 .
In order to increase the overall statistics we performed multiple measurements per configuration. The source positions of these measurements were selected randomly to reduce the autocorrelations. To obtain a better overlap with the ground state we applied Wuppertal smearing [23] in the interpolating current V ν using APE smeared gauge links [24].
For the statistical analysis we generate 1000 bootstrap samples per ensemble using a binsize of 4 to further eliminate autocorrelations. For the purpose of maximizing the statistics of the second moments, we average for each bootstrap sample over all suitable combinations of 2-point functions, all possible momentum directions as well as over forward and backward running states as pointed out in eqs. (3.9) and (3.13). In order to reduce contributions from excited states the choice of the starting point of the fit range is important. As an example, figure 1 demonstrates that, with increasing source-sink distance, the excited states fall below the noise and plateaus of the correlation functions for R ± emerge. The starting time t start is then chosen in such a way that fits with even larger starting times no longer show any systematic trend in the fitted values. Multi-state fits (over larger fit ranges) yield consistent results.

Renormalization
Having computed the bare values of the second DA moments, we are left with the task of renormalizing these bare quantities to obtain results in the standard continuum MS scheme, which is used in the perturbative calculations of the exclusive reactions discussed in the introduction. As already mentioned above, our bare operators are chosen such that there is only mixing between the respective + and − operator multiplets, so we have to determine 2 × 2 mixing matrices such that One then obtains for the second moments of the DAs in the MS scheme a ,

2,MS
with the renormalization factors Z V and Z T of the vector and the tensor currents, respectively. Note that one cannot expect ζ , 22 to be equal to one, since the Leibniz rule holds on the lattice only up to discretization artefacts.

Fit
Fit interval n loops Lattice Table 3. Choices for the fits of the renormalization and mixing coefficients.
We want to evaluate the renormalization and mixing coefficients non-perturbatively on the lattice employing a variant of the RI ′ -MOM scheme, because lattice perturbation theory is not sufficiently reliable. Since forward matrix elements of the + operators vanish in the continuum limit, we cannot work with the momentum geometry of the original RI ′ -MOM scheme but must use a kind of RI ′ -SMOM scheme [25]. We follow exactly the same renormalization procedure as in our investigation of the pion DA [14]. Thus, we need the MS vertex functions of our operators in order to convert the results from our SMOM scheme to the MS scheme. While these are known to two loops in the longitudinal case, i.e., for the operators (2.22a), see ref. [26], as well as for the currents (2.21), see ref. [27], the corresponding results for operators with derivatives involving the matrix σ µν , e.g., the operators (2.22b), are not yet available in the literature. Therefore we discuss the latter case, the so-called transversity operators, in appendix A.
In the end, we determine the matrix Z(a, µ 0 ) (and analogously ζ(a, µ 0 )) at the reference scale µ 0 = 2 GeV by fitting the chirally extrapolated Monte Carlo results Z(a, µ) MC with the expression where the three matrices A i parametrize the lattice artefacts and W (µ, µ 0 ) describes the running of Z in the three loop approximation of continuum perturbation theory.
Ignoring the very small statistical errors, we estimate the much more important systematic uncertainties of Z(a, µ 0 ) by performing a number of fits, where exactly one element of the analysis is varied at a time. More precisely, we choose as representative examples for fit intervals 4 GeV 2 < µ 2 < 100 GeV 2 and 2 GeV 2 < µ 2 < 30 GeV 2 , and we use the expressions for the conversion functions with n loops = 1, 2. For the parametrization of the lattice artefacts we either take the complete expression in eq. (5.6) or we set A 3 = 0. Finally, we consider values for r 0 and r 0 Λ MS corresponding to the results given in ref. [28]. The various fit possibilities are compiled in table 3.
As in the case of the pion DA, the largest effect comes from the variation of n loops . In order to obtain our final numbers for the second moments of the DAs we extract them from the bare data R , ± using each of these sets of values for ζ 11 , ζ 12 and ζ 22 . So we have six results for each of our gauge field ensembles. As our central values we take the results from Fit 1. Defining δ i as the difference between the result obtained with the ζs from Fit i and the result determined with the ζs from Fit 1, we estimate the systematic uncertainties due to the renormalization factors as δ 2 2 + (0.5 · δ 3 ) 2 + δ 2 4 + δ 2 5 + δ 2 6 . Here we have multiplied δ 3 by 1 2, because going from two loops to three or more loops in the perturbative conversion functions is expected to lead to a smaller change than going from one loop to two loops. This should amount to a rather conservative error estimate. The renormalization factors Z V and Z T needed for the evaluation of f ρ and f T ρ , respectively, are calculated in the same way.

Data analysis
From the bare values of f ρ etc. we obtain renormalized results in the MS scheme with the help of our renormalization (and mixing) coefficients on each of our gauge field ensembles. With the range of ensembles available (see table 2) we are able to study the pion mass dependence and, to only a limited extent, volume and discretization effects. Since our lattice spacings do not vary that much, a continuum extrapolation cannot be attempted. Moreover, the impact of the finite lattice size on matrix elements of possibly unstable states is not straightforward. So we take into account results from all lattice spacings and volumes for our final numbers.
Considering the pion mass dependence, we make use of Chiral Perturbation Theory (ChPT) for vector mesons [29][30][31] to obtain the one loop extrapolation formulae for the decay constants Details on the ChPT calculation are given in appendix B. For 2m π < m ρ , i.e., below the decay threshold this infinite-volume calculation yields complex numbers. However, as we neglect instability effects in our lattice computation, which is necessarily done on finite volumes, we use only the real part to fit the mass dependence of our data. The pion decay constant F π is kept fixed at its physical value 92.4 MeV, and the chiral renormalization scale µ χ is chosen to be 775 MeV.
Estimates within ChPT suggest that the third-order term ∝ m 3 π is not negligible for most of our masses (see appendix B). Our data confirm this expectation -the third order term is required in order to fit over the full range of pion masses. Consistent fits are obtained including only second order terms for m π < 300 MeV, however, we have, essentially, only two pion masses in this range. Alternatively, one can ignore the information from ChPT and perform polynomial fits, i.e., drop the logarithmic term in the fit functions (6.1). This  yields very similar results. We expect that a fit including the larger pion masses will yield more reliable numbers than simply taking the values at m π = 150 MeV as our final results because, in particular, the lattice used at this pion mass is relatively small.
In order to get at least some idea of the influence of the instability of the ρ we perform two kinds of fits, including all masses or excluding the results at m π = 150 MeV, which should suffer most from the decay. The resulting ChPT extrapolations for the normalization constants are shown in figure 2. Note that the extrapolated values at the physical point are reasonably consistent with the data at the lowest pion mass.
For the second moments of vector meson distribution amplitudes (see figure 3) no ChPT  calculations are available. It is known that these quantities for pseudoscalar mesons [32,33] and the nucleon [34] do not contain chiral logarithms in leading one-loop order. The reasons are rather generic and may apply to vector mesons as well. Therefore we stick to simple linear fits in m 2 π depicted in figure 3. There is no discernible dependence on the lattice spacing. Errors stemming from the renormalization constants are not included in figures 2-3. We perform an extrapolation for every choice given in table 3 and compute the error of the extrapolated result at the physical point caused by the uncertainties of the renormalization factors from the differences of the extrapolated numbers as indicated at the end of section 5.
Although our data do not allow us to study finite-size and discretization effects systematically, we can make some observations. Considering volume effects, for β = 5.29, κ = 0.13632 we have ensembles with three different volumes at our disposal (m π L = 3.4-6.7). The effects for the decay constants are sizable, see figures 2 and 4. Unlike the well-known cases of pseudoscalar meson and baryon masses, the chiral extrapolations cannot be converted directly to predictions for the leading large-volume behavior. The problematic contributions cancel, however, in the ratio of the decay constants f T ρ f ρ , so that it is straightforward to compute the leading finite-volume corrections for this ratio (see appendix B). It turns out that the corrections are numerically tiny so that from the ChPT  analysis one expects that finite-volume effects for f T ρ f ρ are much smaller than for the couplings themselves. This expectation is in agreement with our data, as shown in figure 4: The finite-volume effects for the ratio f T ρ f ρ (right panel) are considerably smaller than for the vector coupling f ρ itself (left panel). Since in phenomenological studies of hard reactions f ρ will always be set to the experimental value, the ratio f T ρ f ρ , which is not experimentally accessible, is a much more interesting quantity. So we do not perform an infinite-volume extrapolation for f ρ and use this measurement mostly for normalization purposes (e.g. computing the second moments). On the other hand, the observed volume dependence of f T ρ f ρ is small compared to the statistical errors and an infinite-volume extrapolation would not have any significant effect.
One can see from figure 3 that the second moments tend to increase with the spatial volume, however, less significantly than for the normalization constants and the data points have comparatively much larger error bars. As mentioned above, the ChPT analysis of the second moments is not available but the corresponding quantities for stable hadrons have no leading chiral logarithms and a very mild finite-volume dependence. We have checked that excluding the smallest-volume lattice with m π L = 3.4 from the fits does not have any  Table 4. Results in the MS scheme at µ = 2 GeV from the two analysis methods explained in the main text. The numbers in parentheses denote the statistical error and our estimate of the uncertainty introduced by the renormalization procedure. noticeable influence on our results.
Discretization errors are notoriously difficult to control. A certain insight can be obtained looking at the quantities ⟨1 2 ⟩ MS and ⟨1 2 ⟩ MS , which indicate the violation of the Leibniz rule at finite lattice spacing. In the continuum limit they should equal one for all pion masses. Results for all ensembles are plotted in figure 5. Again only statistical errors are shown, the uncertainties resulting from the renormalization coefficients are much smaller. While ⟨1 2 ⟩ MS equals one within the statistical errors with a maximal deviation of about 1%, we observe deviations from one of up to 2% for ⟨1 2 ⟩ MS . Note that these deviations are noticeably smaller than what we found in the case of the pion [14].

Results and conclusion
In table 4 we compare the results of the two kinds of final fits that we have performed. The values in the row labelled "analysis 1" have been obtained by fits to all data points, while the row labelled "analysis 2" contains the results from fits where the data with the smallest pion mass have been excluded. In the case of f ρ , f T ρ , and f T ρ f ρ we have used the fit functions (6.1), whereas the second Gegenbauer moments have been fitted with linear functions of m 2 π . One sees that the results of the two fits are in very good agreement, which may be an indication that ρ-meson decay, ρ → ππ, is not of major importance for the shortdistance quantities that we are considering here. Discretization errors and finite-size effects might be more important, but, unfortunately, cannot be estimated reliably using the set of lattices at our disposal. We expect to be able to quantify the discretization errors using the new N f = 2 + 1 lattice configurations that are generated currently in the framework of the CLS initiative [15].
Comparing to the pion case we observe that for the ρ meson we are able to access the second Gegenbauer moments using momenta with a single non-zero component (see eqs. (3.8) and (3.12)), while we have to consider momenta with two non-vanishing components in order to compute a 2 in the pion. This helps to reduce the statistical noise and the corresponding error.
leading order of perturbation theory with N f = 2. In view of the fact that the systematics are not yet fully controlled, the discrepancies do not look worrying.
In table 5 we compare our main results, i.e., the values of f T ρ f ρ , a 2 , and a 2 , with QCD sum rule estimates and older lattice data. The statistical and renormalization errors of our results have been added in quadrature. Again, the sum rule numbers at µ = 2 GeV have been obtained from the original results at µ = 1 GeV by leading order evolution with N f = 2.
Some of these quantities have already been investigated on the lattice. The BGR collaboration [10] has evaluated the ratio f T ρ f ρ in the quenched approximation with chirally improved fermions at a lattice spacing a = 0.10 fm and found f T ρ f ρ = 0.742(14) at µ = 2 GeV. Further related results have been reported in refs. [8,9,11]. The RBC and UKQCD collaborations [12] have used N f = 2 + 1 domain-wall fermions at a lattice spacing a = 0.114 fm and masses down to m π = 330 MeV to obtain f T ρ f ρ = 0.687 (27) at µ = 2 GeV. In ref. [13] they found ⟨ξ 2 ⟩ = 0.27 (1)(2) at the same scale. Adding the two errors in quadrature and utilizing the relation (2.15) yields a 2 = 0.20 (6).
All existing results are, generally, in good agreement, apart from the ratio of decay constants f T ρ f ρ , which in our case is somewhat smaller than the values obtained in other investigations. This ratio depends strongly on the pion mass, cf. figure 2, and the extrapolation could be affected by the ρ → ππ decay at this level of accuracy. Clarification of this issue by doing a Lüscher-type analysis including four-quark interpolators would be highly desirable since the tensor coupling enters the QCD calculations of the B-decay form factors at large recoil (see, e.g., ref. [17]), where, in some cases, there is a tension with predictions of the Standard Model. Our value for the second Gegenbauer coefficient a 2 is significantly more precise compared to previous results. At this level of accuracy, we start to be sensitive to the difference between the longitudinally and transversely polarized mesons. Our results suggest that a 2 may be slightly larger than a 2 , although the difference is not yet statistically significant. The 20% accuracy for a 2 achieved in our work is interesting for studies of deeply-virtual vector meson production in electron nucleon scattering using the GPD formalism [5]. Such processes will be investigated with high priority at the JLAB 12 GeV upgrade and, in the future, at the EIC.
The work reported here will be continued using CLS N f = 2 + 1 lattice configura-tions [15]. Apart from the study of discretization errors our goal is to consider DAs of the whole SU (3) f meson octet, with emphasis on properties of the K * meson, which is of prime importance for flavor physics. This work is in progress.

Acknowledgments
This work has been supported by the Deutsche Forschungsgemeinschaft (SFB/TRR 55), the Studienstiftung des deutschen Volkes, and the European Union under the Grant Agreement IRG 256594. The ensembles were generated primarily on the QPACE systems of the SFB/TRR 55. The BQCD [35] and CHROMA [36] software packages were used, along with the locally deflated domain decomposition solver implementation of openQCD [37,38]. We thank Benjamin Gläßle for software support. Part of the analysis was also performed on Regensburg's Athene HPC cluster, the Regensburg HPC-cluster iDataCool and computers of various institutions which we acknowledge below. One

Appendices A Transversity operators in the continuum
In this section we review our construction of the continuum Green's functions which will be used for connecting the MS scheme to the RI ′ -SMOM scheme employed on the lattice. The procedure we follow has already been applied to several similar quark bilinear operators [26,27,40] and we will highlight the salient differences for the transversity operators considered here. The notation of this section very much runs parallel to, for instance, ref. [40], to which we refer the interested reader for more background. First, the two classes of operators we are interested in are the flavor non-singlet operators, where the operators with a single derivative have been included for completeness. We define ς µν = 1 2 [γ µ , γ ν ] which is related to σ µν by Our use of ς µν is to retain the same conventions with earlier renormalization of similar operators [26,27,40] and our use of generalized γ-matrices which we discuss later. To define the action of the symbol S, which imposes certain symmetrization and tracelessness conditions, it is best to consider the generalized transversity operators O T µν 1 ...ν i ...νn from which we will focus on the values of n = 2 and 3. Specifically, [41], where the label T includes all possible total derivative operators. When n = 2, for example, then for the first operator of the T 2 sector with again a parallel definition for the total derivative operator [42]. In our construction for the T 3 operators we have taken the convention to include an extra factor of 1 6 in the definition of S. We will use T 2 and T 3 to refer to a sector as well as for the non-total derivative operator of each set. It will be clear from the context which is meant. The labelling for each derivative of a total derivative operator is one ∂ symbol applied to the sector label. In defining the operators we have omitted the explicit flavor indices and note that our perturbative renormalization will be for massless quarks; in other words we are in the chiral limit. The total derivative operators are required since there is operator mixing within each separate sector. It would not usually be necessary to include these but since the Green's functions they are needed for are non-forward matrix elements, then a momentum will flow through the operator insertion and the mixing will be activated. Part of the evaluation of these matrix elements requires the renormalization of the operators. Again our basis choice is partly driven by the need to simplify this aspect. Operators with the same quantum numbers and dimension will mix under renormalization. However, for our choice the mixing matrix will be upper triangular. For instance, we have for l = 2 and 3 where the subscript o denotes the bare operator. Then We use 1 and 2 to label the elements of the T 2 matrix where 1 is the operator T 2 . Similarly 1, 2 and 3 label the T 3 matrix elements which respectively correspond to T 3 , ∂T 3 and ∂∂T 3 . The explicit mixing matrix for the T 2 system has been determined in ref. [42] to three loops in the MS scheme. Prior to the results we present here, the T 3 matrix was known only partially to the same order. Entry (ij) = (11) is the renormalization constant for the operator T 3 itself and the remaining two diagonal entries are the same as the operator T 2 and the tensor current [42][43][44]. In other words the operators of the T 2 system without the total derivatives. In addition the off-diagonal element (ij) = (23) is known purely because the non-zero entries of the final two rows of Z T 3 ij are the non-zero entries of the Z T 2 ij matrix. We have determined the final two off-diagonal elements of Z T 3 ij by renormalizing the operators in a quark 2-point function where the momentum of one of the external quark legs is nullified. In other words there is a non-zero momentum flowing through the inserted operator. This was the method used to determine a similar mixing matrix for the third moment of the usual twist-2 Wilson operators in deep inelastic scattering [26]. However, in ref. [26] it was noted that such a computational setup was not sufficient to determine each of the (ij) = (12) and (ij) = (13) elements separately. To disentangle them an extra piece of information was required. This is achieved here for T 3 by the identity which is straightforward to establish by integration by parts. Thus to deduce these remaining two off-diagonal elements we have applied the Mincer algorithm [45] to the three loop renormalization of the operator T 3 . As the resulting anomalous dimensions for T 2 are given in ref. [42], we record the first row of the three loop anomalous dimension mixing matrix for T 3 which is as the remaining rows are given in ref. [42] where a = g 2 (16π 2 ). Here ζ n is the Riemann zeta function. We note that our anomalous dimensions pass all the usual consistency checks. In particular we derived (A.8) in an arbitrary linear covariant gauge and checked that the gauge parameter cancels as it ought to for gauge invariant operators in the MS scheme.
Having summarized the renormalization of the operators of interest the next stage is to provide the perturbative corrections to the Green's function where the operator is inserted in a quark 2-point function. As we are considering non-forward matrix elements there is a momentum flowing through the operator. More specifically we consider the Green's function ⟨ψ(p)O T l µ 1 ...µ l+1 (−p − q)ψ(q)⟩ for the two cases l = 2 and 3. There are two independent external momenta p and q and we will evaluate the Green's function at the fully symmetric point given by where µ is a mass scale. For this section we will take this scale to be the same mass scale that is used in dimensional regularization in d = 4 − 2 dimensions to ensure the coupling constant is dimensionless in d-dimensions. Therefore, our results for the Green's function will not have any logarithms of mass parameter ratios. As each Green's function has free Lorentz indices we choose to decompose them into a basis of Lorentz tensors denoted by P T 2 (k) µνσ (p, q) and P T 3 (k) µνσρ (p, q). Here T 2 and T 3 indicate the sector as the basis will be the same for the Green's function with the total derivative operators of each sector too. The choice of tensors in each basis is not unique. However, each basis is large due to the number of objects available to build the tensors. These include the momenta p µ and q µ as well as η µν . In addition there are Lorentz tensors built from the γ-matrices. As in previous perturbative evaluations [26,40] we use the generalized γ-matrices of [46][47][48]  for integers n ≥ 0. In the definition an overall factor of 1 n! is understood. These matrices span the spinor space when dimensional regularization is used. As an aside we note that it is in this context that our choice of ς µν in the operator definition fits naturally. The algebra and properties of these matrices is well-established [49,50]. We note one specific property which is important here which is where there is no sum over repeated m or n and I µ 1 ...µmν 1 ...νn is the generalized unit matrix.
The key point is that this trace partitions the space spanned by the tensors in the basis into distinct sectors. As we consider the operators in massless QCD, only Γ (0) , Γ µν (2) and Γ µνσρ (4) will be needed. For T 3 it might be expected that Γ µ 1 ...µ 6 (6) would be required but the symmetrization conditions exclude this γ-matrix from the basis. Finally with these objects we have constructed the tensor basis for each sector. For T 2 that involves 30 tensors consistent with the symmetry properties of the inserted operator. A sample set is presented below. For T 3 there are 42 tensors and for space reasons these as well as the full T 2 set are given in the attached data file. The next step is to compute the coefficients in the decomposition of each Green's function into their respective tensor basis. In other words we need the values of the amplitudes (A.13) with N 2 = 30 and N 3 = 42. The factor C l is a sector specific normalization to account for the differing dimensionalities of the tensor basis and Green's functions for each sector. Thus we have C 2 = −i and C 3 = µ 2 . The algorithm to determine these coefficients has been given in refs. [26,40] for instance. Briefly, to apply the multiloop perturbative integration techniques to find these amplitudes we have to extract scalar Feynman integrals which is achieved by a projection method. The projection matrix, M T l ij , required for each sector is constructed from the respective tensor basis [26,40] as it is the inverse of the matrix (A.14) Due to the size of the matrices, their explicit form is given in the auxiliary data file provided. Nevertheless, the partitioning due to the generalized γ-matrices provides a computational shortcut. Hence we have (A.15) Next we briefly note the practical details of actually carrying out the two loop evaluation of the Green's function which proceeds in an automatic way. The Feynman diagrams are generated using the Qgraf package [51]. These have to be converted to Form [52,53] notation after all the Lorentz and color indices have been included. There are 3 graphs at one loop. At two loops there are 32 graphs for O T 2 µνσ and 37 for O T 3 µνσρ with fewer graphs for total derivative operators. After this the Feynman rules for either operator together with the propagators and vertices are substituted and the various amplitudes are projected out to produce a large number of scalar Feynman integrals that need to be calculated. To achieve this we have used the Laporta algorithm approach [54]. After projection the scalar products of the momenta in the numerators of the integrals are written in terms of the propagators. In addition there may be propagator forms which are not present which are referred to as irreducible. In this format the Laporta algorithm [54] is then applied which uses integration by parts to systematically construct all the algebraic relations between reducible and irreducible scalar integrals for a specific momentum topology. The upshot is that all the required scalar integrals are written in terms of a small basis of master integrals whose expansion is known from direct computation [55][56][57][58]. Therefore, we are able to reduce all the scalar amplitudes to known integrals and hence evaluate them exactly at one and two loops. Whilst this is in essence the Laporta method [54] one has to construct the relations in a practical way. We have chosen to use the Reduze package [59]. Moreover, the output files from the database that Reduze builds is straightforward to interface with the Form modules that constitute the automatic computation. For the two loop calculation we perform here, it transpires that for the Reduze setup there is only one momentum topology at one loop and two at two loops. The latter are the ladder and non-planar topologies. All the Feynman diagrams that we have to compute can be mapped into these three cases. The final stage is to carry out the overall renormalization. This is achieved by computing all the graphs as a function of the bare parameters, such as the coupling constant and gauge parameter, following the procedure introduced in ref. [60] for automatic symbolic manipulation loop calculations. Then the renormalized parameters are introduced via the usual renormalization constant definitions with the operator renormalization constants being extracted at the end to leave the finite expressions for each scalar amplitude.
To allow orientation to the full data available in the attached data file we give a selection of the various amplitudes. We provide these in numerical form for one representative from each Γ (n) -matrix partition for both operators of the T 2 sector. For instance, we have where α is the gauge parameter and the restriction ⋯ stands for evaluation at (A.9) and (A.10). Although we are only interested in the values in the Landau gauge, defined by α = 0, we have performed our computations for arbitrary α. This is mainly as a check on the renormalization of the operators since their anomalous dimensions are independent of α in the MS scheme.
Next we summarize some aspects of the tensor basis and projection matrix for the T 2 sector. Indeed one purpose of this summary is to provide an aid to the understanding of the full information given in the attached data file for both T 2 and T 3 . Due to the size of the bases and matrices we used, a useable electronic format is more appropriate for their representation. First, we present a selection of tensors in the T 2 basis choosing several representatives from each Γ (n) -matrix partition. When one of the external momenta is contracted with a Lorentz index then that momentum appears as an index. For example, for T 2 we have (2)µνσ (p, q) = ς µν q σ + ς µσ q ν + [2ς µq q ν q σ − ς νq q µ q σ − ς σq q µ q ν ] 1 µ 2 , P T 2 (5)µνσ (p, q) = ς νp η µσ + ς σp η µν + [2ς µp q ν q σ + dς νp q µ q σ + ς νp q µ q σ + dς νq p µ q σ + 3ς νq p µ q σ + dς νq p σ q µ + ς νq p σ q µ + dς νq q µ q σ + 2ς νq q µ q σ + dς σp q µ q ν + ς σp q µ q ν + dς σq p µ q ν + 3ς σq p µ q ν + dς σq p ν q µ + ς σq p ν q µ + dς σq q µ q ν + 2ς σq q µ q ν ] 1 µ 2 , We have only shown one tensor from the final partition as the other is given by replacing the uncontracted vector p by q.
For each of the bases we have explicitly constructed the projection matrix coefficients. For T 2 as there are 30 projectors this would correspond to a 30 × 30 matrix where the entries are rational polynomials in d. However, as we are using the generalized basis of γ-matrices in d-dimensions the projector matrix is block diagonal due to the property of (A.12). In other words where the subscript on the block matrices corresponds to the label of the analogous Γ (n)matrix appearing in the projection tensor. Each of these partitions is of different size being respectively 22, 6 and 2 dimensional. Given the size of the first two submatrices it is again not feasible to display all entries. Instead we choose to give a few reference entries to facilitate the extraction of the full matrices from the data file. We have where indices of M T 2 (0) i j range from 1 to 6 and these can be mapped to the labels of the tensor basis by adding 22. Finally, the remaining sector is compact enough to record it completely as Overall the matrix M T 2 is symmetric as is M T 3 . Finally, this information should be sufficient to connect with the full electronic representation for both sectors.

B.1 Effective field theory framework
In the specific framework of Chiral Perturbation Theory (ChPT, see, e.g., refs. [61][62][63]) applied here, the generating functional of all QCD correlators is evaluated by means of a path integral involving an effective low-energy Lagrangian L eff (U, v, a, s, p, . . .) (compare with ref. [61], and eqs. (1) and (2) of ref. [30]), Formally, all QCD Green's functions can be obtained by taking functional derivatives of the generating functional w.r.t. the external (Hermitian) scalar, pseudoscalar, vector, axialvector and antisymmetric tensor source fields s, p, v µ , a µ ,t µν . It should be noted that the tensor structure with an additional γ 5 is not independent due to the identity σ µν γ 5 = i 2 µνρσ σ ρσ . The dots stand for other possible source fields (for example, the coupling to symmetric tensor fields has been considered in ref. [64,65]). The tensor sourcet µν has been incorporated in ref. [66]. The matrix field U collects the pion (Goldstone boson) fields in a convenient way (see below). The effective Lagrangian has to be invariant under local chiral transformations of the Goldstone boson and source fields, and shares all other symmetries of L QCD . A formal proof that low-energy QCD can indeed be analyzed in this way has been given by Leutwyler [62]. Under chiral transformations (L, R) ∈ SU (2) L × SU (2) R , the quark and external source fields transform as The effective Lagrangian and the perturbative series are ordered by a low-energy power counting scheme, counting suppression powers of Goldstone boson momenta and masses (or quark masses). For details and further references, we refer to refs. [61][62][63]. At leading chiral order, the effective Lagrangian describing the interaction of the pseudo-Goldstone bosons (pions) with the external source fields and each other is given by (see ref. [61]) with χ = 2B(s + ip), s = M + δs, where M is the quark mass matrix, and δs the remaining part of s. The brackets ⟨⋯⟩ denote the flavor (or isospin) trace, F is the pion decay constant in the chiral limit (F ≈ 86 MeV), and where j is a channel (particle species) index which labels the specific pion, andλ are the pertaining channel matrices. We write out φ as where the τ a are the Pauli matrices. The matrix field U transforms as U → RU L † under chiral transformations. We also introduce u as the square root of U , u 2 = U , which transforms as u → RuK † = KuL † , thus defining the so-called compensator matrix K = K(L, R, U ) (which is also unitary). Below we shall set the external fields p, a to zero, s = M (the quark mass matrix) and set v µ = v a where we only show the terms needed for our present work (see refs. [61,66,67] for the complete Lagrangian at that order, and eq. (B.7) for the definition of the operators u µ , F ± µν and T ± µν ).

B.2 Chiral Lagrangians for resonances
Explicit vector meson degrees of freedom have been incorporated in the effective Lagrangian of ChPT already in ref. [29,30]. In the following, a "heavy vector meson" framework was set up [31,[68][69][70][71] to deal with problems related to the modified power counting in the extended effective theory, caused by introducing a new heavy mass scale (the vector meson mass in the chiral limit). Today, it is better understood how to deal with such problems in a manifestly Lorentz-invariant way, by employing the freedom of choice of the subtraction scheme for the effective field theory [72][73][74]. Such methods have been applied to the case of heavy meson resonances in refs. [75][76][77][78][79][80][81][82]. We refer to these references for details on the vector meson effective field theory outlined below. Keeping in mind the transformation behavior of the external source fields v µ , a µ and t µν given above, we can write down the following terms describing the interaction of the vector mesons with the external source fields and the pions (compare also the previous references, and ref. [83]): see eq. (B.3) for the channel matricesλ. We have used a large-N c argument here to cast the ρ and ω fields in the matrix form of the last line in eq. (B.7), compare also with eq. (27) of ref. [31]. The dots indicate terms of higher chiral order, terms involving external source fields s, p (which are not needed here), or terms involving more derivatives, which result in contributions of the same form as those resulting from the terms given above, when using the equations of motion or field transformations [84]. The vector field propagator in momentum space is
It easily follows that the form factors f v,t ππ must have the phase δ 1 1 (s) in the elastic region.

B.4 Contributions to ρ matrix elements
Here we use the definitions ⟨0 q τ a 2 γ µ q ρ b (k, λ)⟩ = δ ab m ρ e (λ) µ f ρ √ 2 , (B.14) ⟨0 q τ a 2 σ µν q ρ b (k, λ)⟩ = iδ ab (e (λ) µ k ν − e (λ) 15) and find at the one loop level up to O(p 4 ) Here we have to set k 2 ≡ s equal to the rho pole, k 2 → s pole = m 2 ρ − im ρ Γ ρ [85]. The wave function renormalization factor is derived from the ρ self-energy Π ρ (s), where the contribution of the one loop graphs to the self-energy is given by (compare ref. [80]) of order ∼ 3 GeV −3 . Inserting this estimate, and µ χ = 770 MeV, the third-order term becomes comparable to the leading chiral logarithm for M ≳ 200 MeV, so it might give a non-negligible contribution for most data points. In eq. (B.21), we have written the result for the chiral expansion of Re f T ρ f ρ , which motivates the extrapolation formula (6.1c), while the formulae (6.1a) and (6.1b) result from (B.16) and (B.17), respectively, upon inserting the explicit expressions for the loop functions given below. The cusp effects and imaginary parts in the chiral behavior of the couplings could only be extracted indirectly from the computed correlators, which are real on Euclidean lattices with a finite volume. A more thorough analysis is needed to deal with such complications. It is, however, important to note that the leading non-analytic terms given above are not afflicted by this deficiency. This can be deduced from the fact that they agree with the corresponding results in the heavy vector meson framework [67,71], where the unitarity effects due to the ππ loops are either dropped or derived from contact terms of a non-Hermitean part of the effective Lagrangian (see, e.g., ref. [68]).
In the expression for the ratio given in eq. (B.21), the factors of Z χ ρ and the nonanalytic terms in the loop function I A ππ containing the imaginary part cancel at one-loop order. Due to this simplification, it is straightforward to compute the finite-volume corrections for this ratio. Here, we attempt only an estimate of the leading finite-volume correction, related to the O(M 2 ) 'chiral-log' term contained in the tadpole loop integral I π (compare eq. (B.27) below). According to the standard formalism of ChPT in a finite cubic volume V = L 3 [86], this loop integral is replaced by its finite-volume counterpart Here K 1 (z) is the modified Bessel function of the second kind, which decays exponentially for large positive z, K 1 (z) →

B.5 Loop functions
To render this appendix self-contained, we give the definitions of the loop integrals occuring in the formulae above. The loop integral with two pion propagators is given by , k 2 =∶ s . Note that this integral has an imaginary part ImĪ ππ (s) = −σ 0 (s)θ(s − 4M 2 ) (16π) for real s > 4M 2 . We have also employed the abbreviation and we refer to appendix B of ref. [80] for details on the chiral expansion. We also use where I V is given by the formula for I π with M → m V . Here, the letter V stands for the vector meson running in the loop (ρ, ω, . . .).