$\eta$-$\gamma$ and $\eta'$-$\gamma$ transition form factors in a nonlocal NJL model

We study the $\eta$ and $\eta'$ distribution amplitudes (DAs) in the context of a nonlocal SU(3)_L x SU(3)_R chiral quark model. The corresponding Lagrangian allows to reproduce the phenomenological values of pseudoscalar meson masses and decay constants, as well as the momentum dependence of the quark propagator arising from lattice calculations. It is found that the obtained DAs have two symmetric maxima, which arise from new contributions generated by the nonlocal character of the interactions. These DAs are then applied to the calculation of the $\eta$-$\gamma$ and $\eta'$-$\gamma$ transition form factors. Implications of our results regarding higher twist corrections and/or contributions to the transition form factors originated by gluon-gluon components in the $\eta$ and $\eta'$ mesons are discussed.

and there is no natural way to normalize the gluon DA. The prefactor present in Eq. (7) is just a convention, and other definitions can be found in literature (see the discussion in Ref. [20]). A change in this prefactor can be compensated through a factor into the integrand in the second term of Eq. (1).
In the case of the π → γγ * TFF the situation is simpler, since there is no singlet contribution.
One has where the pion DA is given by B. Neutral pseudoscalar meson distribution amplitudes in a nonlocal NJL model We consider here the meson DAs within a nonlocal NJL (nlNJL) model. The corresponding Euclidean effective action, in the case of SU(3) F flavor symmetry, is given by [30] where ψ(x) is the SU(3) F fermion triplet ψ = (u d s) T , andm = diag(m u , m d , m s ) is the current quark mass matrix. We will work in the isospin symmetry limit, assuming m u = m d . The model includes flavor mixing through the 't Hooft-like term driven by H, where the constants A abc are defined by with a = 0, . . . , 8. The fermion currents are given by where the functions G(z) and F (z) are covariant form factors responsible for the nonlocal character of the interactions. Notice that the relative weight of the interaction driven by j r (x), which leads to quark wave function renormalization, is controlled by the parameter κ. In the mean field approximation (MFA), which will be used here in what follows, the quark propagator for each flavor f = u, d, s can be expressed as where M f (p) and Z(p) stand for momentum dependent effective mass and wave function renormalization (WFR), respectively. One has [30] M f (p) = Z(p) m f +σ f g(p) , where g(p) and f (p) are Fourier transforms of G(z) and F (z), whileσ f andζ are mean field values of scalar fields associated with the corresponding currents in Eq. (12). Details of the procedure carried out to obtain these quantities are given in App. A.
The momentum dependence of the interaction form factors can be now obtained from lattice QCD results. Following the analysis in Ref. [32], the effective mass M u (p) can be written as where f m (p) = 1/ 1 + ( with α = 3/2. From Eqs. (14) one has then α m = (m uζ /κ +σ u )/(1 −ζ/κ). For the wave function renormalization we use the parametrization [33,35] where f z (p) = 1/ 1 + p 2 /Λ 2 Here the new parameter α z is given by α z = −ζ/(κ −ζ). The functions f (p) and g(p) can be now easily obtained from Eqs. (14)(15)(16)(17)(18). As shown in Refs. [33,35], for an adequate choice of parameters these functional forms can reproduce very well the momentum dependence of quark mass and WFR obtained in lattice calculations. We complete the model parameter fixing by taking as phenomenological inputs the values the of the pion, kaon and η ′ masses and the pion weak decay constant [30].
The resulting model parameters are given in Table I.
In the case of the above described nonlocal model, however, the situation is more complicated since the inclusion of gauge interactions implies a change not only in the kinetic piece of the Lagrangian but also in the nonlocal currents appearing in the interaction terms. If x and z denote the space variables in the definitions of the nonlocal currents [see Eq. (12)], one has Here the function W (s, t) is defined by where r runs over an arbitrary path connecting s with t.
This procedure has been analyzed in detail within nlNJL models, in particular regarding the calculation of the pseudoscalar meson decay constants [29,33,38], see Eq. (5). The situation is similar for the case of the bilocal axial current in the definition of the meson DA. In fact, the basic physical idea beyond the factorization of the meson TFF into hard and soft contributions is that for high Q 2 the struck quark loses its high momentum before being able to interact with the remaining quarks and gluons of the hadron (Q 2 ∼ 1 GeV 2 implies a time scale of the order of 10 −24 s).
Therefore, the nonlocal interaction does not see the struck quark but only the quarks in the hadron before and after the photon absorption-emission process. This can be effectively implemented by introducing an external fictitious probe carrying the adequate quantum numbers, which in our case is a gauge axial field (a similar situation has been studied in the case of the pion parton distribution, see Refs. [33,34]). Thus, the axial vertex in Eq. (4) will become dressed by the nonlocal interaction, irrespective of whether the quark current is a local or a bilocal one (as in this case).
The steps to be followed in the explicit calculation of the quark DA within the nlNJL model are detailed in Appendix A. We quote here the resulting expression. In the flavor basis (i.e. q i = q ℓ , q s ) where g M qq stands for an effective quark-meson coupling constant [see Eq. (59) in Appendix A] and the integration variables are related to the meson and quark Euclidean four-momentum P and k, respectively. Considering the light front variables in the frame where the transverse component P T vanishes, the invariants k 2 and k · P can be written in terms of the variables w and k T as It is convenient to separate the integrand in Eq. (22) into two pieces, The explicit expressions for these functions are where Here α − g and α + f depend in general on the integration path in Eq. (21). For a straight line path one has where g ′ (k) ≡ dg(k)/dk 2 , and same for f ′ (k).
It is important to mention that, even when our effective model leads to an adequate phenomenological pattern for low energy meson phenomenology, there are some differences between model predictions and phenomenological values of the η and η ′ decay constants (see Table XII in App. A).
In our numerical calculations, when evaluating the η and η ′ DAs from Eq. (22) we will take the values of f i M arising from our model, in order to guarantee the proper normalization condition Eq. (6). On the other hand, we will use the phenomenological values for f ℓ M or f i M when evaluating the flavor mixing leading to the quark DAs, Eq. (3).
Once the latter are known at a given µ 0 scale, their evolution up to a higher scale µ can be obtained from perturbative QCD. In order to study this evolution it is convenient to expand the DAs in series of Gegenbauer polynomials: where i = 0, 8, and we have now explicitly denoted the µ dependence of the DAs. Notice that only n-even terms contribute to the sums, due to the symmetric (antisymmetric) behavior of the quark M (x, µ) (i = 0, 8) satisfy the sum rule Eq. (6), the first coefficients a M n (µ 0 ) = 144 (2n + 5) (n + 1) (n + 2) (n + 3)(n + 4) Notice that Eq. (30) holds either if one is working in the flavor basis (i = ℓ, s) or in the SU(3) F basis (i = 0, 8). At the LO the Gegenbauer polynomials are eingenfunctions of the evolution kernel, therefore a M n coefficients of different order n do not mix with each other [21]. On the other hand, as stated, QCD evolution equations mix the gluon and singlet quark components for n ≥ 2. The evolution of these coefficients up to a scale µ is given by (see Refs. [20,21]) Here β 0 = 11 − 2 n f /3, n f being the number of active flavors at the scale of the process (in our case we take n f = 4), and γ (±) n are the eigenvalues of the anomalous dimension matrix γ n . These are given by where the (LO) matrix elements of γ n read γ gq n = 36 (n + 1) (n + 2) , The coefficients ρ (+) n and ρ (−) n , which weight the presence of quarks in the gluon DA and gluons in the singlet quark DA, respectively, are given by Finally, the evolution of the strong coupling constant α s at the LO is given by with Λ = 0.224 GeV.
In Table II we quote the values of the anomalous dimensions for the first values of n. Already for n = 2 it is seen that γ (+) and γ (−) are close to γ qq and γ gg , respectively, and the differences become even smaller for larger n. The numerical values for ρ (±) and the product ρ for the first values of n are given in Table III.  n , γ qq n and γ gg n coefficients, and asymptotic values.  In the case of the distribution amplitude φ M , at the LO the evolution is just governed by the anomalous dimension γ qq n . One has We will also take into account the effect of NLO corrections to the DAs. In general, at the NLO the evolution equations for different coefficients a (q i ) M n get mixed, and the pattern becomes more complicated. We will consider the NLO evolution for the octet component (see discussion in the next section). The corresponding coefficients evolve according to [39] Explicit expressions for the renormalization factors E NLO n (µ, µ 0 ), as well as for the off-diagonal mixing coefficients d k n (µ, µ 0 ) in the MS scheme are collected in Appendix B. Usually the calculation of a few coefficients a (q i ) M n (µ) is sufficient to get a good estimate of the πDA at the scale µ from Eq. (29).

A. Quark DAs
From the numerical evaluation of the integrals in Eq. (22) one can obtain the quark DAs for π, η and η ′ mesons within the above described three-flavor nlNJL model. The corresponding curves are displayed in Fig. 1, where for comparison we also include the asymptotic limit φ asym (x) = 6x(1 − x).
The dotted lines correspond in all cases to the asymptotic limit φ asym (x) = 6 x(1 − x).
As stated in the previous section, our calculations have been performed in Euclidean space. The consistency of our procedure has been discussed in Ref. [13], where the pion DA and TFF are analyzed within a two-flavor version of the model. Our main test in this sense is the verification that the sum rule Eq. (6) is satisfied. In the case of η and η ′ mesons, however, the stringency of this test becomes weakened owing to the numerical uncertainties in the calculations. In fact, when computing the integrals in Eq. (22) one has to take into account that the functions F i (w, x, k T ) show cuts in the complex w plane that require a deformation of the integration paths (see e.g. the discussion in App. B of Ref. [40]). In addition, depending on the value of x these functions have poles that also need to be compensated numerically. The normalization of quark DAs obtained from our calculations in the present model are 1.0004 for φ π (x), 0.9989 and 0.9753 for φ (q ℓ ) η (x) and φ (qs) η (x), respectively, and 1.027 and 0.870 for φ (q ℓ ) η ′ (x) and φ (qs) η ′ (x), respectively. It is worth mentioning that the effect of poles and cuts gets increased for higher quark and meson masses, therefore numerical uncertainties are particularly large in the case of the s quark DA in the η ′ meson, where we find the largest departure from the normalization condition (the effect of the error in the determination of the η ′ DAs is discussed above). All quark DAs shown in Fig. 1 have been renormalized so that they satisfy the sum rule.
When using a quark model to describe the deep structure of hadrons it is crucial to establish the chosen scale µ 0 that will be associated with the results provided by the model. In our case the scale should be the same as that used in lattice calculations, since we have taken into account lattice results in order to fix the shape of the form factors in the quark propagators. Thus, from Ref. [31] we take µ 0 = 3 GeV, which is a rather large value in comparison with the scale µ 0 ∼ 0.5 − 1 GeV usually adopted in quark model calculations.
In the left panel of Fig. 1 we show the pion DA. Our result is pretty similar to that obtained within the two-flavor nonlocal NJL model studied in Ref. [13]. Notice that this might not have been the case, since the change from a two-flavor model to the present three-flavor model implies a refitting of all model parameters. By looking at the DAs in Fig. 1 it is seen that in all cases the curves have two symmetric maxima. This is also shown by the π DAs calculated in Refs. [41,42], albeit in our case the two maxima are much closer to x = 0.5. In fact, in the nlNJL model this feature arises from the term F (2) i (w, x, k T ) [see Eqs. (23)(24)(25)], which is a genuine nonlocal contribution. Now, by comparing the curves in the different panels of Fig. 1 one can see the effect of meson and quark masses in the behavior of the DAs. As expected, the π DA at µ 0 = 3 GeV is relatively close to the asymptotic limit φ asym (x) = 6x(1 − x). This also holds for the u (or d) quark DA in the case of the η meson, Finally, in the case of the η ′ meson (right panel in Fig. 1 η ′ (x) lie far from the asymptotic limit. Another important feature, common to all obtained DAs, is that they go to zero rather fast near the points x = 0 and x = 1, supporting the idea of suppression of the kinematic end points [43,44].
Let us consider the QCD evolution of the DAs. We recall that we are working within a quark model in which there is no gluon content. Moreover, according to the numerical values of the ρ (±) n coefficients in Table III    Our results for the case of the η meson can be compared with those obtained within the (local) Nambu−Jona-Lasinio model in Ref. [28], where only the η meson case is analyzed, since the η ′ turns out to be unbounded. It is seen that the shapes of the η DAs are quite different from those obtained in the present model, showing only one central maximum (we recall that the origin of the two-maxima behavior shown in Fig. 1 arises from the purely nonlocal contribution). As expected, the differences in the shapes are translated to the coefficients of the Gegenbauer expansion: the first coefficients

B. TFFs without gluons
In this subsection we present the results obtained within our approach for the pseudoscalar meson-γ TFFs. In fact, we have modified the expression on the right hand side of Eq. (1) by adding subleading terms in an expansion in inverse powers of Q 2 . This procedure has been already used in Refs. [13,15,28,45] in order to account for contributions coming e.g. from higher twist operators. Here we propose to include two terms in this expansion. In addition, let us neglect for now the gluon contribution to the TTFs. This is consistent with a description of mesons within the nlNJL, which has no gluon content. We have in this way In accordance with our approximation of neglecting gluon contributions, we will use octet evolution for the whole quark DAs Φ Finally, the short-dashed curves correspond to what we call the "asymptotic behavior", obtained In Table VII we quote the values of c and d arising from our fits, together with the number of experimental data considered in each case and the corresponding χ 2 values. For comparison we also include the χ 2 obtained when we take c = d = 0. By looking at the χ 2 values it is seen that the introduction of higher twist corrections through the c and d terms leads to a significant improvement in the theoretical description of the data for both the π-γ and η ′ -γ TFFs, while the improvement is not so important in the case of the η-γ TFF. In this regard, notice that the better quality of the fits is basically dominated by the description of the low virtuality region (which has less impact in the case of the η ′ -γ TFF owing to the wide dispersion of the data). In fact, by comparing the solid and dashed curves in the figure we observe that in the case of the π-γ and η-γ TFFs the differences are ruled by the behaviors at Q 2 3 GeV 2 , while for the η ′ -γ TFF there is a more steady deviation which covers a region up to Q 2 ∼ 10 GeV 2 . Moreover, from Table VII it is seen that the signs of c and d are the same for π-γ and η-γ TFFs, whereas they are opposite to those obtained from the fit to η ′ -γ TFF data. This could be related with the octet character of the π and the prevailingly octet character of the η, which contrast with the mostly singlet character of the η ′ .  Table VII. Short-dashed lines show the NLO asymptotic QCD behavior (see text), and dotted lines indicate the QCD asymptotic limits. In the case of the η ′ -γ TFF, the gray region corresponds to a change in a (q i ) η ′ n coefficients within a range of 15%. Notice that in all graphs we have used a logarithmic scale for Q 2 . Now, while higher twist effects influence the low Q 2 region of the TFFs, it is interesting to analyze the high virtuality region from the point of view of QCD, comparing our results with the asymptotic QCD behavior and the asymptotic limit of the TFFs. From the graphs in Fig. 2 it is seen that in all cases the introduction of NLO corrections to the parton level subprocess amplitudes T qq (while keeping the asymptotic limit for the DAs) has a negative contribution to the TFFs. In addition, it is seen that in all cases the results obtained within the nlNJL model approximate experimental data from below.
Let us comment separately the situation for each meson. In the case of the pion, the experimental data seem to cross the asymptotic limit at some Q 2 value between ∼ 10 − 20 GeV 2 , hence the NLO corrections go in the wrong direction. This is a well-known problem that we have already discussed in the context of the two-flavor version of the nlNJL model in Ref. [13]. In fact, the puzzling pion data can be described by some models based on flat DAs and some cutoff in the parton amplitudes [15,45,46].
In the case of the η-γ TFF, even if experimental data for Q 2 > 10 GeV 2 seem to follow the asymptotic behavior, the trend of the data shows that it is not unlikely that the TFF crosses the QCD asymptotic limit for higher Q 2 [28]. In any case, according to present experimental results, it can be said that our model provides a good description of the TFF.
Finally, for the η ′ -γ TFF the experimental data lie clearly below the asymptotic behavior, and quite far from the asymptotic QCD limit. Once again the results obtained within the nlNJL model are shown to be in good agreement with the data. Given the uncertainty in the numerical calculations for the η ′ DA discussed in the previous subsection, we have studied in this case the stability of our results against some variation in the coefficients of the Gegenbauer expansion of the quark DAs. In order to get an estimation of the error we have considered the η ′ -γ TFF for c = d = 0, allowing for a change in a (q i ) η ′ n coefficients (n ≥ 2) within a 15% range. The corresponding range obtained for the TFF is shown by the small gray area in Fig. 2. In general we can state that this error does not affect qualitatively our results.

IV. THE EFFECT OF GLUONS
In this section we discuss the possible effect of the presence of gluon components in the DAs.
According to the discussion in Sec. II.C, it is natural to carry out our analyses using the octet-singlet basis for the quark distribution amplitudes. At any scale µ, we can use Eqs. (28) to obtain the octet and singlet quark DAs from the flavor ones, and analogous expressions can be written for the coefficients of the Gegenbauer expansion. Let us assume that we know the flavor DAs or, equivalently, the coefficients of the Gegenbauer expansion at some scaleμ 0 . For the octet and singlet Gegenbauer coefficients we have At LO, the evolution fromμ 0 up to a higher scale µ is obtained from Eqs. (37) and (32). In particular, for quark singlet and gluon coefficients one has hence In fact, to the order we are working, we should compute the NLO evolution of the DAs. At NLO the evolution of a (q i ) M n (µ) coefficients for different order n become mixed, as one can see from Eq. (38) for the case of octet components. However, the impact of NLO corrections to the evolution of these coefficients is not significant in comparison with the corresponding corrections for the subprocess amplitudes T qq and T gg given in Eqs. (2). Indeed, the most important effect on the DAs when going from LO to NLO evolution comes from the change in the strong coupling constant, α s (µ).
Thus we adopt the following prescription: we consider the NLO corrections for T qq (x, Q 2 , µ 2 ) and T gg (x, Q 2 , µ 2 ) [given by Eqs. (2)] together with Eqs. (37) and (44)(45) for the evolution of the octet and singlet DAs, respectively. In all these equations we take the NLO running equations for the strong coupling constant α s , given by Eq. (64) in App. B. In order to test the validity of this prescription, let us study the case of the octet DAs. In Tables VIII, IX     the NLO evolution equation (64) for α s . From the values in the tables one can conclude that the "mixed" approximation can be used to estimate NLO calculations with reasonably good accuracy.
We consider here two different ways of estimating the effect of gluons in the DAs. Our first analysis is based on the fact that in general one assumes that the scale at which standard quark modelswith no gluon content-can be used to describe hadron physics lies around µ ∼ 0.5 − 1 GeV.
Thus, we evolve the quark DAs obtained within the nlNJL from our input scale µ 0 = 3 GeV to a lower energy scale, which we choose to beμ 0 = 0.5 GeV, and at this lower scale we impose the condition of no gluons. Then, for higher values of µ, we allow gluon contributions to be generated through the evolution equations, which mix quark singlet and gluon components of the DAs. For the second analysis, once again we proceed by using the nlNJL quark model parametrization in order to calculate the coefficients a Taking into account the results of our first analysis (discussed in the previous subsection), in which we obtain χ 2 /(number of points) = 1.33, it is seen that then = 2 fit shows no gain of quality in the description of the experimental data. In addition, then = 4 fit leads to a strong cancellation between the n = 2 and n = 4 gluon coefficients. There is no physical reason for this cancellation, therefore we interpret this result as a spurious solution. Thus we conclude that there is no evidence of a significant presence of gluons in the case of the η meson.
For the η ′ meson then = 2 fit leads to a Contrary to the case of the η-γ TFF, now we observe that there is a significant gain of quality in the description of the experimental values in comparison with the results from our first analysis and with those quoted in Sec. III. We recall that the latter, obtained under the assumption of no gluon contributions to the η ′ DA, lead to a fit of η ′ -γ TFF with χ 2 /n = 2.5 (see Table VII). Moreover, although then = 4 fit has one more free parameter with respect to the casen = 2, the theoretical description of the η ′ -γ TFF is approximately the same in both cases. Our result is shown by the dashed line in Fig. 3 (n = 2 andn = 4 fits are indistinguishable). For comparison we also include in the figure the previous NLO result with no gluon contribution (full line, indetermination indicated by the grey band), the "asymptotic behavior", according to the definition in Sec. III (short-dashed line), and the asymptotic Q 2 → ∞ value (dotted line). Our analysis shows that the gluon contribution is sizable in the case of the η ′ meson. From the figure it is seen that in the low virtuality region the difference between our NLO calculation and the asymptotic behavior is similar to the difference between the present fit and the NLO result. In fact, the result obtained after considering the fitted gluon contributions to the η ′ DA is comparable to that arising from the inclusion of higher twist contributions, discussed in the previous section.
Finally, it is interesting to compare our results with those obtained in Refs. [19][20][21]. The authors of these articles perform model independent fits of the η-γ and η ′ -γ TFFs, considering only the n = 2 coefficients of the Gegenbauer expansions of the DAs. Moreover, they assume meson independence of the quark and gluon DAs, i.e. they take φ η ′ . In this way they end up with only three free parameters, namely the coefficients a η (′) 2 (1 GeV) = 19 ± 5, as well as the results in Eq. (63) of Ref. [19], which translated to our notation lead to a (q 0 ) η (′) 2 (1 GeV) = −0.12 ± 0.11 and a (g) η (′) 2 (1 GeV) = 18.2 ± 4.5. It is worth noticing that our results do not support the hypothesis of meson independence of quark and gluon DAs, in fact, we find significative differences between them.
Nevertheless, it is seen that the values obtained from our analysis are consistent with the above results. Indeed, considering Eq. (3), and taking values of meson decay constants from Table XII, it is seen that that the coefficients of φ η (′) 2 in Ref. [21] should be compared with our result in Table IX, a (q 8 ) η2 = −0.14, while the results for a (q 0 ) η (′) 2 and a (g) η (′) 2 in Ref. [21] and Ref. [19] are to be compared with our values a (q 0 ) η ′ 2 = −0.06 and a (g) η ′ 2 = 11, see Eq. (49). Taking into account the theoretical and experimental uncertainties, we conclude that the values quoted in Refs. [21] and [19] are compatible with each other and with our results.

V. CONCLUSIONS
In this work we have evaluated the quark DAs for the η and η ′ mesons and the associated η-γ and η ′ -γ TFFs within the framework of a nonlocal Nambu-Jona-Lasinio model. This approach, which has been shown to provide a successful description of various meson observables [30,35], has been previously considered in Ref. [13] for the study of the π meson DA and the π-γ TFF. Since the theoretical framework satisfies all basic symmetry requirements (i.e. chiral, Poincaré and local electromagnetic gauge invariances), the quark DAs turn out to be naturally normalized within this scheme.
One of the main ingredients of our model is the quark propagator, which by construction shows a momentum dependence consistent with lattice QCD results. The calculated quark DAs have to be therefore associated to the momentum scale of lattice data, namely 3 GeV [31]. In general, the comparison of any observable related to the quark DAs (as e.g. the M-γ TFF) with experimental data will require a perturbative evolution of the results obtained at this reference scale. Here we have carried out this evolution up to NLO accuracy in α s , neglecting the mixing between Gegenbauer coefficients of different orders for the singlet quark and gluon DAs.
From the obtained quark DAs at the scale of 3 GeV we observe the following features: (i) whereas our πDA is not far from the asymptotic distribution φ asym (x) = 6x (1 − x), the η and η ′ quark DAs move away from the asymptotic behavior, the departure being larger the larger the meson mass is; (ii) all DAs show two maxima, and this structure arises from the nonlocal genuine contributions in Eq. (25); (iii) in all cases the DAs go to zero rather fast near x = 0 and x = 1, supporting the idea of suppression of the kinematic end points [43,44]. Another outcome of our results is that when the DAs are expanded in Gegenbauer polynomials we find that the absolute values of the corresponding coefficients decrease rather slowly with n, in contrast with usual assumptions.
Concerning the evaluation of the M-γ TFFs, we have found that in general NLO corrections lead to a suppression of Q 2 F (Q 2 ). Although this represents a problem regarding the explanation of the already challenging experimental scenario for the π-γ TFF, the corrections go in the right direction in the case of the η-γ and η ′ -γ TFFs. An important difference between the case of the π meson and those of the η and η ′ mesons is that η and η ′ states can include a gluon-gluon component. In this regard, we have firstly performed an analysis in which these gluon components have been neglected for all Q 2 values, while higher twist corrections have been taken into account by adding 1/Q 2 and 1/Q 4 terms to the dominant twist 2 contribution provided by the DAs. Then we have fitted these contributions to M-γ TFF data. From our results it is seen that the effect of higher twist corrections is more important for the π-γ and η-γ TFFs, particularly for Q 2 3 GeV 2 . Moreover, it is seen that the signs of the corresponding contributions are the same in both cases. Conversely, for the η ′ -γ TFF, contrary to what it should be expected, the higher twist corrections appear to be less concentrated in the low virtuality region.
Finally, we have investigated the effect of two-gluon components of the η and η ′ mesons to leadingtwist accuracy, considering NLO perturbative QCD and neglecting the mixing between Gegenbauer coefficients of different orders. From the numerical analysis it is found that the evolution equations do not generate an appreciable contribution if we assume that meson DAs include no gluons at low virtuality. On the other hand, if we allow for the presence of gluon-gluon components in the η and η ′ DAs at low momentum scales, it is seen that the experimental data for the corresponding As discussed in Sec. IV.B, these results are found to be compatible with previous fits for Gegenbauer coefficients quoted in Refs. [19][20][21]. In this way, from our analysis we conclude that π-γ and η-γ TFFs are more sensible to corrections coming from higher twist effects, while the experimental data on the η ′ -γ TFF points to the presence of a significant gluon-gluon component in the η ′ state.