Pion wave function from lattice QCD vs. chiral quark models

We analyze the equal-time Bethe-Salpeter quark wave function of the pion obtained from a quenched lattice QCD calculation with delocalized quark interpolators. We find that the result agrees remarkably well with the predictions of the Nambu--Jona-Lasinio model in all channels. We choose the quenched lattice QCD, since it is closer to the large-Nc limit of the Nambu--Jona-Lasinio model. We also show how transversity information, relevant for the light-cone physics, can be obtained from our equal-time rest-frame lattice calculations.


Introduction
In this paper we examine the pion quark wave functions in the Nambu-Jona-Lasinio (NJL) model and in quenched lattice QCD, and confront both results.
Hadronic wave functions encode important information on bound states in strong interaction physics; in particular, they provide the amplitude for a composite hadron to have quarks in a given momentum state or, equivalently, at a certain spacetime distance. For systems with heavy quarks nonrelativistic quantum mechanics applies and particle number is conserved, thus much of our understanding is directly based on wave functions. However, as a matter of principle the wave functions cannot be directly measured experimentally and one must instead resort to form factors, decay widths, or momentum distributions. For light quark systems particle creation may occur, demanding a fieldtheoretic framework where further complications arise; even in the simplest meson case, relativistic invariance requires that one uses the conventional Bethe-Salpeter amplitude [1] with fixed number of quark field operators, a reminiscent of the approximated parton picture point of view, emphasized by the light-cone approaches [2,3,4]. Color gauge invariance requires that one additionally includes link operators [5,6]. For the pion, the spontaneously broken chiral symmetry is a basic dynamical ingre-dient in the determination of its nonperturbative quark structure. It appears via the pertinent axial Ward-Takahashi identities [1,7] (for a review see, e.g., [8] and references therein). These important constraints are implemented in relativistic fieldtheoretic chiral quark models, such as the NJL model [9]. The regularization needs to be carefully handled (for a relevant review see, e.g., [10]).
On the other hand, lattice QCD solves the bound state problem in a fundamental way. It is thus possible to make a first-principle non-perturbative determination but at the expense of breaking continuum symmetries, such as the Lorentz invariance and, quite often, chiral symmetry, due to the finite lattice spacing. After a pioneering study [11] the hadronic wave functions have been analyzed on the lattice on a number of occasions [12,13,14,15,16,17,18,19,20]. The axial Ward-Takahashi identities can be exactly implemented on the discrete Euclidean lattice as shown by Ginsparg and Wilson [21] (see Ref. [22,23] for a recent practical implementation; here we use the same method), enabling realistically small pion masses.

Bethe-Salpeter amplitude
Since we attempt a comparison between a lattice calculation, where only gauge invariant objects are defined, and a quark model calculation with no explicit mention of the color gauge symmetry, some remarks delineating the scope and meaning of such a comparison are in order before presenting the actual calculations.
The Bethe-Salpeter vertex (or the quarkantiquark wave function) of the pion is given by where q(x) is the spinor field operator carrying flavor and color, |π b (q) is the pion state with isospin b and on-shell four-momentum q, q 2 = m 2 π , and p denotes the momentum of the quark field after the Fourier transform. While chiral quark model calculations are naturally formulated in the momentum space, the basic objects in the Euclidean lattice calculations are correlation functions in the coordinate space defined in Sect. 4, which are gauge and renormalization group invariant at all Euclidean times.
The most general form of the quark-antiquark correlator allowed by the symmetries [24,25] has the structure where the wave functions Ψ a , a = P, A, T , depend on the Lorentz-invariant variables x 2 , x · q, and q 2 = m 2 π . The LHS is the inverse Fourier transform of the vertex function in the Bethe-Salpeter amplitude (1), which is finite and undergoes x-independent multiplicative renormalization. Thus, the ratios Ψ(x)/Ψ(0) become cut-off independent as the cutoff is removed, which on the lattice means the lattice spacing a → 0. In other words, we are studying the renormalization-group-invariant object.
The definition (1) is satisfactory for chiral quark models. In QCD, however, it is appropriate only in the fixed point Fock-Schwinger gauge, x µ A a µ (x) = 0 (the special case being the light-cone gauge n µ A a µ (x) = 0), where the standard derivatives, ∂ µ , and the covariant derivatives, D µ = ∂ µ + igA µ , coincide. On the lattice, fixing the gauge has the problem of the Gribov copies, as there exists no complete gauge fixing for nonabelian nonperturbative theories. Moreover, Elizur's theorem prevents non-vanishing vacuum expectation values of gauge variant operators. Non-gauge-invariant bilinear operators are made gauge invariant by joining them with a link operator, however path dependence sets in. 1 Furthermore, gluons carry momentum in the pion and any gauge invariant but path-dependent definition will yield different results (see Ref. [27] for a discussion of various possibilities). For definiteness, we choose the straight line prescription for the path. We also undertake a smearing procedure of the link, as described in Sect. 4.
Another important issue concerns the comparison with either the quenched or dynamical results, with the full inclusion of the fermion determinant. The NJL Lagrangian models to the one-quark-loop level correspond to a large-N c approximation. On the other hand, at large N c the fermion determinant is suppressed in QCD, explaining why mesons are stable in that limit [28,29]. The quenched approximation contains all the leading-N c and a piece of the subleading in N c contributions, which is actually suppressed for heavy quarks. That means that pion loops are 1/N c -suppressed, although not all of the 1/N c contributions come from pion loops [30]. Moreover, besides including the fermion determinant, one should also consider higher Fock states (qq) 2 , etc., in the wave function when comparing to the dynamical lattice results.
In the present paper we consider the quenched lattice results and restrict to the NJL model, as the simplest prototype of a chiral quark model. Our study can be viewed as a useful quantitative test of pion wave functions which can be carried out for other chiral quark models at any level of sophistication. We recall here an early comparison of the instantonmodel hadronic wave functions to lattice results [31].
Taking the charged pion for definiteness of the notation, multiplying (3) with appropriate Dirac and isospin matrices, and taking the traces yields the relations directly used in our evaluation: Our choice of the kinematics is q = (m π , 0) and x = (0, r). In order to extract the instant-form wave functions Ψ a from the lattice we will consider in Sect. 4 the matrix elements with the explicit form of the vertices Γ P = iγ 5 , Γ A = iγ 5 γ 0 /m π and Γ T = γ 5 σ 0i r i /(m π r 2 ).

Wave functions from the NJL model
A convenient way to determine the pion wave function in chiral quark models is by exploiting the axial Ward-Takahashi identity (see, e.g., [10] for the details), relating the quark propagator S(p) and the irreducible vertex function, Γ µ,a A (p + q, p), corresponding to the axial current. The spontaneous breaking of the chiral symmetry generates a constituent quark mass M given by the gap equation, yielding S(p) = i/(/ p − M − m). The pion wave function is extracted from the pion pole of Γ µ,a A in the form of an unamputated vertex function [10], The pion-quark coupling constant, g πqq , satisfies the Goldberger-Treiman relation, g πqq = M/f π , with f π denoting the pion decay constant. The approach is non-renormalizable and requires introducing a finite cut-off which should fulfill a number of requirements concerning the gauge and Lorentz invariance as well as causality. This is actually crucial for many applications, including the determination of Chiral Effective Lagrangeans [32,33], Parton Distribution Functions [34,35], dispersion relations for two point functions [36,37], Parton Distribution Amplitudes and the light-cone wave functions [38,39,40,41], Generalized Parton Distributions [42,43], Generalized Form Factors [44], the photon Distribution Amplitude [45], or the pion-photon Transition Distribution Amplitude [46,47,48,49,50]. For practical calculations the so-called bosonized form is more convenient and we refer to [10] for further details.
We apply the simplest twice-subtracted version of the Pauli-Villars regularization. For an observable A it amounts to the replacement M 2 → M 2 + Λ 2 , followed by the subtraction The model has three parameters, which can be traded for f π , m π , and M . At a fixed value of M we determine m = m 0 and Λ by fixing m π and f π [10] to their physical values m phys π = 139 MeV and f π = 93 MeV. In this work we use M = 300 MeV, Λ = 790 MeV, and m 0 = 8.2 MeV.
The lattice simulations are performed at m π > m phys π , hence we need to increase accordingly the pion mass in the model. At not-too-large values of m π this may be conveniently achieved via the Gell-Mann-Oakes-Renner relation, m 2 π ∝ m, from where m = m 0 (m π /m phys π ) 2 . This value of m is actually taken for the lattice values of m π . Now we evaluate the quantities of our interest, Ψ a (r) ≡ Ψ a (0, −r 2 ). The calculation proceeds according to the diagrams of Fig. 1. Standard Feynman rules yield Ψ a (r) = − d 3 p (2π) 3 e ip·r dp 0 (2π) Tr[Γ a S p g πqq γ 5 S p−q ].
Then we perform the integration over p 0 and carry out the Fourier-Bessel transform over p, with p = |p|. The result is The resulting ratio Ψ a (r)/Ψ a (0), obtained numerically from Eqs. (8), is plotted in Fig. 2. In the chiral limit (m π = 0, m = 0) one can carry out the integration in Eq. (8) analytically, which yields where K 0 and K 1 are the modified Bessel functions. This leads to the following asymptotic behavior at r → ∞: We note an exponential fall-off and a longer tail in the A channel than in the P and T channels. One may also compute the two-dimensional Fourier-Bessel transform of Eq. (9), passing from r to the transverse momentum k T , which then yields the transverse-momentum light-cone wave functions integrated over α. We find in the chiral limit Note that while in the coordinate representation

Wave function from the lattice
We now turn to the quenched lattice calculation of Ψ a (r), defined as (4) and illustrated in Fig. 1(a). The interpolator annihilates quark and anti-quark at given time t and at distance r apart. The Lorentz invariance ensures that Ψ a (r) is independent of the arguments y and t. The gauge link P[G] from y to y + r ensures the gauge invariance and it depends on the choice of the path, as discussed in Sect. 2. We choose a straight path between y and y + r. Since a straight path is unique only along the lattice directions r = N ae i , we evaluate Ψ a (r) only for integer multiples of lattice spacing a, i.e., r = N a. Then P[G] is a product of N gauge links. For a = P, A we use in fact the interpolator (12) averaged over six lattice points at the distance r = N a.
The correlation function is a basic object on the Euclidean lattice, which allows one to extract the pion mass m π , as well as its coupling to a given interpolating composite operator 0|O r a |π . For our purpose we compute a correlation function with a pseudoscalar point source O r=0 P at time 0 and a delocalized sink with a = P , A, or T , O r a at the imaginary time t = −iτ : where we project on the total momentum q = 0. The quantities related to the pion can be extracted after the complete set of physical states is inserted into (13), yielding The ground-state pion at rest dominates at large τ , From a single exponential fit to the lattice correlation functions at large τ we extract m π and w r a for a range r = [0, a, .., 4a]. We use several current quark masses m, which corresponds to m π in the range 345 − 740 MeV.
The extracted pion mass does not numerically depend on r nor the choice of the channel P , A, or T , as expected. The ratios Ψ a (r)/Ψ a (0) for P and A channels are extracted from the identity and plotted in Fig. 2. The wave function Ψ T can not be evaluated at r = 0 on the lattice and we normalize it arbitrarily, such that the NJL model results and lattice values at r = a agree (cf. Fig. 2). Finally, we provide some details of the lattice simulation. We use quenched lattice QCD, since it is closer to the large-N c limit of the NJL model, as explained in Sect. 2. We use 100 gauge configurations generated by the Lüscher-Weisz gauge action [51], as described in [52].
We use so called HYP smearing on the gauge configurations [53]. This replaces a gauge link between the neighboring points on the lattice with a "fat" link. The fat link is a sum of links within hypercubes attached to the original link only. Since this smearing is local, it does not destroy short distance quantities, while it smooths out the large local fluctuations of the gauge field.
The lattice spacing a ≃ 0.148 fm is determined from the Sommer parameter, and the lattice volume is 16 3 ×32. We employ the chirally improved valence quarks [22,23,52], which have good chiral properties.

NJL vs. quenched lattice
Our results are presented in Fig. 2, with the lines showing the model and the points the lattice calculation. In the P and A channels we normalize Ψ a (r) to its value at the origin, while for the T channel we normalize to the model prediction at r = a. We note that the results of the NJL model agree remarkably well with the lattice determination in all three channels. The dependence on m π is very weak. The insensitivity of the lattice wave functions to the value of the pion mass confirms previous findings [14,15,19]. As expected from Eq. (9), the pseudoscalar and tensor wave functions are identical within the error P P P P bars, while the axial wave function is found to be wider. This agrees with the different asymptotic behavior of Eq. (10).
We have also checked that the model results are almost independent of M when this parameter is in the reasonable range, M = 250 − 350 MeV.

Rest-frame kinematics and transversity
We end this paper with general remarks relating the equal-time wave functions presented above to the light-cone wave functions in the impactparameter representation, also known as transver-sity wave functions. Formally, we can write the transformation where α is the Feynman parameter [24]. The presence of a gauge link operator between the quarks is implicit. Comparison to (3) yields by construction As a matter of principle, all scalars Ψ a in Eq. (17) depend on the scalar variables x 2 , q 2 , and x · q, hence we are free to choose any reference frame to evaluate Ψ a . In the rest-frame, or equal-time (ET) kinematics, used in the previous sections of the paper, we take x = (0, r) and q = (m π , 0), whence x 2 = −r 2 and x · q = 0. In the infinite-momentumframe kinematics (q 0 , q) = lim qz →∞ ( m 2 π + q 2 z , q z ) and on the light cone (LC), where x + = 0, one has x · q = q + x − , x 2 = −r 2 . The parameter α acquires the meaning of the light-cone momentum fraction of pion carried by one of the quarks. By comparing the two calculations we find Ψ ET a (0, −r 2 ) = 1 0 dαe i(2α−1)q + x − Φ LC a (α, −r 2 ). For the chosen kinematics q + x − = q · x = 0, hence In other words, our rest frame calculation allows for a direct determination of the transverse-coordinate (impact parameter) dependence of the light-cone wave-function integrated over the α parameter. A similar property for the Generalized Parton Distributions was suggested for the nucleon [54] and the pion [43]. The relation between the equal-time and light-front wave functions has been analyzed recently [55] in the momentum space, where the transversity relation (19) cannot be explicitly seen. It would also be useful to verify the equal timelight cone transversity connection on the lattice. While there exist transverse lattice calculations [56], their focus was on the Distribution Amplitude, Ψ(α, 0) = ϕ(α), leaving out the transverse dependence. Some results where also presented in [57].
The transversity relation (19) is explicitly verified for the NJL light-cone wave function [39].

Conclusions
Here are our main points: -The leading-N c chiral quark model interpretation of the quenched lattice data is not only qualitatively correct, but also remarkably accurate. This is yet another manifestation of the fact that the spontaneously broken chiral symmetry is the key dynamical factor in the pion dynamics. The quality of the agreement suggests that the 1/N c contributions in the quenched smeared lattice simulations are small. -The dependence of both the lattice and model results on the value of the pion mass is very small in the broad tested range m π = 345 − 740 MeV. -The pseudoscalar and tensor quenched pion wave functions are equal on the lattice within the error bars. In the model they are equal in the chiral limit and nearly equal for the used values of m π . -The asymptotic fall-off of the pion wave functions in the model is exponential, ∼ exp(−M r)/r p , with p = 3/2 in the pseudoscalar and tensor channels, while p = 1/2 in the axial channel. This qualitatively complies to the lattice data, where the axial wave function exhibits a longer tail. -A general result, which originates from the Lorentz invariance, concerns the utility of the equal-time rest-frame smeared lattice simulations to determine the transversity information relevant for the light-cone physics. The integrated lightcone (infinite-momentum) pion wave functions in the impact-parameter space, 1 0 dαΦ LC a (α, b), coincide with our equal-time rest-frame wave functions evaluated at the same quark-antiquark separation, Ψ ET a (0, −r 2 ) | r=b .