A Novel Algorithm for Extracting the Parton Distribution Amplitude from the Euclidean Bethe-Salpeter Wave Function

We propose a new numerical method to compute parton distribution amplitude(PDA) from the Euclidean Bethe-Salpeter wave function. The essential step is to extract the weight function in the Nakanishi representation of the Bethe-Salpeter wave function in Euclidean space, which is an ill-posed inversion problem, via the maximum entropy method(MEM). The Nakanishi weight function as well as the corresponding light-front PDA can be well determined. We confirm the previous works on PDA computation therein the different method has been performed.


Introduction.
The parton distribution amplitudes (PDAs) of mesons, demonstrating the meson structure on the light front, play an essential role in the hard exclusive processes. The cross sections of the processes can be written as the convolution of the hard-scattering kernel which can be computed perturbatively, and the PDAs of hadron involved [1]. The leading twist parton distribution amplitudes of mesons are defined by integrating out the transverse momentum k ⊥ from the light front wave function, which then can be obtained through projecting meson's Bethe-Salpeter wave function onto the light-front.
Although some efforts have been made to calculate Bethe-Salpeter equation(BSE) directly in Minkowski space (see, e.g., Ref. [2]) with the simple scattering kernel, many BSE calculations are still carried out in Euclidean space which are much easier to compute. The challenge in the Euclidean scheme is how to project the discrete Euclidean wave function data on the light-front to get light-front quantities. The Nakanishi representation of wave function provides a natural way to solve this problem. One can transform this challenging question to the proposition that if it is possible to compute weight function of Nakanishi representation when one have an appropriate solution of BSE in Euclidean space. Nakanishi representation was proposed in Ref. [3] to parameterize the relativistic twoparticle bound state in Minkowski space. Although there lacking a proof of uniqueness of the weight function in the nonperturbative case in Euclidean space, we suppose the wave function can still be parameterized by the similar form in Euclidean space as where k 2 > 0 is the space-like momentum and P 2 = −M 2 bs with M bs the bound state mass and M is an infrared regulated scale. The weight function g(γ, z) is a two-dimensional function in real space. The corresponding leading twist two-particle lightfront parton distribution can be defined as where n is the light-like vector n 2 = 0. We neglect the possible spin structure for simplification at the moment. With the help of m-order moment defined as 1 0 dxx m ϕ(x), one can get the one-dimensional integral representation of the light-front PDA in Euclidean space as where N is the corresponding normalization constant to ensure 1 0 dxϕ(x) = 1. One can plot the light-front PDA with the complete knowledge of weight function in the Nankanishi representation in hand . The main aim of this paper is then to provide a practical algorithm to construct the weight function with the numerical data of Φ(k, P). We would like to be modest to emphasize that we do not have the ability to prove that whether the weight function is unique definite in a more fundamental view at present. The techniques for projecting a realistic pion's Bethe-Salpeter (BS) amplitude in Euclidean space onto the light-front have been pioneered recently [4], therein a special Nakanishi representation of each scalar component of the BS amplitude was parameterized to produce the corresponding Euclidean functions and the Mellin moments of PDA can be computed directly, which was implemented to reconstruct the PDA.
Besides, this Euclidean Nakanishi representation has also been extended to evaluate the pion elastic form factor [5] and pion valence parton distribution [6]. This technique has been proved successfully in some senses of the integral of k ⊥ physics(PDA for example), practically at least.
On the other hand, for the case of heavy meson system, owing to a damping influence from the large quark mass, Ref. [7] proposes a 'brute-force' approach to calculate the Mellin moments viz. direct integration using interpolations of the numerical solutions for the quark propagator and Bethe-Salpeter amplitude. A damping factor 1/(1 + k 2 r 2 ) m has been introduced for each Mellin moment to reduce the oscillation problem in the integration. Such procedure works well for the low order of moments, but uncertainty increases progressively for the higher order ones. The limit number of Mellin moments always brings large uncertainty in reconstructing PDA. For example, the PDA near x = 1 behavior depends on the higher order moments' behavior. However, Ref. [7] ignored this issue and supposed that the heavy meson PDA exponentially damped at the end-point.
In this letter we provide another numerical method to compute the PDA that is independent of the approach given in Refs. [4] and [7]. The procedure is straightforward: we first extract the weight function of Nakanishi representation (Eq. (1)) via the maximum entropy method (MEM) and then produce the PDA with Eq. (3).

Maximum Entropy Method.
As pointed by the authors of Ref. [8] that the extracting g(γ, z) from the Φ(k, P) is a typical ill-posed inversion problem. They regulate the linear system with a diagonal term ǫI and show the unstability of solution with different value of ǫ. They appeal that more efficient methods are required to obtain a stable and unique solution in the Euclidean space. In this Letter we will claim that the maximum entropy method (MEM) is an appropriate algorithm to solve this problem. In the following we will firstly illustrate this method by an analytical model of weight function and then produce PDAs of π-meson and η c -meson.
The MEM [9][10][11] is an approach that can be used to solve an ill-posed inversion problem, in which the number of data points is much smaller than the number of degrees of freedom available to the function whose reconstruction is sought. Its basis is Bayes' theorem in probability theory [12], which states the probability of an event "A", given that a condition "B" is satisfied: where, within the sample space, P(B|A) is the probability that events of type "A" satisfy the condition "B" (usually denoted as likelihood function); P(A) is the total probability that event "A" can occur (commonly referred to prior probability); P(B) is the total probability that condition "B" is satisfied (playing the role of normalisation constant).
In making use of the MEM to reconstruct the spectral density function, one works with the conditional probability that a spectral function g(γ, z) corresponds to a correlation function Φ(k E ): where M represents the set of all definitions and prior knowledge of the spectral function. According to the central limit theorem, the likelihood functional is usually taken as: where Z L is a normalisation factor; {Φ data (k i,E ), i = 1, . . . , N data } are provided with dynamical approaches and {Φ ρ (k i,E ), i = 1, . . . , N data } are obtained from Eq. (1) using any given model The central feature of the MEM is the prior probability, which can be expressed here in terms of the spectral entropy as where Z S is a normalisation factor, α is a positive-definite scaling factor, and the exponent involves the Shannon-Jaynes entropy [13][14][15] The quantity m(γ, z) is the "default model" of the spectral function, which is usually chosen to be a uniform distribution so as to avoid assumptions about the structure of the spectral density function; viz., A MEM result for g(γ, z) is considered reliable if it does not depend on the choices for the m 0 and Λ. By adding the entropy functional, the solution of the spectral function g is unique, at least the most probable. Thus, with this procedure, people can extract the spectral function from the integration. The application of MEM to extract the quark spectral density at finite temperature to explore the properties of stronginteraction matter successfully (see, e.g. Refs. [16][17][18][19]) show that the MEM is a quite promising method to get the spectral density function. To check further the efficiency of the MEM in exploring the weight function in the Nakanishi representation of the BS wave function, we take the Nakanishi weight function model quoted from Ref. [8] g(γ, z) = e −(γ+1)/(1−z 2 ) , with which we can give the variation behavior of the weight function with respect to the parameters γ and/or z easily. We can also get the g(γ, z) with the MEM described above where the {Φ data (k i,E ), i = 1, . . . , N data } being provided by solving the BSE for mesons and the Dyson-Schwinger equation (DSE) for quark propagator (the calculations will be described in next Section). We use this model to create the mock data of BS amplitude Φ(k i , cos(θ j )) with cos θ = k·P M bs √ k 2 by numerically 2 integrating the Nakanishi representation. The comparison between the numerical result of the Nakanishi weight function we obtained and the demonstration of the above analytical expression is shown in Fig. 1. It is evident that the difference between the numerical result and the analytical one is modest and the point-to-point comparison displays quite good convergence and stability. We can conclude, by this experience, MEM is an appropriate approach to extract the weight function that can solve the problem pointed out in Ref. [8].
3. PDAs of πand η c -mesons in Realistic Model. The success of extracting the weight function via MEM encourage us to compute PDA from an realistic BS wave function. In this section we take the numerical data of pesudoscalar mesons (π and η c ) wave functions that are calculated within the rainbowladder truncation of the BSE and the DSE [20][21][22][23][24][25][26] to achieve our aim.
The meson wave function will be studied by solving the ho-mogenous BSE in ladder truncation where k and P are the qq state's relative and total momenta, respectively, q ± = q ± P/2. The notation stands for a Poincaré invariant regularization of the integral, with Λ the regularization mass-scale. The regularization can be removed at the end of all calculations by taking the limit Λ → ∞. D f αβ represent the free gluon propagator, G denotes the effective interaction and S the quark propagator. This equation has solutions at discrete values of P 2 = −m 2 H , where m H is the meson mass. The equation also determines completely the BS amplitude Γ(k; P) together with the the appropriate normalization condition for the bound states. The normalization of the meson's BS amplitude is usually taken with the condition [27] where N c = 3 is the color number and N J = 2J + 1 is the number of the polarization directions of a meson with angular momentum J. The rainbow truncated DSE for the quark propagator in Euclidean space reads where Z 2 and Z 4 are the wave functions and mass renormalization constant, respectively, m q (µ) is the current quark mass at space-like renormalization point µ. We perform a flavorindependent renormalization scheme as explained in Ref. [6] to define the Z 2 and Z 4 at µ. These equations are consistent and coupled with an effective coupling function G(s) [28], for which we employ the infrared constant Ansatz [29], The first term characterized by the parameters ω and D determines the intermediate-momentum part of the interaction. The second term describes the ultraviolet part and produces the correct one-loop perturbative QCD limit. The leading twist PDA of pseudoscalar meson can be defined by where we have introduced the BS wave function χ 0 − = S Γ 0 − S and the lepton decay constant f 0 − . Two Dirac covariant components contribute to the leading twist PDA, i.e., χ(k; P) = γ · Pχ 1 (k; P) + γ · kχ 2 (k; P) + · · · , where χ 1,2 (k; P) depend on (k 2 ; k · P) and can be determined numerically. As usual, we suppose they can be represented by the Nakanishi form in Eq. (1) and the corresponding weight function would be extracted by the MEM. We can then obtain the PDA with Eq. (3). We have carried out calculations by choosing the same parameters ω = 0.5 GeV and Dω = (0.87 GeV) 3 as the same as those in Ref. [4] to produce the π-meson PDA; and ω = 0.8 GeV and Dω = (0.7 GeV) 3 (the corresponding leptonic decay constant of η c is about 0.28 GeV) to calculate the η c -meson PDA. The presently obtained results of the πand η c -meson PDAs via MEM and the comparison with previous results are illustrated in Fig. 2.
The Fig. 2 shows apparently that the PDAs of the pseudoscalar mesons π and η c presently obtained via the MEM and DSE approach of QCD match the previous results given in the same dynamical method in the valence region very well and the slight difference in the middle region of x is tolerable. Such a good agreement confirms the previous results on one hand and, on the other hand, indicates that the MEM is efficient to determine the PDA.

Summary and Remarks.
In this Letter we propose a practical algorithm to determine the PDA of mesons in the framework of Bethe-Salpeter equation and Dyson-Schwinger equation approach of QCD. The key point of our new algorithm is implementing the MEM to extract the weight function of the Nakanishi representation of the mesons Bethe-Salpeter wave function. The merit of the MEM is that one neither needs to rely on the limit knowledge of the Chebyshev moments of the mesons' Bethe-Salpeter amplitude to parameterize the Nakanishi weight function (like previous π case) by special form, nor has to be restricted by the limit number of Mellin moments (like previous η c case) to suppose some special forms of PDA. The potential advantage of MEM can be applied to find the light-front wave function of meson when one has the Bethe-Salpeter wave function in hand, which we will leave it for future work. In addition, we should also admit some weaknesses of the MEM, for example, much uncertainty would be encountered when the physical weight function is not positive definite that might be the case of excited state's PDA. We also run into some difficulty when the Bethe-Salpeter wave function is not monotonous. However the equivalence of the three methods mentioned above allows us to choose appropriate one to analyze the PDA case-by-case.