Distribution Amplitudes of Heavy Mesons and Quarkonia on the Light Front

The ladder kernel of the Bethe-Salpeter equation is amended by introducing a different flavor dependence of the dressing functions in the heavy-quark sector. Compared with earlier work this allows for the simultaneous calculation of the mass spectrum and leptonic decay constants of light pseudoscalar mesons, the $D_u$, $D_s$, $B_u$, $B_s$ and $B_c$ mesons and the heavy quarkonia $\eta_c$ and $\eta_b$ within the same framework at a physical pion mass. The corresponding Bethe-Salpeter amplitudes are projected onto the light front and we reconstruct the distribution amplitudes of the mesons in the full theory. A comparison with the first inverse moment of the heavy meson distribution amplitude in heavy quark effective theory is made.


Introduction
The introduction of hadronic light-cone distribution amplitudes (LCDA) dates back to the seminal works on hard exclusive reactions in perturbative QCD [1][2][3][4][5]. These nonperturbative and scale-dependent functions can be understood as the closest relative of quantum mechanical wave functions in quantum field theory. They describe the longitudinal momentum distribution of valence quarks in a hadron in the limit of negligible transverse momentum, here given by the leading Fock-state contribution to its light-front wave function, the so-called leading-twist LCDA. In particular, the light-front formulation of a wave function allows for a probability interpretation of partons not readily accessible in an infinite-body field theory, since particle number is conserved in this frame. In other words, φ(x, µ) expresses the light-front fraction of the hadron's momentum carried by a valence quark.
The slowly establishing consensus, after decade-long controversies about the shape of pion's distribution amplitude, points at a function φ π (x, µ) that is a concave function symmetric about x = 1/2 and much broader than the asymptotic distribution φ(x, µ) µ→∞ − −−− → 6x(1 − x) [19,39], where x is the longitudinal light-front momentum fraction and µ is the renormalization scale. For heavier quarkonia, such as the η c and η b , the distribution amplitudes appear to be increasingly more localized in x ∈ [0, 1] with narrower width and a convex-concave functional behavior [40][41][42]. In the infinite-mass limit these distribution amplitudes tend towards a δ-like function, though this limit is far from being reached at the bottommass scale. The transition from concavely shaped distribution amplitudes into convex-concave ones occurs between the strange and charm quark masses, a mass-scale region known for the onset of important flavor-symmetry breaking effects [43][44][45][46][47].
With respect to factorization approaches in weak heavymeson decays, distribution amplitudes of heavy mesons defined in heavy quark effective theory [48] were for the longest time based on models whose functional form in a given limit is dictated by QCD sum rules [48][49][50][51], guided by an operator product expansion [52] or obtained from a combination of dispersion relations and light-cone QCD sum rules [53]. Additional model approaches exist, see Refs. [54][55][56][57] for instance. A heavy-light LCDA was also extracted from the extrapolation of Bethe-Salpeter amplitudes calculated with an unphysical pion mass [58].
Herein we re-appreciate earlier work on heavy-light mesons and quarkonia [59][60][61][62][63][64][65] within a continuum approach to two-point and four-point functions whose salient features will be summarized in Section 2. The crucial difference in the present approach is the flavor-dependence of the interaction in the ladder truncation of the Bethe-Salpeter kernel, as we effectively take into account that the quark-gluon vertex dressing has a different impact for a light quark than for a charm or bottom quark. In general, D and B mesons are of particular interest as they offer a rich laboratory to study two limiting mass-scale sectors of QCD with associated emergent approximate symmetries: chiral symmetry in the sector of light quarks where m q Λ QCD and heavy quark symmetry for masses m q Λ QCD [45,66]. In Ref. [67] we applied these considerations to a contact-interaction model to obtain the mass spectrum and decay constants of D mesons. It turned out that introducing different effective couplings, due to unlike dressing effects for light and heavy quarks, in the ladder truncation of the Bethe-Salpeter kernel significantly improved the description of experimental results. This was also observed in Ref. [68].
This logic was subsequently applied to the Bethe-Salpeter equation (BSE) [69,70] with a flavor-dependent infrared component of the interaction model introduced in Ref. [71]. We employ the combined approach of the Dyson-Schwinger equation (DSE) for the quark and BSE with a flavor-dependent, slightly modified interaction to first compute the mass spectrum and weak decay constants of the pseudoscalar π, K, D, D s , B, B s and B c mesons and η c and η b quarkonia in Section 2. In doing so, we fix the light and heavy quark flavors uniquely and solve the DSE on the complex plane using Cauchy's theorem [72,73]. The resulting nonperturbative propagators are then inserted consistently and simultaneously in the BSE for the aforementioned mesons. Since the resulting masses and decay constants are obtained without any extrapolations of eigenvalues or masses, we also obtain the corresponding LCDA at a physical mass with appropriate projections of the Bethe-Salpeter amplitudes on the light front in Sections 3 and 4. These distributions amplitudes are then used to compute their first inverse moment which we compare with values in heavy quark effective theory in Section 5. We wrap up with a brief conclusion in Section 6.

Pseudoscalar Bound States
The calculation of a meson's distribution amplitude requires the knowledge of the wave function of this bound state. We do so regarding the mesons as a continuum bound-state problem described by the homogeneous BSE in leading symmetry-preserving truncation. The solution of this eigenvalue problem yields the mass and the Bethe-Salpeter amplitude (BSA) of the meson which can be projected on the light-front to extract a distribution amplitude. It also allows to obtain the leptonic decay constant which directly tests the wave function normalization of the meson. The main ingredients of the BSE's kernel are the dressed-quark propagators and the dressed effective gluon interaction which are described in the next sections.

Dressed-Quark Propagators
The dressed propagators are solutions of the quark's gap equation which can be obtained from the appropriate DSE for a given flavor. The DSE describes the two-point Green function in terms of a non-linear tower of coupled integral equations, each involving other Green functions, most prominent amongst them the dressed-gluon propagator and the quark-gluon vertex. The DSE for a quark of flavor f reads [74,75] where m bm f is the bare current-quark mass, Z f 1 (µ, Λ) and Z f 2 (µ, Λ) are the vertex and wave-function renormalization constants at the renormalization point µ, respectively. The integral is over the dressed-quark propagator S f (k), the dressed-gluon propagator D µν (q) with momentum q = k − p and the quark-gluon vertex, Γ a µ (k, p) = 1 2 λ a Γ µ (k, p), where the color SU(3) matrices λ a are in the fundamental representation; Λ is a Poincaré-invariant regularization scale, chosen such that Λ µ. The most general Poincaré covariant form of the solutions to Eq. (1) is written in terms of scalar and vector contributions: The scalar and vector dressing functions are σ f s (p 2 ) and σ f v (p 2 ), respectively, whereas Z f (p 2 ) defines the quark's wave function and M f (p 2 ) = B f (p 2 )/A f (p 2 ) is the running mass function. In a subtractive renormalization scheme the two renormalization conditions, are imposed, where m f (µ) is the renormalized currentquark mass related to the bare mass by, 1 Henceforth we employ a Euclidean metric which implies: {γµ, γν } = 2δµν ; γ † µ = γµ; γ5 = γ4γ1γ2γ3, tr[γ4γµγν γργσ] = −4 µνρσ ; σµν = (i/2)[γµ, γν ]; a · b = 4 i=1 aibi; and for a timelike vector Pµ we have ⇒ P 2 < 0. and Z f 4 (µ, Λ) is the renormalization constant associated with the Lagrangian's mass term.
The rainbow-ladder (RL) truncation of the integral equation (1) and of the BSE kernel has proven to be a robust and successful symmetry-preserving approximation of the full tower of equations in QCD when it comes to the description of light ground-state mesons in the isospinnonzero pseudoscalar and in the vector channels as well as of the N , N * and ∆ baryons. The rainbow truncation of the DSE is given by the prescription, where we work in Landau gauge in which the free gluon propagator is transverse [74,76], and G(q 2 ) in an effective interaction model of the gluon and vertex dressing. In essence, the complexity of the nonperturbative quark-gluon vertex is reduced to the one leading Dirac term, where an Abelianized Ward-Green-Takahashi identity, is enforced which leads to Z f 1 = Z f 2 [74]. In perturbation theory this is tantamount to neglecting the contributions of the three-gluon vertex to Γ µ (k, p) and obviously a drastic simplification of the Slavnov-Taylor identity for the quark-gluon vertex, as it implies equality of the renormalization constants of the ghost-gluon vertex and ghost wave function:Z 1 =Z 3 [77][78][79][80][81][82][83]. Furthermore, with the ansatz (6) we introduce an additional factor Z f 2 to ensure multiplicative renormalizability of Eq. (1) and thus the renormalization-point independence [84] of the mass function M f (p 2 ): With this, the constants, Z f 2 (µ, Λ) and Z f 4 (µ, Λ), are determined using Eqs. (3) and (4), respectively, which leads to a non-linear coupled renormalization condition [65]. Their values at µ = 2 GeV are listed in Tab. 1.
The ansatz in Eq. (6) implies that a single dressing function G(q 2 ) describes the joint effects of vertex and gluon dressing in the DSE. While this truncation is effective and successful for light hadrons for reasons elucidated in Ref. [85], the dynamics in open-flavor mesons is dominated by a wide array of energy scales. In particular, the nonperturbative interactions of a light quark and a charm or bottom quark with a gluon cannot be assumed to be similar and therefore be described by equal dressing functions. And whilst the truncation of Eq. (6) preserves the axialvector Ward-Green-Takahashi identity (WTI) and therefore chiral symmetry, it is clear that the identity (8) for a bare vertex, can only hold approximately when This is, at best, the case for very heavy quarks as dynamical chiral symmetry breaking (DCSB) contributes little to their mass functions. Yet, this is far from true for the charm quark [45] and for lighter quarks dressing effects are important.
In order to introduce the flavor asymmetry in the antiquark-quark interaction kernel of the D and B mesons, we assume an explicit flavor dependence in the interaction G f (q 2 ) and denote it by a subscript. Herein, we will use G u (q 2 ) = G d (q 2 ) = G s (q 2 ) = G c (q 2 ) = G b (q 2 ). The dressing function G f (q 2 ) is modeled after Ref. [71] and consists of the sum of an ansatz in the infrared region, which dominates for |k| < Λ QCD and is suppressed at large momenta, and a second term that implements the regular continuation of the perturbative QCD coupling and dominates large momenta, where we deliberately absorb a factor 1/q 2 from the gluon propagator (7) in the definition. The expressions for both terms are given by, with γ m = 12/(33 − 2N f ) being the anomalous dimension and N f the active flavor number, Λ QCD = 0.234 GeV, τ = e 2 −1, F(q 2 ) = [1−exp(−q 2 /4m 2 t )]/q 2 and m t = 0.5 GeV. The ansatz in Eq. (12) can be parametrized by where m 2 g (q 2 ) is an effective gluon mass that vanishes in the ultraviolet and M g is a mass scale [86]. The low- momentum component leads to an infrared massive and finite interaction consistent with modern DSE and lattice-QCD results and is responsible for DCSB. Hence, the flavor dependence is explicit in the infrared component of the interaction via the constant, where we use equal values of κ f and ω f for f = u, d, s and different ones for κ c , ω c and κ b , ω b ; note that κ f is in unit of GeV. For comparison, we plot G f (q 2 )/q 2 in Fig. 1, from which it is clear that the interaction is strongly attenuated in the heavy sector and the interaction probes more the light quarks in the infrared domain. The interaction strengths of the heavy quarks overlap within the uncertainty bands due to ∆ω f discussed in Section 2.3, albeit that of the charm quark is slightly more suppressed contrary to expectation. Indeed, G b (q 2 )/q 2 can be made weaker than G c (q 2 )/q 2 with readjustments of ω b and κ b , while keeping the mass spectrum of the B, B s , B c and η b virtually the same. However, the decay constants of the B mesons suffer a decrease of 20 − 30%. With this interaction ansatz, we solve the DSE for each quark flavor at space-like momenta, p 2 > 0, using the parameters reported in Tab. 1. These parameters have been chosen so to reproduce the masses and decay constants studied herein. More precisely, m u (µ) = m d (µ) and κ u and ω u are set with the pion's mass and decay constant, m s (µ) is fixed likewise with the kaon using κ s = κ u , ω s = ω u , and these same light-quark parameters are employed in the heavier mesons. Similarly, the values of m c (µ), κ c and ω c are chosen so to reproduce m D and f D , whereas m η b and f η b are obtained from adjusting m b (µ), κ b and ω b . We fix the quark masses at µ = 19 GeV and then evolve them to a scale of 2 GeV at which we compute the quark propagators in the complex plane, as will be discussed shortly in Section 2.3.

Bound-State Equation
The axialvector WTI is crucial to satisfy the chiral properties of the Goldstone bosons of QCD and to guarantee that the pion is massless in the chiral limit. The identity is derived from chiral transformations and reads [87], where Γ f g 5µ (k; P ) and Γ f g 5 (k; P ) respectively denote the axialvector and pseudoscalar vertices for two quark flavors, f and g, and P is the total four-momentum of the meson, P 2 = −m 2 M . The short notation for the quark momenta, k η = k + ηP and kη = k −ηP , defines momentum-fraction parameters, η +η = 1, η ∈ [0, 1].
In order to constrain the Bethe-Salpeter kernel by the quark propagators S f (k), the interaction and the ansatz for the quark-gluon vertex (9), one inserts the DSE (1) as well as the axialvector and pseudoscalar vertices given by, in which K f g (q, k; P ) is the fully-amputated quark-antiquark scattering kernel and the Dirac-and color-matrix indices are implicit, into Eq. (15) which leads to the rela- where we define: Closely inspecting both sides of Eq. (18) one realizes that, in a RL truncation with a flavor-dependent interaction, the kernel K f g (q, k; P ) on the left-hand side must express an average of the interactions. Note that in the limit ∆ f µν (l) = ∆ g µν (l) the identity Eq. (18) is satisfied by the usual RL kernel, In a consistent ansatz for K f g (q, k; P ) that satisfies Eq. (18), it can be shown [70] that the kernel behaves for large momenta q → ∞ as, whereas in the infrared limit this becomes, In both cases the kernel tends to an average of interaction functions, in the latter case weighted with flavored quarkdressing functions.
In this light, we choose the ansatz for the kernel, combining the wave-function renormalization constant of both quarks where the averaged interaction in the low-momentum domain is described by: We insert this ansatz in the homogeneous BSE, and obtain Poincaré-invariant solutions which are the BSAs, Γ f g M (k, P ), in the pseudoscalar channel J P C = 0 −+ . They can be expanded in a non-orthogonal base with respect to the Dirac trace: We remind, though, that mesons with unequal quarks, such as the kaon, D and B mesons, are not eigenstates of the charge-conjugation operator defined as, Thus,Γ M (k, P ) = λ c Γ M (k, P ) does not imply λ c = ±1 for their charge parity. For mesons made of valence quarks with equal current mass and J P C = 0 −+ , the constraint that the Dirac base in (27) satisfies λ c = +1 requires the scalar amplitudes to be even under k · P → −k · P . Each amplitude, in which F 0,1 i (k, P ) are even under k · P → −k · P . As a consequence, the neutral pion and quarkonia have the property F 1 (k, P ) ≡ 0, yet F 1 (k, P ) will contribute to flavored mesons.
We apply the Nakanishi normalization condition [88,89] which makes use of the eigenvalue trajectory λ(P 2 ) of the BSE, Fig. 2. The integration domain defined by the BSE for a meson of mass mM described by a parabola in the complex q 2 η -plane.
to normalize the BSA and verify the normalization with the usual canonical method: The normalization is required to calculate the weak decay constant of the pseudoscalar meson: defines the Bethe-Salpeter wave function. As already noted, the quark momenta, k η = k + ηP and kη = k −ηP , define momentum-fraction parameters; no observables can depend on them owing to Poincaré covariance.
The weak decay constant may also be inferred from the Gell-Mann-Oakes-Renner (GMOR) relation which is just a different expression of the axialvector WTI that describes the axialvector-current conservation in the chiral limit; see for instance Refs. [59,69] for details of the calculation. Comparing the decay constant obtained with Eq. (32) and with the GMOR relation provides us an additional check of the kernel in Eq. (23) and we find variations for f D and f B of the order of 3 %.

Numerical Results on the Complex Plane
The numerical solution of the BSE (26) implies the knowledge of the quark propagator, and likewise for S f (qη). In Euclidean space, the arguments q 2 η and q 2 η define parabolas on the complex plane,  where z = q ·P/|q||P |, −1 ≤ z ≤ +1, is an angle. In Fig. 2, we illustrate a typical situation encountered in numerical studies of the DSE with an external time-like momentum. In this figure, the points (−η 2 m 2 M , 0) and (0, ±2η 2 m 2 M ) are respectively the intersection points with the real and imaginary axis defined by, and similarly forη. We note that the size of the parabolic domain is determined by the meson mass m M and the parabola is fully defined once m M and η are known. Since the entire interior of the parabola is sampled in the numerical integration of the BSE, singularities inside this domain must be avoided. For the simple case of the pion, where m M < 1 GeV, these complex-conjugate singularities will be outside the parabolic region, as illustrated with the symbol "×" in Fig. 2 2 . In this case, the propagator functions σ f s,v are analytic and Cauchy's integral theorem can be applied. On the other hand, with increasing meson mass, such as for the D where m D > 1 GeV, singularities may lie within the integration domain and the parameters η andη can be chosen to adapt the parabola size and shift external momentum from one constituent propagator to the other. In doing so, there is a limiting bound-state mass for which the integration domain starts to overlap these singularities.
The numerical application of Cauchy's integral theorem we employ is explained in detail in Ref. [73] and our parametrization of the parabola contour is described in Ref. [59]. In particular, we use a distribution of momenta on the contour that is skewed towards the vertex of the parabola. As an example, we plot the real part of σ u s (q 2 η ) on the complex plane in the left-hand pannel of Fig. 3, where the maximal hadron mass that can be reached is m max M = 0.2 GeV > m π . Clearly, the dressing function is analytical in this complex domain.
In case of the B mesons we consider, their large masses do not allow for an optimized η andη pair that produces a parabola free of singularities in the b-quark propagator and in the light(er)-quark propagator. We therefore  resort to a complex-conjugate pole representation of the b-propagator for these mesons. The BSA of the η b , on the other hand, is obtained with the propagator solution on the complex plane, as the conjugate-complex singularities remain outside the parabolic region which can be inferred from the left-hand panel of Fig. 4. Hence in order to compute the static properties of ground-state B u,s,c mesons, we combine two approaches. Due to the issue of unavoidable singularities in the b-quark propagator on the complex plane, we implement a complex conjugate pole (ccp) parametrization for the heavy-quark propagator, while for the u, s and c quarks we use their solutions on the complex plane. The ccp parametrization is given by, m f k and z f k being complex numbers. These parameters are fitted to the DSE solution (2) for N = 2 on the real spacelike axis p 2 ∈ [0, ∞), and the thus obtained ccp representation is then analytically extended to complex momenta. Since we make use of the this representation to calculate the Mellin moments (45) of the LCDA, we list their parameters in Tab. 2 for all flavors.
Obviously, we want to make sure that these parametrizations constitute a realistic reproduction of the dressing functions on the complex plane. To this end, we define the function, where the superscripts cp and 2ccp denote respectively numerical solutions on the complex plane and solutions using Eq. (38) with two complex conjugate poles and the fitted parameters in Tab. 2. As can be seen in Figs. (3) and (4), the deviations ∆σ s (q 2 ) are noticeable near the vertex of the parabola, yet the scale of these variations is dwarfed by the magnitude of σ cp s (q 2 ). We also note that the weak decay constant of the pion calculated with the 2ccp approach differs by only 3% from the cp result; similar observations hold for the kaon and D mesons and the η b mass is almost equal with either method, while the decay constant differs by 6%. We conclude that the use of the 2ccp bottom-propagator is a reliable approach and gives us confidence to calculate the B meson's static properties.
Our results for the masses and leptonic decay constants of the ground-state pseudoscalar mesons are listed in Tab. 3, from which it becomes clear that they are in very good agreement with experimental data when available or lattice-QCD results otherwise. In this table, we exclude the B mesons, the reason for which is that the above mentioned hybrid approach is employed. The results for the B mesons as well as for the η b using the hybrid approach are found in Tab. 4. The masses of these mesons are in excellent agreement with experimental values, while our decay constants compare reasonably well with simulations of lattice-QCD.
The theoretical uncertainties are obtained as follows: in adjusting the dressing function of the interaction (12) in the light-meson sector, we set the scale with the pion and kaon masses and weak decay constants. As well known, these observables are rather insensitive to a range ω ± ∆ω and we set an upper and lower limit, 0.45 ≤ ω u,s ≤ 0.55 GeV, about the central value ω u,s = 0.5 GeV, as depicted by the uncertainty bands in Fig. 1. Having introduced this uncertainty in the light sector, the repercussions are immediate in computing the properties of the D u and B u ; their mass uncertainties are due to the lowenergy scale insensitivity of ω u,s . Likewise, we observe the sensitivity of D u to variations of 10% in ω c and this yields and error estimate for the D s , η c and B c . Finally, we fix the bottom quark at the η b mass scale and check the combined effect of permissible variations of ω b in the BSE of the B s , B c and η b that ensure the η b mass stays within 1% of its central mass value.
We remind that these results are not achieved without the implementation of the flavor dependence of the interaction in Eqs. (12) and (25). The flavor dependence is important to accommodate the fact that heavy quarks probe shorter distances than the light quarks at the corresponding quark-gluon vertices, thereby implying a smaller coupling strength for heavy quarks.

Distribution Amplitudes
A unique leading-twist LCDA exists for any pseudoscalar meson M with total momentum P and is defined in QCD via a meson-to-vacuum matrix element of a nonlocal antiquark-quark light-ray operator as, 0|q f (y 2 n)W [y 2 n, y 1 n] γ · n γ 5 q g (y 1 n)|M (P ) where f M is the weak decay constant of the pseudoscalar meson, n is an auxiliary light-like four-vector with n 2 = 0, x = k z /P z is the momentum fraction of the quark in the infinite-momentum frame withx = 1 − x, y 1 and y 2 are real numbers, and W [y 2 n, y 1 n] is a light-like Wilson line connecting the quark fields q f and q g to produce gauge-invariant quantities for any choice of y 2 and y 1 . The momentum-space distribution amplitude φ M (x, µ) is the Fourier transformed distributionφ M (y, µ) in coordinate space.
In principle, φ M (x, µ) is directly accessible from the light-front wave function [93] by integrating over the meson's transverse momentum, where ψ M (x, k ⊥ ) is the Fourier transform of the positiveenergy projection of the Bethe-Salpeter wave function evaluated at equal time, y + = y 3 + y 0 = 0, in coordinate space [5].
In calculating the BSE in momentum space, however, we work in Euclidean space amenable to numerical calculations. We therefore do not have direct access to the light-front wave function ψ M (x, k ⊥ ). A method to eschew the calculation of ψ M (x, k ⊥ ) consists of computing instead Mellin moments, using Eq. (40) and then reconstructing the LCDA from these moments. In particular, the zeroth moment serves to normalize the distribution amplitude and we choose, In order to make use of Eq. (40), we need to Fourier transform the matrix element to momentum space where after appropriate use of the LSZ reduction formula it can be expressed as the light-front projection of the Bethe-Salpeter wave function χ M (k η , kη), The moments x m and therefore reconstructed distribution amplitudes are valid at a given scale at which the BSA was calculated. All our results are given for a fixed scale: µ = 2 GeV. To conclude this section, we note again that the definition in Eq. (40) contains a Wilson line between the points y 1 and y 2 . In light-cone gauge this operator is trivial: W [y 2 n, y 1 n] ≡ 1, though the implementation of this gauge in numerical approaches to bound-state equations is currently impracticable. On the other hand, it has been argued within a nonperturbative instanton vacuum approach [94] that the contribution of this gauge link to the leading twist-2 quark operators is suppressed. With this in mind, we omit the contribution of the link W [y 2 n, y 1 n] and postpone a nonperturbative approach to the matrix element.

Pion Distribution Amplitude
In order to reconstruct φ M (x, µ) from the moments we write it in terms of Gegenbauer polynomials, C α n (2x − 1), of order α which form a complete orthonormal set on x ∈ [0, 1] with respect to the measure [x(1 − x)] α−1/2 . As argued in Ref. [9], the common projection of φ M (x, µ) on a C 3/2 n [n = 0, ..., ∞] basis comes at the cost of a large number of terms in the Gegenbauer expansion. It turns out to be more economic to consider α itself a parameter which allows to limit the expansion to two terms for the pion, wherex = 1 − x and the normalization is given by, and proceed as follows: we reconstruct the LCDA by minimizing the function, with the moments x m obtained by means of Eq. (45) and x m rec. using the definition (42) and the expansion of Eq. (46). We remind that in the asymptotic limit the LCDA tends to [4]:

Kaon Distribution Amplitude
Flavored mesons like the kaon are composed of valence quarks with different masses and are not eigenstates of charge conjugation. This quark-mass asymmetry reflects in the distribution amplitudes: φ K (x) = φ K (1 − x). In order to adapt the method described above to reconstruct the LCDA to unequal-mass mesons, we define moments in terms of the difference of the momentum fractions denoted by, and define the moments, We thus reconstruct the kaon's LCDA with the parity decomposition, where we employ one and two Gegenbauer polynomials, respectively, in the even and odd components, and N (α) and N (β) are both as in Eq. (47). The even and odd components of the distribution amplitudes are then determined independently by separately minimizing, where the reconstructed moments ξ m E,O rec. are obtained with the distribution amplitudes in Eqs. (54a) and (54b).

Heavy Mesons and Quarkonia
The D and B mesons and heavy quarkonia are treated similarly, yet we employ a different functional form for φ rec. H given by [42], where the normalization is, using the definition of the error function Erf(x) = 2 √ π x 0 dt e t 2 : The reason for this choice is that the Gegenbauer procedure sketched above is appropriate for broader and concave amplitudes, whereas a distribution amplitude with a convex-concave behavior of functions reminiscent of the δ-function in the infinite-mass limit is more appropriately described by Eq. (57). This functional form of the distribution amplitude for heavy quarkonia is also found in the application of the maximum entropy method to extract the Nakanishi weight function of the quarkonia's Bethe-Salpeter wave function [95]. We verify the validity of our reconstruction with the simple polynomial ansatz, φ rec.
H (x, µ) = N (α, β)x α (1−x) β and observe that over the entire range, x ∈ [0, 1], the LCDAs reconstructed either way are but indistinguishable. The uncertainty in reconstructing the LCDA is therefore much smaller than that due to the model parameter ω f . On the other hand, using the separation in even and odd components with Eqs. (54a) and (54b) in case of the D and B mesons requires the computation of a large number of Mellin moments to fix their coefficients. The larger moments suffer numerical instabilities for these heavy-light mesons and we thus prefer the representation in Eq. (57).
We reconstruct the LCDA as in Section 3.1 by minimizing, with x m rec. calculated as described before and making use of Eq. (57).

Mellin Moments
The Mellin moments x m are integrals over a BSA and quark propagators. We follow Ref. [9] in using Nakanishitype representations of the scalar BSA amplitudes F i = E M , F M , G M , H M , and likewise use the 2ccp propagators (38) which allows us to represent the moments in Eq. (45) by Feynman integrals. The amplitudes F i for equal-valence quark mesons are therefore parametrized by, with the definitions, , and the spectral density is given by, The scalar amplitude H(k, P ) is negligibly small, has little impact, and is thus neglected. We do not fit the amplitudes directly but rather the Chebyshev moments F m i (k, P ) of the expansion, where z p = k ·P/|k||P | and we typically use m = 4 Chebyshev polynomials U m (z p ). The fit parameters for mesons with equal valence-quark masses, namely pi, η c and η b , are tabulated in Appendix A. We here report the first four Mellin moments computed with Eq. (45) combining the 2ccp parametrization for the quark propagators in Tab. 2 and the Nakanishi representation of the BSAs introduced above.
x m π 0.500 0.318±0.008 0.228±0.006 x m ηc 0.500 0.273±0.001 0.160±0.001 x m η b 0.500 0.262±0.001 0.144±0.001 In the case of flavored mesons with unequal quark masses, a satisfactory representation of the numerical solutions to the scalar function of the BSAs is [96], The sets of parameters that fit the BSAs of flavored mesons are listed in Appendix A. We calculated the Mellin moments of the kaon ξ m K using Eq. (45), where we remind that the moments are in terms of ξ = 2x − 1, as we separated the even and odd moments to reconstruct the LCDA. The first three moments are: On the other hand, we compute the moments x m M of the D and B mesons since we do not separate the even and odd components of the distribution amplitude by means of Gegenbauer polynomials. In principle, this can also be achieved for, yet while we managed to obtain a reasonable fit for the D and D s , the solutions for the B mesons are numerically unstable. This is not unexpected, as the distribution amplitudes of the heavy-flavored mesons reveal a pronounced asymmetry with respect to x → (1 − x) not easily reproduced with just a few Gegenbauer polynomials. The BSA of the D and B mesons are fitted to the Nakanishi-like representation in Eq. (67) and also tabulated in Appendix A. To reconstruct the heavy-meson LCDA only these three moments are needed.

Reconstructed Distribution Amplitudes
We are now able to present numerical results for the reconstructed LCDAs of the pseudoscalar mesons discussed in Section 3. The economic form of Eq. (46) limited to two terms in the Gegenbauer expansion can be fitted with m max = 50 moments, x m π , and yields the parameters: The heavy quarkonia and heavy-light mesons, for which we use an exponential parametrization for the distribution amplitude (57) and m max = 3 moments, are described with the following parameter set: The theoretical uncertainties are due to ω f variations as described in Section 2.3.
In the left panel of Fig. 5 we observe that φ π (x, µ), is concave, symmetric and much broader than the asymptotic limit ϕ asy (x) as a consequence of DCSB. The symmetric shape of the pion's LCDA is precisely due to the fact that this meson is made up of two quarks of the same flavor, each carrying the same amount of momentum fraction of the bound state on the light front. On the other hand, φ K (x, µ) turns out to be equally concave, yet its functional form is characterized by an asymmetric shift toward a peak at x = 0.61. This is a clear sign of dynamical SU(3) flavor-symmetry breaking, where the heaviest valence quark inside the kaon carries a greater amount of the meson momentum. In this case we have M E u /M E s = 0.73. Moving our attention to mesons with larger asymmetries, the right panel of Fig. 5 shows that the LCDAs of the D u and D s are not anymore concave as a function of x, rather their functional form is convex-concave. The heavier charm carries most of the fraction of the meson's momentum. Moreover, φ Du (x, µ) is slightly more asymmetric and peaks higher than φ Ds (x, µ) which is due to the fact that the mass difference between the strange and charm quarks is smaller, i.e.: M E u /M E c = 0.30 and M E s /M E c = 0.42. This stands in contrast to the LCDA of the η c which is symmetric about the mid-way point x = 1/2, though much more sharply peaked than the asymptotic limit. Its behavior as a function of x can be described by convex-concave-convex .
In Fig. 6 we present a comparison of the LCDAs of the charmonium, bottonium and the different D and B mesons. We note that the B u and B s distributions are extremely asymmetric and that the heavy valence quark inside the B u and B s carries almost all of the meson's momentum. The maxima of φ Bu (x, µ) and φ Bs (x, µ) are located at x = 0.92 and x = 0.90, respectively, whereas those of φ Du (x, µ) and φ Ds (x, µ) are at x = 0.76 and x = 0.63. The situation of the B c is somewhere in between the lighter B u and B s and the quarkonia, and we observe that its LCDA is less dislocated from x > 1/2. The maximum of  32. Finally, we note that φ η b (x, µ) is, as expected, narrower than φ ηc (x, µ).

Matching to Heavy Quark Effective Theory
The matrix element in Eq. (40) imply quark fields in the full theory and therefore give rise to distribution amplitudes in QCD which are not directly related to those in Heavy Quark Effective Field Theory (HQET), e.g. in QCD factorization applied to weak decays of B mesons [30][31][32][33][34][35][36][37][38]. In HQET, the LCDA φ B (x, µ) is defined by [48], where h v is the heavy-quark field in the effective theory. In the heavy-quark limit, m Q → ∞, the velocity of the heavy quark is almost unaffected by the interactions since ∆v = ∆p/m Q . The interaction with a light quark alters its on-shell four-momentum, p µ = m Q v µ , to an off-shell momentum p µ = m Q v µ + k µ , where k ∼ Λ QCD is the residual momentum. In this limit, the heavy-quark propagator reads at leading order, where v µ is a time-like unit vector, v 2 = 1, for instance v µ = (1, 0) in the heavy quark's rest frame. This propagator must then be inserted in the bound-state equation (26) and the resulting BSA is projected on the light front. We hold off this calculation for the time being and turn our attention instead to the inverse moment of the heavy-meson distribution amplitude, λ H (µ), defined by,  [97], which to leading order in α s (µ) is given by, Here, we use for the running coupling at leading order: In Tab. 5 we list the inverse moments obtained with our LCDAs of the D and B mesons and the corresponding values in HQET at µ = 2 GeV. We remind, however, that the heavy-quark expansion is not reliable in case of charmed mesons and the values for λ(µ) are only presented for completeness. For comparison, QCD sum rules predict λ QCD B (1 GeV) = 0.460 ± 0.110 GeV [49] and λ QCD B = 0.460 ± 0.160 GeV (no scale given) [51], a model LCDA in HQET leads to λ HQET B (2 GeV) = 0.58 ± 0.04 GeV [52], a DSE-BSE approach finds λ QCD B (2 GeV) = 0.54 ± 0.03 GeV [58], whereas a range of 0.2 GeV ≤ λ HQET B (1 GeV) ≤ 0.5 GeV is considered in an analysis of relevant form factors in the decay B → γ ν [53].

Final Remarks
Based on earlier insights in a contact-interaction model of QCD [67], we modify the ladder truncation of the Bethe-Salpeter kernel to take into account the different impact of vertex dressing in case of light and heavy quarks. The usual ladder truncation works very well for heavy quarkonia, yet earlier calculations [59] demonstrated that treating the charm and light quark on equal footing leads to issues with the hermiticity of the interaction kernel for heavy-light mesons.
In essence, we keep the light-meson ladder kernel which preserves the axialvector WTI unchanged, but modify the dressing function of the charm and bottom quark with the ansatz in Eq. (23). This prescription comes at the cost of introducing new interaction parameters, ω c , κ c , ω b and κ b . Nonetheless, it is a justified price to pay not only for its phenomenological success of yielding masses and weak decay constant in very good agreement with experiment, but also due to theoretical considerations. Indeed, in highly asymmetricQq bound states dynamical effects cannot cancel each other to produce a symmetric dressing of both quark-gluon vertices in the BSE, and the interaction strength in the infrared region is strongly suppressed for the charm and bottom quarks. For self-consistency, we verify the values of weak decay constants of the D and B mesons with the GMOR relation, which is an expression of the WTI.
With these results, we project the Bethe-Salpeter amplitudes of the pseudoscalar mesons on the light front and compute moments of the corresponding LCDA. The pion and kaon are reconstructed from these moments with a Gegenbauer expansion, whereas we employ an exponential form of the LCDA for the heavy quarkonia and heavy-light mesons. The latter assumption can be related to the Nakanishi weight function of the Bethe-Salpeter wave function by means of the maximum entropy method, though the almost identical LCDA can be reconstructed from a simple polynomial ansatz.
We stress that our results cannot be obtained without the modified flavor dependence in the heavy-quark sector, in particular numerical calculations of the quark dressing function on the complex plane become feasible as singularities are avoided and our results are valid without any extrapolations. The distribution amplitudes we compute follow the expected pattern, i.e. the pion distribution amplitude is a concave function, much broader than the asymptotic one. The same is observed for the kaon which in addition is not symmetric about the midpoint x = 1/2, a visual expression of SU(3) flavor breaking due to DCSB, and this asymmetry is growing with increasing mass of the heavier quark. The distribution amplitudes of D and B mesons describe a convex-concave function, whereas for the η c and η b the symmetric distribution amplitude is of convex-concave-convex form which tends to a Dirac δ function in the infinite-mass limit.
Eventually, for applications in heavy quark effective theory, these heavy-meson distribution amplitudes must be obtained from Bethe-Salpeter amplitudes in a heavyquark expansion of the charm-or bottom-quark propagator. We have postponed this task for now, but calculated the inverse moments of the heavy-meson distribution which can be related to those in HQET. Of course, the computation of the LCDA in HQET, in particular for B mesons, will be of great interest in reassessing branching fractions of semi-leptonic and non-leptonic decays in factorization approaches.

A Bethe-Salpeter Amplitude Parameters
We here collect all the parameters relative to the BSAs of the pion, η c and η b , and separately the parameters of all flavored mesons, namely the D u , D s , B u , B s and B c . The corresponding parametrizations are found in Section 3.4. Note that these parametrizations correspond to fits to the unnormalized BSAs,Γ f g M (k, P ), and we relate them to the normalized BSA (27) by the normalization, where N M is obtained with Eq. (30) or Eq. (31). Equivalently, we may calculate the Mellin moments with the unnormalized scalar amplitudes F i (k, P ) and apply the condition (43): x 0 ≡ 1.