Comparison of two Minkowski-space approaches to heavy quarkonia

In this work we compare mass spectra and decay constants obtained from two recent, independent, and fully relativistic approaches to the quarkonium bound-state problem: the Basis Light-Front Quantization approach, where light-front wave functions are naturally formulated; and, the Covariant Spectator Theory (CST), based on a reorganization of the Bethe–Salpeter equation. Even though conceptually different, both solutions are obtained in Minkowski space. Comparisons of decay constants for more than ten states of charmonium and bottomonium show favorable agreement between the two approaches as well as with experiment where available. We also apply the Brodsky–Huang–Lepage prescription to convert the CST amplitudes into functions of light-front variables. This provides an ideal opportunity to investigate the similarities and differences at the level of the wave functions. Several qualitative features are observed in remarkable agreement between the two approaches even for the rarely addressed excited states. Leading-twist distribution amplitudes as well as parton distribution functions of heavy quarkonia are also analyzed.


Introduction
In the last decade the renaissance of interest in quarkonium systems has been driven by the discovery of new particles such as the X (3872) state. Ever since, quarkonium spectroscopy has enjoyed an intensive flow of new results provided by B-and charm-factories such as Belle, BaBar and BES III. Also, several experiments at hadron machines such as the LHC can now investigate quarkonium produced promptly in high-energy hadronic collisions, in addition to charmonium produced in B-decays. Data samples with unprecedented statistics are now available and additional data a e-mail: sofia.leitao@tecnico.ulisboa.pt are anticipated with the advent of SuperKEKB, the new B factory at KEK. (A detailed review on the experimental status can be found in [1] and the references therein.) Theoretically, this situation represents an exciting opportunity where QCD-inspired models in parallel with lattice calculations can be extensively tested. In this line of investigation, two recent models of QCD, the Covariant Spectator Theory (CST) and the Basis Light-Front Quantization (BLFQ) obtained a successful description of heavy quarkonia below open flavor thresholds.
In this work we apply improved versions of the models explored in Refs. [2] for CST and [3] for BLFQ respectively, and we extend their range of results with new sets of predictions for mass spectra, decay constants and light-front distributions. We find that both approaches consistently succeed in describing the experimental data, despite the limitations intrinsic to each model.
In CST, a quasi-potential equation is obtained by reorganizing the Bethe-Salpeter equation (BSE) and solving for a given kernel. On the other hand, in BLFQ, solutions are obtained by diagonalizing an effective QCD Hamiltonian. While the goals of both approaches are the same-to formulate a successful relativistic model of the mesons in terms of quark-antiquark degrees of freedom-they are distinct, especially in the way certain features of QCD are implemented. This motivates a detailed comparison of these approaches as well as the resulting observables. We focus on heavy quarkonia where we have the advantage of making useful comparisons at the nonrelativistic limit [4,5].
We emphasize at the outset that in both approaches quarkonium is treated non-perturbatively as a relativistic bound state. Both approaches include a one-gluon exchange interaction, known to be essential for a proper description of the low-lying heavy quarkonia [5,6], as well as a confining interaction.
Another feature that both models possess is the fact that they are formulated directly in Minkowski space-time, making them complementary to other approaches such as Euclidean Dyson-Schwinger Equations (DSE) [7][8][9][10][11][12][13][14][15] and Lattice QCD [16]. In fact, the shared Minkowskian nature of CST and BLFQ invites a detailed comparison. Here we pursue this comparison as follows: first, we take the results for quarkonia of the two approaches and extend them by providing new sets of physical observables obtained with improvements in each approach and presented here in parallel for convenience; secondly, we compare these results and try to address the natural question-is it possible to quantify their similarities and differences not only at the level of observables but also in terms of quantities of interest for further applications, e.g. light-front wave functions?
Though it is well known that establishing a proper connection between Bethe-Salpeter amplitudes or any of its threedimensional reductions, with light-front wave functions is a non-trivial problem, attempts to solve or at least bridge the two approaches have been developed in multiple contexts [17][18][19][20][21][22][23][24][25][26].
Here, and to answer the previous question, we have adopted the widely used Brodsky-Huang-Lepage prescription [27] (cf. [28,29]) to express the CST amplitudes in terms of light-cone variables and compare these amplitudes to lightfront wave functions calculated from BLFQ.
Notwithstanding formal difficulties and caveats, we produce here good agreement between results of the two approaches not only for the lowest states but also for the higher radial and angular excited states. In addition we show that, contrary to what intuition might dictate, a simple map (described below), applied to the CST amplitudes allows us to capture all the qualitative features of genuine light-front wave functions.
We also calculate decay constants using the mapped CST amplitudes within the Hamiltonian light-front formalism and compare them with those calculated directly within the CST approach. The results are consistent, differing roughly by less than 2% for bottomonium and 10% for charmonium.
This good agreement further motivated us to use these mapped CST amplitudes and the BLFQ light-front wave functions for calculating other relevant quantities. Most importantly, we highlight the parton distribution amplitudes (PDAs), whose precise knowledge is crucial for the study of a panoply of processes such as quarkonia production at high-energies [28,30], J/ψ + η c pair production in e + e − annihilation [31], B c → η c transitions [32], decays of heavy S-wave quarkonia into lighter vector mesons [33]; deeply virtual quarkonia production [34] and Higgs boson decays into quarkonia [35].
We also calculate heavy quarkonia parton distribution functions (PDFs) and provide their moments heretofore less studied in the literature. This paper is organized as follows: Sect. 2 introduces the notation, definitions and the general formalism employed in this work; Sect. 3 is devoted to results and discussion and finally, in Sect. 4, we present a brief summary, our conclusions and outlook for further research.

Formalism
In this section we introduce the two models used to study heavy quarkonium that we compare in this work. Motivations, goals and specific details of each approach can be found, respectively, for CST in Refs. [36][37][38][39] and for BLFQ in Refs. [40][41][42][43] and the references therein. Here we aim for a brief self-contained description, accompanied by the definitions of the relevant quantities under discussion in the following sections: decay constants, CST amplitudes, lightfront wave functions (LFWFs), parton distribution amplitudes and parton distribution functions.

Basis light-front quantization for quarkonium
In the light-front Hamiltonian approach, quarkonium is described by state vectors |ψ h . These state vectors can be obtained by diagonalizing the light-front quantized Hamiltonian operatorP − , where P = (P − , P + , P ⊥ ) is the four-momentum; J and m J are the total angular momentum and the magnetic projection, respectively. For a four-vector v, the light-front variables are defined as The eigenvalue equation is usually rewritten as, and so it is easy to identify P + P − − P 2 ⊥ = P μ P μ ≡ H lc as the invariant mass squared operator, also known as the "light-cone Hamiltonian".
In the BLFQ approach [3], an effective Hamiltonian is adopted based on light-front holographic QCD [44]: where m is the mass of the quark; x = p + 1 /P + is the longitudinal momentum fraction of the quark, and is the "soft-wall" light-front holography [44] in the transverse direction. r ⊥ = r 1⊥ − r 2⊥ is the transverse separation of the quark and the antiquark, while is the longitudinal confining potential introduced in Ref. [3]. The partial derivative ∂ x is taken with respect to the holo- The strength of both V T and V L depends on a confinement parameter κ. We have for the one-gluon exchange term with the average four-momentum of the exchanged gluon and C F = 4/3. As an extension to the work presented in [3], instead of using a fixed value for α s , a running coupling is used here (see Ref. [45] for details).

Light-front wave functions
Light-front wave functions are defined from the Fock space expansion of the state vector in Eq. (1). For example, the Fock space representation of a quarkonium state reads Here the coefficients ψ ss/ h (k ⊥ , x) are the valence sector LFWFs with s (s) representing the spin of the quark (antiquark). N c = 3 is the number of colors, k ⊥ ≡ p ⊥ − x P ⊥ is the relative transverse momentum, and |0 is the Fock space vacuum. 1 The quark and antiquark creation operators b † and d † satisfy the canonical anti-commutation relations, The state vector is normalized according to a one-particle state [cf. Eq. (9)]: Then the normalization of the LFWFs reads s,s Note that the state vectors of different particles, e.g., J/ψ and ψ , are orthogonal.

Covariant spectator theory for quarkonium
In CST, quarkonium is described as a relativistic system of a quark and antiquark, bound together by a QCD-inspired interaction. The CST equation can be derived from the BSE. For a bound state of total four-momentum P coupled to a quark with momentum p 1 = p + 1 2 P and an antiquark with momentum − p 2 = −p + 1 2 P, the BSE reads where S(k i ) is in principle the dressed quark propagator. However, in this work, S has been replaced by a propagator with fixed quark mass. The idea of CST is to approximate this equation by keeping in the k 0 -contour integration only the contribution from the positive-energy pole of one quark propagator (for details of this prescription see [36,38,48,49]). This leads to the so-called one-channel spectator equation (1CSE), given by where Θ describe the momentum dependence of the kernel, m is the constituent mass of the quarks, and Note that in the context of CST a "ˆ" over a four-momentum indicates that the particle is on-mass-shell, and that we use k to indicate three-momenta. It is worth mentioning that Eq. (13) retains from the BSE four important properties: manifest covariance, cluster separability, and the correct onebody and nonrelativistic limits. The kernel we employed in the 1CSE consists of a covariant generalization of the linear (lin) confining potential used in Ref. [50], a one-gluon exchange (OGE), and a constant (C) interaction: The Lorentz structure of the confining kernel is flexible: the mixing parameter y allows one to dial between a scalar-pluspseudoscalar structure, which preserves chiral symmetry as shown in Ref. [51], and a vector structure, while leaving the nonrelativistic limit unchanged. An analysis of the asymptotic behavior of large momenta |k| shows that we need to regularize the kernel. We use Pauli-Villars regularization for both the linear and the OGE parts, which yields one additional parameter, the cut-off parameter Λ. We found that Λ = 2m is a reasonable choice.
The momentum-dependent terms of the interaction kernel are where q ( ) = p − k ( ) .
As an extension to the model described in Ref. [2], we now allow the constituent masses to be free parameters. We also let the parameter y to be determined by the fit, in order to examine how much the Lorentz structure of the kernel is constrained by the mass spectra. These parameters are fitted exclusively and simultaneously to bb and cc states.
These amplitudes can be expanded in a very useful basis, wherek ≡ k/|k|, χ λ are two-component helicity spinors and the K ρ j (k) operators are 2 × 2 matrices that depend on the total angular momentum J and the parity P of the meson we study. A list with all K ρ j (k) operators used in this work can be found in Ref. [52].
In terms of these wave function components the pseudoscalar and scalar mesons are normalized according to Similarly, for vector and axial-vector mesons the normalization condition is where M is the mass of the bound state, N c is the number of colors and ψ s , ψ d , ψ p s and ψ p t refer to S-, D-, singlet and triplet P-waves, respectively. Both the total angular momentum J and the parity P are exact quantum numbers of the CST solutions. Expressed as in (20), the CST equation is transformed into a system of coupled partial-wave equations, where each partial wave has a definite orbital angular momentum and total spin. Thus it is straightforward to identify the angular momentum content of each state which is useful when comparing with experiment. On the other hand, the 1CSE solutions do not have a definite charge-conjugation parity, C (for an explicit calculation see [38,48]). However, this can be remedied when the appropriately symmetrized contributions of all four poles of the two quark propagators are taken into account in the k 0 contour integration, and a coupled four-channel equation is solved instead (cf. the discussion in Ref. [2]).

Brodsky-Huang-Lepage prescription
Having specified our models, we now describe our method for converting the CST amplitudes of Eq. (20) into LFWFs.
The covariance of the CST equations allows us to evaluate the longitudinal momentum fraction x of the on-shell quark in the rest-frame (x is an invariant under longitudinal boosts). From the CST kinematics described in Sect. 2.3, the relative four-momentum explicitly reads k = (k 0 , k) = (k 0 , k ⊥ , k 3 ), k ⊥ = |k ⊥ |, and where E k is the on-shell energy (14). In the rest-frame the x variable should in principle be identified as [53] x = k + Consequently, From Eq. (24) one verifies that The latter limit poses a difficulty because x can be outside the region 0 ≤ x ≤ 1. How to properly deal with this issue certainly requires further investigation and it is beyond the scope of this article. Recent work has been done in that direction, investigating the formal relation between the light-cone and CST box diagrams for a scalar theory [53]. In any case, we expect that the contribution to the wave function for values of x > 1 should be small, and that it vanishes exactly in the nonrelativistic limit [17]. We circumvent this difficulty by adopting the Brodsky-Huang-Lepage (BHL) prescription [27], where x is automatically limited between 0 and 1, and investigate to what extent such a prescription gives reasonable results.
For the equal mass case of quarkonium m q = mq = m, the BHL prescription provides From Eq. (27) it is straightforward to derive We thus identify the "CST LFWFs" in terms of the CST amplitudes of Eq. (19) as up to some normalization factors, where s 1 and s 2 are the quark spin projections in z-direction. In Fig. (1) we compare the CST amplitudes for one of the dominant wave function components of J/ψ after using the The CST amplitude on the left panel spreads beyond the physical region 0 ≤ x ≤ 1, but for x > 1 it is fairly small for the illustrative case of charmonium, where relativity is no longer negligible. On the right panel, the wave function mapped using the BHL prescription is symmetric with respect to x = 0.5 and is restricted to 0 ≤ x ≤ 1, consistent with longitudinal light-front momentum conservation.

Definition of physical observables and distribution functions
Both the LFWFs and the CST amplitudes allow us to calculate a variety of observables. But with the LFWFs obtained from the BHL mapping, one also gains direct access to quantities such as light-cone distributions, whose extraction is not as straightforward in approaches relying on Euclidean formulations. In this section we apply the LFWFs to the calculation of decay constants and leading-twist parton distribution amplitudes and parton distribution functions.

Decay constants
Decay constants are very important quantities to probe shortrange physics. In practice they will be sensitive to the effective short-range potential. This implies that, for any realistic model of quarkonia, having a correct implementation of the one-gluon exchange interaction is essential for a good description of the decay constants.
In the absence of a proper renormalization procedure, decay constants could develop dependence on the regularization scheme adopted. By construction, any UV regulator estimate within BLFQ is tied to the basis truncation N max (see Ref. [3] for the detailed definition and for more technical details). In fact, previous studies indicate that the cut-off scale is very well approximated by μ ≡ κ √ N max . On the other hand, in CST there is no dependence on any basis, but α s has been kept fixed in the CST calculations. Furthermore, in CST the regularization of the integral over k in Eq. (13) is governed by the Pauli-Villars cut-off parameter Λ and for that reason it will be taken as the CST estimate for the UV regulator. Later we will come back to this point when analyzing the results obtained for the parton distribution functions. The choice of μ ≈ 1.7m in BLFQ and Λ ≈ 2m in CST, permits a good description of the decay constants, with models just fixed by spectroscopy. For the remainder of the work in BLFQ results this scale cut-off is ensured by choosing N max = 32 for bottomonium and N max = 8 for charmonium, making the scales of the two approaches comparable.
In general, the decay constants for pseudoscalar (P), axialvector (A), scalar (S) and vector (V) mesons are defined, respectively, by the matrix elements where P μ is the total momentum of the meson and the polarization vectors are For pseudoscalar and vector states, BLFQ decay constants are determined as follows: and for CST we have wherek = |k|/(E k + m) and the wave function components ψ(|k|) are implicitly defined in (20).

Leading-twist parton distribution amplitudes
In this work we calculate leading-twist parton distribution amplitudes for pseudoscalar φ P (x) and longitudinally polarized vector φ || V (x) mesons. They are determined through the LFWFs as where f P,V are the previously defined decay constants and μ is related to the renormalization scale or UV cut-off scale. The PDAs defined here satisfy the normalization condition Moments of these distributions are given by In NRQCD, moments are related to the r.m.s. relative velocity of the valence quarks [28] through

Parton distribution functions
Similar to PDAs, parton distribution functions depend on the UV cut-off scale as well. They can be accessed by From the previous discussion, it is known that both approaches already have built-in regulators and for that reason there is no need for a hard cut-off, so we conveniently extend these integrals to infinity in Eqs. (40) and (44) and drop the reference to μ.
The moments of these distributions are given by

Results and discussion
We now proceed to the analysis of the results. This section is organized as follows: first we present a comparison of some relevant physical observables, namely, in Sect. 3.1 the mass spectra, and in Sect. 3.2 the decay constants. In Sect. 3.2, also an important consistency check is performed by computing decay constants directly within the CST approach and comparing to those from the light-front formalism using the CST mapped LFWFs. A discussion as regards the LFWFs themselves is given in Sect. 3.3. We conclude with a calculation of several PDAs and PDFs in Sect. 3.4 for pseudoscalar and vector states.

Mass spectra
The ability to reproduce the mass spectroscopy is a first test of meson models. It is known from the recent work in Refs. [2,3] that both approaches perform well in this regard.
Here, we report on the updated spectra using the latest improvements of each model. A summary of the parameters used is given in Tables 1 and 2. In the case of CST, the tensor mesons (J ≥ 2) have not yet been calculated. Also as discussed in Ref. [2], the axialvector states do not have a definite C-parity. For that reason, they were not included in the fit. We observe a larger deviation from the data for χ c1 and h c than for other states. Both approaches predict consistently similar results for the yet unobserved states such as 3 1 S 0 and 1 3 D 1 and 2 3 D 1 in the bottomonium spectrum. Table 1 BLFQ model parameters. N f is the number of flavors and μ g is the gluon mass used to regularize the integrable Coulomb singularity. κ is the same as in Eqs. (4) and (5), and m is the constituent quark mass. n states is the number of states used in the fit.   Fig. 3 Comparison between the decay constants determined with CST-LFWFs and with BLFQ-LFWFs as well as with other methods and with experiment (PDG). The lattice results are from Refs. [56][57][58][59] and the DSE results are from Ref. [8] Finally, in order to quantify the agreement of the predictions with the experimental measurements, we determined the root-mean-square difference from the measured states below threshold and shown in blue in Fig. 2
In BLFQ the results were obtained with N max = 8 for charmonium and N max = 32 for bottomonium, in order to guarantee the aforementioned UV cut-off of 1.7m. The numerical uncertainty is estimated by varying the scales: for charmonium Δf = 2| f max=8 − f max=16 |; for bottomonium Δf = 4| f max=24 − f max=32 |. As estimates for the numerical uncertainties of the CST calculations we took the difference between the results obtained with 12 and 8 splines.
Overall, both CST and BLFQ results reproduce the data quite well and are consistent with other approaches, where available. In the CST approach the pseudoscalar decay constants are closer to experiment than the vector meson decay constants. It is worth emphasizing that these results are pure predictions in the sense that none of them were included in the fits. For that reason the agreement with experimental data for decay constants (where available) could be improved by incorporating those data in the model parameter fits. The only notable discrepancy comes from the D-wave vector mesons ψ(3770) and 1 3 D 1 where the decay constants are very small in both approaches, as is also observed in some other approaches but not all (see Ref. [60] and the references therein). In principle, these decay constants are only non-vanishing due to the mixing with the S-wave in the nonrelativistic case. However, ψ(3770) is just above the opencharm threshold. It will be interesting to resolve the theoretical speculations about 1 3 D 1 through experimental measurements.
Next we discuss an important test of the map adopted and described in Sect. 2.5. The CST solutions as functions of k as in Eq. (20) are first expressed in terms of x and k ⊥ using Eq. (27), and then we use the light-front definitions of Eqs. (36) and (37) to recalculate the decay constants.
These new "CST-mapped" results (indicated by "CSTmap") are presented in Fig. 3 as well. For bottomonium states, shown in the right panel of Fig. 3, the difference between the two sets of calculations is smaller than the numerical errors. For charmonium, on the left panel, the decay constants determined with the CST-LFWFs are slightly larger. The absolute differences δ are listed in Table  3. The largest deviation (disregarding the D-wave vector states which have tiny decay constants) does not surpass 30 MeV (about 10% in charmonium) and 9 MeV (2% in bottomonium), confirming that the BHL prescription we use works better as one approaches the nonrelativistic limit. Nevertheless, this test provides a reasonable justification for the procedure we follow in CST to obtain heavy quarkonia LFWFs, which we will review in more detail in the next section.  Having new sets of light-front wave functions for quarkonia derived from the CST approach opens the door for several calculations. As already mentioned, the CST equation solved in this work does not respect charge-conjugation symmetry, and thus the CST wave functions do not have a definite C parity. A direct comparison with the BLFQ solutions, which do have definite C parity, allows for a better identification of the axial-vector states obtained from CST. From the BLFQ side, a direct comparison with LFWFs from a different approach also offers benefits. As mentioned earlier, in BLFQ the inevitable basis truncation breaks the rotational symmetry. The total angular momentum J is not well defined and the state identification is based on spectroscopy with the help of P, C, etc. Comparing with CST results, for which J is an exact quantum number, gives guidance to validate this identification. With their rich radial and angular structure, the bottomonium vector meson LFWFs for instance provide a non-trivial test of the methods for identifying J in the BLFQ results.
We investigated LFWFs of all states below open flavor thresholds and with J < 2 (cf. Fig. 2) and for all nonvanishing spin configurations. The obtained wave functions exhibit a close correspondence between CST and BLFQ in their dominant structures for all states and spin alignments.
To visualize the rich structures of the wave functions, we adopt the scheme of Ref. [61]. We note that for a particular polarization λ and spin alignment ss, the LFWFs can be expressed as where k ⊥ = |k ⊥ | and φ = arg k ⊥ . This is valid because the orbital angular momentum projection m = λ − s −s is definite (λ ≡ m J ). In order to visualize these wave functions, we drop the phase exp(im φ), while retaining the relative sign exp(im π) = (−1) m for negative values of k ⊥ . More precisely we plot This scheme essentially takes a slice of the 3D wave function ψ ss (k ⊥ , x) at k y = 0. Let us begin the discussion of the LFWFs with the interesting case of the vector bb because from all the systems this is the one with the largest number of states below its open flavor threshold, the B B threshold. These systems admit a mixture of S-and D-wave components (as long as there is a tensor force). In Fig. 4   and in both, BLFQ-LFWFs and CST-LFWFs, an increasing number of nodes in both transverse (k ⊥ ) and longitudinal (x) directions is observed. As a consequence, and for our particular choice of the coordinate range, a nesting ring pattern emerges. This is consistent with the nonrelativistic interpretation, where the radial excitation is homogeneous in all three directions. Prior to the map described in Eq. (27), CST amplitudes expressed as functions of k show precisely this behavior (see Fig. 3 of Ref. [52]). The 1 3 D 1 wave function resembles the shape of the spherical harmonic Y 0 2 (k). The same happens for 2 3 D 1 , where the complicated inner structure is also compatible with a Y 0 2 (k) but now with an extra node in both k ⊥ and x.
In Fig. 5, in addition to the dominant triplet component, other subdominant components of purely relativistic origin are shown. Here, significant differences appear between the CST and the BLFQ LFWFs. While in BLFQ there is only one ψ ↓↓ component, in CST two extra components compatible with a quantum number = 1 and with spin alignments ψ ↓↓ and spin singlet ψ (↑↓−↓↑) appear and are presented in the last row of Fig. 5. These components emerge from the CST amplitude's ψ ++ component [cf. Eq. (20)] and are absent in BLFQ because there positive-energy and negative-energy states do not mix.
In order to demonstrate how the comparison of the mapped CST with the BLFQ LWFWs can help in the identification of states we give an example: in Fig. 6 we show the dominant components (always labeled "dominant") of the two lowest mass solutions of the bottomonium 1 + states in both formalisms. We observe that the mapped spin triplet CST LFWF resembles very closely the dominant triplet BLFQ component ψ λ=1 ↑↑ (on the left) with J PC = 1 ++ , whereas the spin singlet CST LFWF has its correspondence in the dominant singlet BLFQ component ψ λ=1 (↑↓−↓↑) (on the right) with J PC = 1 +− . Therefore we can conclude that these two CST LFWFs describe mostly the χ b1 (1P) and the h b1 (1P) states, respectively. However, a comparison of the subdominant components shows differences: for instance, one of the BLFQ components of the χ b1 (1P) state, ψ λ=1 ↓↓ , exhibits F-wave features, which, in principle, violates the angular momenta addition |L − S| ≤ J . Such F-wave contributions are not present in the CST LFWFs for states with J = 1. On the other hand, the CST LFWF with the same spin alignment, ψ λ=1 ↓↓ , is a D-wave, which in BLFQ solutions only appears for the h b1 (1P) state. Nevertheless, in the near future, we will solve the more complicated CST equations with chargeconjugation symmetry and some of the observed differences are expected to disappear.

Leading-twist parton distribution amplitudes and parton distribution functions
In this section we present the results for the leading-twist parton distribution amplitudes of 1 S 0 and 3 S 1 states and compare them with other results from the literature. For completeness and in order to compare with other approaches, the moments are displayed in Table 4. In addition, estimates for the root-mean-square velocity of the constituents are displayed in Table 5. In Fig. 7, distribution amplitudes are given for pseudoscalar and vector ground states and their corresponding first radial excitations.
Several global features can be observed in both approaches, such as the number of maxima and minima and the consistent broader curves for charmonium than for bottomonium. We notice, however, that the CST distribution amplitudes exhibit larger tails than the BLFQ distribution amplitudes where a gaussian behavior causes a stronger falloff. This indicates that in the CST approach the constituent quarks are noticeably more relativistic as suggested by the estimate of v 2 (cf. Table 5).  Table 4 Moments ξ n of leading-twist PDAs for pseudoscalar and vector states with longitudinal polarization. The BLFQ results are given at N max = 8 for charmonium and N max = 32 for bottomonium. The other moments correspond to the results from NRQCD [29], QCD sum rules [28] and DSE [15]. The pQCD asymptotic value for ξ n is 3/((n + 1)(n + 3))  In Fig. 8, the longitudinal distribution amplitude φ || (x) of J/ψ is shown for CST and BLFQ, together with the pQCD asymptotic limit given by 6x (1 − x), the light-front holographic QCD applied to charmonium [63], labeled here as AdS/QCD (ii) . For completeness, we also present the longitudinal PDA of the pion calculated from the light-front holo-graphic model of Brodsky and de Téramond [62], labeled AdS/QCD (i) , and in which the PDA equals (8/π ) √ x(1 − x). The CST and BLFQ PDAs are narrower than the latter PDA for the pion, but visibly broader than the prediction from AdS/QCD (ii) for the corresponding J/ψ state. A significant difference is the fact that in this approach, contrary to CST and BLFQ, the one-gluon exchange interaction is absent. This interaction modifies the short-range behavior, to which the PDAs are particularly sensitive to.
The CST and BLFQ PDA curves also differ from the pQCD limit, in principle valid only at μ → ∞. Nevertheless, it is reassuring to confirm that in the CST approach, by increasing the cut-off parameter Λ, indeed the PDAs get a broader shape and smoothly approximate the asymptotic perturbative QCD limit, as shown in Fig. 9 for both η c and η b states. In BLFQ, preliminary studies show that the pseudoscalar PDAs also approach the pQCD asymptotics as the cut-off scale increases, although larger μ calculations are needed to confirm this point (cf. Fig. 9).
Also, connected to the scale dependence, it is interesting to note that approaches with a smaller cut-off, typically of the order of the constituent quark mass or even smaller (DSE results), lead naturally to distributions with lower moments, as shown in Table 4. CST and BLFQ, both have a larger cut-off, roughly of 2m (cf. Sect. 2.6), resulting in larger and comparable moments.
We also show the parton distribution functions in Fig. 10. Once again the results of the two approaches are consistent. In particular the first moments (displayed in Table 6) are very similar. On the other hand, the CST PDFs in Fig. 10 tend to display more pronounced wavy structures than the BLFQ PDFs for the η b (2S) and η c (2S).  Fig. 8 Comparison of CST and BLFQ longitudinal leading-twist distribution amplitude of J/ψ with respect to pQCD asymp prediction [62]; AdS/QCD (i) of Brodsky and de Téramond [44] and AdS/QCD (ii) [63]

Summary and conclusions
In this work we explore similarities and differences of two relativistic models of quarkonium. The successful results, reported in [2,3] were improved, as confirmed by a refined spectroscopy, and extended with a new calculation of decay constants. We find our results for observables to be in good agreement with each other and with several theoretical approaches including lattice methods. Beyond the level of observables, we used a map that allowed us to compare CST amplitudes with light-front wave functions and we observed a remarkable agreement between them. We used these lightfront wave functions to calculate parton distribution ampli-tudes and parton distribution functions. Our results again appear consistent with each other.
Noticeably, a first general conclusion is that in both approaches fixing the models solely to the mass spectra is sufficient to guarantee a reasonable overall description of the decay constants, whose precise description is known to be particularly challenging.
Since both approaches are established in Minkowski space, we could apply and test the BHL prescription, allowing us to perform not only a benchmark comparison between the two obtained LFWFs, but also to have direct access to quantities such as light-cone distributions. These are extracted here in a straightforward way, which is not possible in approaches relying on Euclidean formulations of quantum field theories.
The combined analysis of the two approaches enriches our understanding and the predictive power of each model alone, providing robustness tests to each other. For instance, while in BLFQ the angular momentum is not a good quantum number, the agreement with CST suggests that the used method to extract the angular momentum of each state is indeed reliable. Also, for the present 1CSE solutions, which do not possess a definite charge-conjugation parity, the comparison provides a way to judge the deviations in observables due to that violation.
All these conclusions point towards the need of comparisons between different approaches like the one in this paper. This way one obtains control on model dependencies and isolates method-independent features in the spectrum and production processes. This may be an advantage in the future, for instance when investigating the existence of exotic mesons.  Despite the formal difficulties, bottom up approaches shed light on difficult questions, such as how one may bridge different approaches and combine knowledge. This research has raised many questions, some of them in need of further investigation. One of them is certainly whether or not analogous linkages could be developed for lighter systems, which is planned for a future work.