Virtuality Distributions in Application to $\gamma \gamma^* \to \pi^0$ Transition Form Factor at Handbag Level

We outline basics of a new approach to transverse momentum dependence in hard processes. As an illustration, we consider hard exclusive transition process $\gamma^* \gamma \to \pi^0$ at the handbag level. Our starting point is coordinate representation for matrix elements of operators (in the simplest case, bilocal ${\cal O}(0,z)$) describing a hadron with momentum $p$. Treated as functions of $(pz)$ and $z^2$, they are parametrized through virtuality distribution amplitudes (VDA) $\Phi (x, \sigma)$, with $x$ being Fourier-conjugate to $(pz)$ and $\sigma$ Laplace-conjugate to $z^2$. For intervals with $z^+=0$, we introduce the transverse momentum distribution amplitude (TMDA) $\Psi (x, k_\perp)$, and write it in terms of VDA $\Phi (x, \sigma)$. The results of covariant calculations, written in terms of $\Phi (x, \sigma)$ are converted into expressions involving $\Psi (x, k_\perp)$. Starting with scalar toy models, we extend the analysis onto the case of spin-1/2 quarks and QCD. We propose simple models for soft VDAs/TMDAs, and use them for comparison of handbag results with experimental (BaBar and BELLE) data on the pion transition form factor. We also discuss how one can generate high-$k_\perp$ tails from primordial soft distributions.


Introduction
Analysis of effects due to parton transverse momentum is an important direction in modern studies of hadronic structure. The main effort is to use the transverse-momentum dependence of inclusive processes, such as semi-inclusive deep inelastic scattering (SIDIS) and Drell-Yan pair production, describing their cross sections in terms of transverse momentum dependent distributions (TMDs) f (x, k ⊥ ) [1], which are generalizations of the usual 1-dimensional parton densities f (x). The latter describe the distribution in the fraction x of the longitudinal hadron momentum carried by a parton.
Within the operator product expansion approach (OPE) f (x) is defined [2,3] as a function whose x n moments are proportional to matrix elements of twist-2 operators containing n derivatives (D + ) n in the "longitudinal plus" direction. Analogously, the (k 2 ⊥ ) l moments of TMDs correspond to matrix elements of operators containing the derivative (D 2 ⊥ ) l in the transverse direction. However, in a usual twist decomposition of an original bilocal operatorψ(0) . . . ψ(z) one deals with Lorentz invariant traces ofψ(0)D µ . . . D µ ψ(0) type that correspond to parton distributions in virtuality k 2 rather than transverse momentum k 2 ⊥ . Since D 2 ⊥ is a part of D 2 , it is natural to expect that distributions in transverse momentum are related to distributions in virtuality.
Our goal is to investigate the relationship between distributions in virtuality and distributions in transverse momentum. We find it simpler to start the study with exclusive processes. This allows to avoid complications specific to inclusive processes (like unitarity cuts, fragmentation functions, etc). Also, among hard exclusive reactions, we choose the simplest process of γ * γ → π 0 transition that involves just one hadron. Furthermore, this process was studied both in the light-front formalism [4] and in the covariant OPE approach [5,6,7,8].
In our OPE-type analysis of the γ * γ → π 0 process performed in the present paper, we encounter a TMD-like object, the transverse momentum dependent distribution amplitude (TMDA) Ψ(x, k ⊥ ) that is a 3-dimensional generalization of the pion distribution amplitude ϕ π (x) [9,10,11,12]. The definition of TMDA Ψ(x, k ⊥ ) is similar to that of TMDs f (x, k ⊥ ), and also to that of the pion wave function used in the standard light-front formalism (see, e.g. [4]).
We start in the next section with an analysis of a scalar handbag diagram. We discuss the structure of the relevant bilocal matrix element p|φ(0)φ(z)|0 as a function of Lorentz invariants (pz) and z 2 and introduce the virtuality distribution amplitude (VDA) Φ(x, σ), the basic object of our approach. It describes the distribution of quarks in the pion both in the longitudinal momentum (the variable x is conjugate to (pz)) and in virtuality (the variable σ is conjugate to z 2 ).
The main features of our approach remain intact for the case of spinor quarks, and they also hold when the gluons are treated as gauge particles (Abelian or non-Abelian). For these reasons, we introduce the basic elements of the VDA approach using the simplest scalar example. In particular, we show that the covariantly defined VDA Φ(x, σ) has a simple connection to the impact parameter distribution amplitude ϕ(x, z ⊥ ) (IDA) defined for a spacelike interval z = {z − , z ⊥ }. Then we define TMDA Ψ(x, k ⊥ ) as a Fourier transform of IDA ϕ(x, z ⊥ ).
In Sect. 3, we consider general modifications that appear for spin-1/2 quarks, and then show that the structure of the results does not change if one switches further to gauge theories. Using the parametrization in terms of VDA Φ(x, σ), we calculate the handbag diagram and express the result in terms of the TMDA Ψ(x, k ⊥ ). In Sect. 4, we formulate a few simple models for soft TMDAs, and in Sect. 5 we analyze the application of these models to the pion transition form factor. In QCD, the quark-gluon interactions generate a hard ∼ 1/k 2 ⊥ tail for TMDAs. The basic elements of generating hard tails from soft primordial TMDAs are illustrated on scalar examples in Sect. 6. Our conclusions and directions of further applications of the VDA approach are discussed in Sect. 7.

Choosing representation
We start with analysis of general features of the handbag contribution to the γ * γ → π 0 form factor (see Fig. 1). In the momentum representation, the hadron structure is described by the hadron-parton blob f (k, p), which by Lorentz invariance depends on 3 variables: two parton virtuaities k 2 , (p − k) 2 , and the invariant mass p 2 , which is k-independent. None of these invariants is convenient for extraction of the basic parton variable x which is usually defined as the ratio k + /p + of the "plus" light cone components.
A standard way to extract k + is to incorporate a light-like vector n ≡ n − , which is additional to variables k, p of the hadronparton blob. Using the Sudakov parametrization k = xp+ηn+k ⊥ [13] for the γγ * π 0 amplitude, it is convenient to take n = q , the momentum of the real photon. Another convention was used in the light-front approach of Ref. [4], where n is not directly related to the momenta involved in the process (in particular, q = p + n q 2 ⊥ /(2pn) − q ⊥ in that definition). Switching to the coordinate representation, one deals with the variable z that is Fourier conjugate to k (see Fig.2). The blob now depends on p and z, and by Lorentz invariance is a function of (pz) and z 2 . Then the parton fraction x may be defined just as a variable that is Fourier conjugate to (pz). There is no need to have an external vector like n in such a definition, which is truly process-independent and involves only a minimal set of vectors p, z describing the hadron state under study.
In the coordinate representation (see Fig. 2) we have for a scalar handbag diagram, where D c (z) = −i/4π 2 z 2 is the scalar massless propagator, q is the momentum of the initial virtual "photon" (q 2 ≡ −Q 2 ) given by q = p − q , with p being the momentum of the final "pion".

Twist decomposition of the bilocal operator
The pion structure is described by the matrix element p|φ(0)φ(z)|0 . To parametrize it, one may to wish to start with the Taylor expansion The next step is to to write the tensor z µ 1 . . . z µ n as a sum of products of powers of z 2 and symmetric-traceless combinations {. . . z µ i . . . z µ j . . .} satisfying g µ i µ j {. . . z µ i . . . z µ j . . .} = 0. Using the notation {z∂} n ≡ {z µ 1 . . . z µ n } ∂ µ 1 . . . ∂ µ n for products of traceless tensors, we obtain The operators containing powers of ∂ 2 have higher twist, and their contribution to the light-cone expansion is accompanied by powers of z 2 . Considering the lowest-twist term, one can write matrix elements of the local operators in terms of the coefficients A n . To perform summation over n, one can introduce the twist-2 distribution amplitude (DA) ϕ(x) as a function whosex n moments (x ≡ 1 − x) are related to A n by The calculation of the sum is complicated by traceless combinations {z µ 1 . . . z µ n }. To this end, one can use the inverse expansion {zp} n = (zp) n − 1 4 (n − 1) z 2 p 2 (zp) n−2 + . . . to obtain the parameterization of the matrix element. The plane wave factor e ix(pz) has a natural interpretation that the parton created at point z carries the fractionxp of the pion momentum p. The (z 2 ) k terms brought in by the re-expansion of {zp} n are accompanied by (p 2 ) k factors. These terms are purely kinematical, and may be summed explicitly, leading to an analog of Nachtmann [14,2] corrections. To concentrate on dynamical effects, one may take p 2 = 0, in which case {zp} n = (zp) n , and summation over n is straightforward.

Virtuality distributions
For a light-like momentum p, the O(z 2 ) terms in Eq. (7) only come from the (z 2 ) l terms of the original expansion (3) for φ(z). These terms are accompanied by matrix elements of higher-twist operators The derivation above assumes that the matrix elements of operators containing high powers of ∂ 2 are finite, with their size characterized by some scale Λ, which has an obvious meaning of typical virtuality of parton fields inside the hadron. For l = 0, the coefficients A k0 ≡ A k reduce to those defining the twist-2 DA ϕ(x). In general, for each particular l, we can define the coefficients A kl to be proportional to thex k moments of appropriate functions ϕ l (x), and arrive at parametrization in terms of the bilocal function B(x, z 2 /4).
In fact, using the α-representation as outlined in Refs. [15,16,17], it can be demonstrated that the contribution of any Feynman diagram to p|φ(0)φ(z)|0 can be represented as where d is the space-time dimension, P(c.c.) is the relevant product of the coupling constants, L is the number of loops of the diagram, l is the number of its internal lines, and A(α), B(α), C(α) are positive functions of the α σ -parameters of the diagram. Thus, we can write the matrix element in the form without any assumptions about regularity of the z 2 → 0 limit. There is also no need to assume that p 2 = 0. In a formal Taylor expansion of φ(0)φ(z) above, the z 2 factors are accompanied by ∂ 2 , thus the variable σ is related to ∂ 2 , i.e., parton virtuality. For this reason, we will refer to the spectral function Φ(x, σ) as the virtuality distribution amplitude (VDA). One should be aware, though, that while the virtuality k 2 may be both positive and negative, σ is a positive parameter.
The VDA Φ(x, σ) is related to the bilocal function B(x, z 2 /4) by a Laplace-type representation According to Eq. (10), B(x, z 2 /4) is a function of z 2 − i .

Transverse momentum distributions
The parton picture implies a frame in which p has no transverse component, a large component in the "plus" direction and a small component in the "minus" direction, so as p → p + when p 2 → 0. In the latter limit, only the "minus" component z − is essential in the product (pz). In general, without assuming p 2 = 0, we can take a space-like separation z having z − and z ⊥ components only (i.e., z + = 0), and introduce the impact pa- Note that B(x, β) is defined both for positive and negative β, while ϕ(x, z ⊥ ) corresponds to negative β only. The IDA function may be also treated as a Fourier transform of the transverse momentum dependent distribution amplitude (TMDA) Ψ(x, k ⊥ ). TMDA can be written in terms of VDA as Actually TMDA depends on k 2 ⊥ only: Ψ(x, k ⊥ ) = ψ(x, k 2 ⊥ )/π. This relation is quite general in the sense that it holds even if the z 2 → 0 limit for the matrix element p|φ(0)φ(z)|0 of the bilocal operator is singular. However, if this limit is regular, then the coefficients of the Taylor series of B(x, β) in β are given by σ moments of the VDA Φ(x, σ), which, in turn, are related to k 2 ⊥ moments of TMDA Ψ(x, k ⊥ ), namely For these moments to be finite, Φ(x, σ) should vanish for large σ faster than any power of 1/σ. The functions having this property will be sometimes referred to as "soft". The lowest moment of VDA Φ(x, σ) gives the usual twist-2 distribution amplitude ϕ(x). The reduction relation for TMDA Ψ(x, k ⊥ ) is equivalent to the reduction relation ϕ( For negative β = −z 2 ⊥ /4, this formula may be also obtained by performing the angular integration in Eq. (14). If β is positive (then we can write β = |z| 2 /4), one may understand Eq. (19) as where I 0 is the modified Bessel function. Thus, in some cases, the bilocal function B(x, β) may be expressed in terms of TMDA Ψ(x, k ⊥ ) both for spacelike and timelike values of β. This fact suggests that it might be sufficient to know just the "spacelike" function Ψ(x, k ⊥ ) to describe all virtuality effects. In our actual calculations, we have no need to separate integrations over spacelike and timelike z and use Eqs. (19), (20) and/or the presumptions (like Taylor expansion in z 2 ) on which they are based. All the coordinate integrations are done covariantly, and furthermore, without any assumptions whether the z 2 → 0 limit for B(x, z 2 /4) is regular or not. The results are expressed in terms of VDA Φ(x, σ), and then Eq. (15) is incorporated that relates Φ(x, σ) to the TMDA Ψ(x, k ⊥ ). Our final expressions are given in terms of TMDA.

Handbag diagram in VDA representation
The strategy outlined above may be illustrated by calculation of the Compton amplitude. The starting expression is given by When the partons have zero virtuality, i.e. if all σ l moments with l ≥ 1 vanish, we have B(x, z 2 /4) ⇒ ϕ(x). In this case, only the twist-2 operators contribute, and we have The xxp 2 term here produces target-mass corrections analogous to those analyzed in Refs. [14,2]. For positive Q 2 and p 2 , the denominator has no singularities, and i may be omitted. Also, given the smallness of the pion mass, we will neglect p 2 in what follows. If necessary, the target mass corrections may be reconstructed by the appropriately made changes xQ 2 → xQ 2 + xxm 2 π in the formulas given below. To include the virtuality corrections, we use the VDA parametrization (11) that gives The first term in the brackets does not depend on σ and produces the twist-2 expression (22). Using the TMDA/VDA relation (15), one may rewrite Eq. (23) in terms of TMDA as

TMDA and light-front wave function
One could notice that the definition of TMDA is similar to that for the light-front wave function used in Ref. [4]. The difference is that our approach is based on covariant calculations, with Ψ(x, k ⊥ ) appearing at final stages through its relation (15) to a covariant function Φ(x, σ) which is a generalization of the DA ϕ(x). For this reason, we shall continue to use the "distribution amplitude" terminology for the VDA-related functions. In particular, Ψ(x, k ⊥ ) will be referred to as TMDA.

Non-gauge case
When quarks have spin 1/2, the handbag diagram for the pion transition form factor is given by where S c (z) = / z/2π 2 (z 2 ) 2 is the propagator for a massless fermion. Using that the antisymmetric part of γ ν / z γ µ is iz β µναβ γ 5 γ α and writing the twist-2 part of the matrix element one obtains the same formula as in the scalar model considered above. The higher-twist operators may be included through the VDA parametrization in which we omitted the terms proportional to z α , since they disappear after convolution with z β µναβ . All the formulas relating VDA with IDA and TMDA that were derived in the scalar case are valid without changes. However, the spinor propagator has a different ∼ / z/(z 2 ) 2 functional form. Using it, we obtain For large Q 2 , Eq. (30) contains a power-like ∼ 1/Q 4 correction that corresponds to the twist-4ψγ 5 γ α D 2 ψ operators. Though they are accompanied by a z 2 factor, the latter does not completely cancel the 1/z 4 singularity of the spinor propagator, and this contribution has a "visible" Λ 2 /Q 4 behavior 1 . The remaining term contains contributions "invisible" in the OPE. In terms of TMDA we have With the help of Eq. (31), one can calculate form factor in various models for the TMDA Ψ(x, κ ⊥ ). Before analyzing form factor in some simple models for VDA, we show in the next section that switching to gauge theories does not bring any further changes in general equations.
An important property of this operator is that the Taylor expansion for O α (0, z; A) has precisely the same structure as that for the originalψ(0)γ 5 γ α ψ(z) operator, with the only change that one should use covariant derivatives D µ = ∂ µ − igA µ instead of the ordinary ∂ µ ones: This means that one can introduce the parametrization that accumulates information about higher twist terms in VDA Φ(x, σ) and proceed exactly like in a non-gauge case.

Modeling soft transverse momentum dependence
To analyze predictions of the VDA approach for F(Q 2 ), we will use Eq. (31) incorporating there particular models for VDA Φ(x, σ) or TMDA Ψ(x, k ⊥ ). For illustration, we consider first the simplest case of factorized models Φ(x, σ) = ϕ(x) Φ(σ) or in which the x-dependence and k ⊥ -dependence appear in separate factors, and later consider some non-factorized models.

Gaussian model
It is popular to assume a Gaussian dependence on k ⊥ , In the impact parameter space, one gets IDA that also has a Gaussian dependence on z ⊥ . Writing we see that the integral involves both positive and negative σ, i.e. formally ϕ G (x, z ⊥ ) cannot be written in the VDA representation (15). However, the form factor formula in terms of TMDA (31) shows no peculiarities in case of the Gaussian ansatz, so we will use this model because of its calculational simplicity.

Simple non-Gaussian models
One may also argue that the Gaussian ansatz (39) has too fast a fall-off for large z ⊥ . For comparison, the propagator of a massive particle falls off as ∼ e −|z|m for large space-like distances. At small z 2 , however, the free particle propagator has 1/z 2 singularity while we want a model for p|φ(0)φ(z)|0 that is finite at z = 0. The simplest way is to add a constant term (−4/Λ 2 ) to z 2 in the VDA representation (11). So, we take as a model for VDA. The sign of the Λ 2 term is fixed from the requirement that (4/Λ 2 − z 2 ) −1 should not have singularities for space-like z 2 . The normalization factor p(Λ, m) is given by

m = 0 models
To concentrate on the effects of introducing Λ, let us take m = 0, i.e. consider The bilocal matrix element in this case is given by which corresponds to for IDA. Note that the z 2 ⊥ term of the z ⊥ expansion of ϕ(x, z ⊥ ) in this model was adjusted to coincide with that of the exponential model, so that Λ 2 has the same meaning of the scale of φ∂ 2 φ operator. The TMDA for this ansatz is given by It has a logarithmic singularity for small k ⊥ that reflects a slow ∼ 1/z 2 ⊥ fall-off of ϕ s (x, z ⊥ ) for large z ⊥ . The integrated TMDA that enters the form factor formula (31) is given by It is also possible to calculate explicitly the next k ⊥ integral involved there, see Eq. (56) below.

m 0 model
The model with nonzero mass-like term corresponds to the function that is finite for k ⊥ = 0 in accordance with the fact that the IDA in this case has ∼ e −m|z ⊥ | fall-off for large z ⊥ . For small z ⊥ , we

Modeling transition form factor
Let us now use these models to calculate the transition form factor with the help of Eq.(31).

Gaussian model
In case of the Gaussian model (38), we have For large Q 2 , Eq. (52) displays the power-like twist-4 contribution and the term that falls faster than any power of 1/Q 2 . Note that the x-integral for the purely twist-2 contribution converges if the pion DA ϕ(x) vanishes as any positive power x α for x → 0, while the total integral in Eq. (52) converges even for singular DAs ϕ(x) ∼ x −1+α with arbitrarily small α. Furthermore, the formal Q 2 → 0 limit is finite: where we have used the normalization condition Note that F(Q 2 ) is finite for Q 2 = 0 in any model with finite Ψ(x, k ⊥ = 0). According to Eq. (31), one has then

m = 0 model
Using the non-Gaussian model (44) with m = 0 gives The size of the twist-4 term here is given by the confinement scale Λ 2 , just as in the Gaussian model.

m 0 model
Turning to the m 0 model, we have Now, the size of the twist-4 power correction depends on the interplay of the confinement scale Λ and mass-type scale m.

Comparison with data
In QCD, the twist-2 approximation for F(Q 2 ) in the leading (zeroth) order in α s is Taking the value of I(Q 2 ) ≡ Q 2 F(Q 2 )/ f π from the data gives information about the shape of the pion DA. In particular, for DAs of ϕ r (x) ∼ (xx) r type, one has I LOpQCD r (Q 2 ) = 1 + 2/r, i.e. I as (Q 2 ) = 3 for the "asymptotic" wave function ϕ as (x) = 6 f π xx.
The most recent data [24,25] still show a Q 2 variation of I(Q 2 ) (see Figs. 4,5), especially in case of BaBar data [24] which contain several points with I(Q 2 ) values well above 3. It was argued [26,27] that BaBar data indicate that the pion DA is close to a flat function ϕ flat (x) = f π . The latter corresponds to r = 0, and pQCD gives I flat = ∞. As shown in Ref. [26], inclusion of transverse momentum dependence of the pion wave function in the light-front formula of Ref. [4] (see also [8]) eliminates the divergence at x = 0, and one can produce a curve that fits the BaBar data. Similar curves may be obtained within the VDA approach described in the present paper. In Fig. 4, we compare BaBar data with model curves corresponding to flat DA ϕ(x) = f π and two types of transverse momentum distributions. First, we take the Gaussian model of Eq. (52). A curve closely following the data is obtained for a value of Λ 2 = 0.35 GeV 2 which is larger than the standard estimate Λ 2 = 0.2 GeV 2 [18] for the matrix element of thē ψγ 5 γ α D 2 ψ operator. However, the higher-order pQCD corrections are known [28] to shrink the z ⊥ width of the IDA ϕ(x, z ⊥ ), effectively increasing the observed Λ 2 compared to the primordial value of Λ 2 . For illustration, we also take the non-Gaussian m = 0 model of Eq. (56), to check what happens in case of unrealistically slow ∼ 1/z 2 ⊥ decrease for large z ⊥ . Still, if we take a larger value of Λ 2 = 0.6 GeV 2 , this model produces practically the same curve as the Λ 2 = 0.35 GeV 2 Gaussian model. Data from BELLE [25] give lower values for I, suggesting a non-flat DA. In Fig. 5, we show the curves corresponding to ϕ(x) ∼ f π (xx) 0.4 DA. If we take the Gaussian model (52), a good eye-ball fit to data is produced if we take Λ 2 = 0.3 GeV 2 . Practically the same curve is obtained in the non-Gaussian m = 0 model of Eq. (56) for Λ 2 = 0.4 GeV 2 . Again, a VDA-based analysis of the higher-order Sudakov effects [28] is needed to extract the value of Λ in the primordial TMDA.

Modeling hard tail
Higher-order pQCD corrections also modify the large-k ⊥ behavior of TMDA, producing a hard ∼ 1/k 2 ⊥ tail. Matching the soft and hard parts of transverse momentum distributions is a very important problem in their studies. Below, on simple scalar model examples we illustrate the basics of using the VDA approach for generation of hard tail terms from original purely soft distributions. going out of the point z 1 gives (after integration over z 1 ) for the analog of VDA. For TMDA, this yields It has a hard powerlike 1/k 2 ⊥ tail for large k ⊥ . Making a formal d 2 k ⊥ integration to produce DA, one faces in this case a logarithmic divergence. In the impact parameter space, we have ϕ 2 (x, z ⊥ ) = 2K 0 (z ⊥ m 2 − xxp 2 ), a function with a logarithmic singularity for z ⊥ = 0, which is another manifestation of the divergence of the k ⊥ integral for Ψ 2 (x, k ⊥ ). In fact, the function Ψ point (x, k ⊥ ) has a "bound state" pole in p 2 at the location given by a well-known light-front combination (k 2 ⊥ + m 2 )/xx. However, we see no reasons to expect that in general TMDAs Ψ(x, k ⊥ ) depend on k ⊥ through k 2 ⊥ /xx.

Hard exchange model
A more complicated toy model involves two currents carrying momenta yp and (1 − y)p ≡ȳp at locations z 1 and z 2 , respectively (see Fig. 7). With an exchange interaction described as an analog of VDA, where g is the coupling constant. For p 2 = 0, the β-integral gives a well-known combination that is a part of ERBL [10,4] evolution kernel. For an analog of TMDA in the p 2 = 0 limit, we have A further step is a superposition model in which the yp,ȳp states enter with the weight ϕ 0 (y), a "primordial" distribution amplitude. Then the model TMDA is given by a convolution The integral producing DA in this case converges to give ϕ conv (x) = g 2 /(8π 2 m 2 ) δϕ(x), where that has the meaning of a correction to ϕ 0 (x) generated by the simplest exchange interaction. However, the k 2 ⊥ moment, and all higher k 2 ⊥ moments of Ψ conv (x, k ⊥ ) diverge, which is reflected by ln z 2 ⊥ terms in the expansion of the relevant IDA since aK 1 (a) = 1 + aI 1 (a) ln a + analytic terms for small a.

Generating hard tail
Thus, an exchange of a "gluon" has converted a superposition of collinear (to p) "quark" states into a state that has 1/(k 2 ⊥ + m 2 ) 2 dependence on the transverse momentum k ⊥ . We may also assume that the initial fields at z 1 and z 2 are described by some "primordial" bilocal function B 0 (y, (z 1 − z 2 ) 2 /4) corresponding to a soft TMDA Ψ 0 (y, k ⊥ ) ≡ ψ 0 (x, k 2 ⊥ )/π (see Fig.  8). To concentrate on virtuality effects induced by B 0 , we use The term in square brackets may be written as where ϕ 0 (y) is the primordial distribution Ψ 0 (y, k ⊥ ) integrated over all the transverse momentum plane. Hence, for large k ⊥ , the leading 1/k 4 ⊥ term is determined by the DA ϕ 0 (y) only. A particular shape of the k ⊥ -dependence of the soft TMDA Ψ 0 (y, k ⊥ ) affects only the subleading ∼ [V ⊗ ψ 0 ](x, k 2 ⊥ )/k 2 ⊥ term. The form of k ⊥ dependence of Ψ 0 (y, k ⊥ ) is also essential for the behavior of Ψ B 0 (x, k ⊥ ) term at small k ⊥ . In particular, we have which gives, e.g., ϕ 0 (y)/Λ 2 in the Gaussian model (38).

Hard tail for spin-1/2 quarks
In case of spin-1/2 quarks interacting via a (pseudo)scalar gluon field ("Yukawa" gluon model), Eq. (67) is modified by an extra k 2 ⊥ factor coming from the numerator spinor trace, which leads to 1/(k 2 ⊥ + m 2 ) dependence in the correction (64) to the TMDA. It results in a ϕ conv Y (x, z ⊥ ) ∼ K 0 (mz ⊥ ) δϕ(x) term for the IDA. The logarithmic divergence for z ⊥ = 0 of this outcome corresponds to evolution of the DA. In the B 0 model, we have (switching to ϕ(x, z ⊥ ) → ϕ(x, z 2 ⊥ ) in our notations below) Substituting formally ϕ 0 (y, z 2 ⊥ ) by ϕ 0 (y) in the z 2 ⊥ → 0 limit, we get a logarithmically divergent integral over ν. In fact, for a function ϕ 0 (y, z 2 ⊥ ) that rapidly decreases when z 2 ⊥ 1/Λ 2 , one gets ln(z 2 ⊥ Λ 2 ) as a factor accompanying the convolution of V(x, y) and ϕ 0 (y). Hence, the pion size cut-off contained in the primordial distribution provides the scale in log(z 2 ⊥ ), and we may keep the hard quark propagators massless. This cut-off also results in a finite value of Ψ B 0 Y in the formal k ⊥ → 0 limit: Thus, the Ψ conv Y (x, k ⊥ ) ∼ δϕ(x)/k 2 ⊥ singularity of the "collinear model" Ψ 0 (y, k ⊥ ) = ϕ 0 (y) δ(k 2 ⊥ )/π converts into a constant 1/Λ 2 in the Gaussian model. Note also that the overall factor in Eq. (71) then contains the x-independent integral of ϕ 0 (y), i.e. f π , rather than the convolution δϕ(x) as in Eq. (64).
Concluding, we emphasize that the VDA approach provides an unambiguous prescription of generating hard-tail terms like Ψ B 0 (x, k ⊥ ) from a soft primordial distribution Ψ 0 (y, k ⊥ ). A subject for future studies is to use this strategy for building hard tail models in case of QCD.

Summary and outlook
In the present paper, we outlined a new approach to transverse momentum dependence in hard processes. Its starting point, just like in the OPE formalism, is the use of coordinate representation. At handbag level, the structure of a hadron with momentum p is described by a matrix element of the bilocal operator O(0, z), treated as a function of (pz) and z 2 . It is parametrized through a virtuality distribution Φ(x, σ), in which the variable x is Fourier-conjugate to (pz), and has the usual meaning of a parton momentum fraction. Another parameter, σ, is conjugate to z 2 through an analog of Laplace transform.
Projecting O(0, z) onto a spacelike interval with z + = 0, we introduce transverse momentum distributions Ψ(x, k ⊥ ) and show that they can be written in terms of virtuality distributions Φ(x, σ). This fact opens the possibility to convert the results of covariant calculations, written in terms of Φ(x, σ), into expressions involving Ψ(x, k ⊥ ). This procedure being a crucial feature of our approach, is illustrated in the present paper by its application to hard exclusive transition process γ * γ → π 0 at the handbag level (which is analogous to the 2-body Fock state approximation). Starting with scalar toy models, we then extend the analysis onto the case of spin-1/2 quarks and vector gluons.
We propose a few simple models for soft VDAs/TMDAs, and use them for comparison of VDA results with experimental (BaBar and BELLE) data on the pion transition form factor.
A natural next step is going beyond the handbag approximation. In QCD, an important feature is that quark-gluon interactions generate a hard ∼ 1/k 2 ⊥ tail for TMDAs. To demonstrate the capabilities of the VDA approach in this direction, we describe the basic elements of generating hard tails from soft primordial TMDAs.
Another direction for future studies is an extension of the VDA approach onto inclusive reactions, such as Drell-Yan and SIDIS processes. In particular, we envisage building VDAbased models for soft parts of TMDs that would have a non-Gaussian behavior at large k ⊥ (the need for such models was recently emphasized by several authors [29,30,31]). The VDA approach would also allow to self-consistently generate hard tails from these soft TMDs.