Structure of Heavy Mesons in the Light-Front Quark Model

We investigate the structure of ground-state heavy mesons within the light-front quark model, utilizing wave functions derived from the Single Gaussian Ansatz (SGA) and the Gaussian Expansion Method (GEM). By performing a $\chi^2$ fit to static properties such as mass spectra and decay constants, we determine the model parameters for each approach. We then compare the impacts of both methods on the light-front wave functions and structural observables. Our analysis reveals significant differences in the distribution amplitudes (DAs) $\phi_{2;M}(x)$ near the endpoints, with GEM showing enhanced amplitudes and correct asymptotic behavior $\phi_{2;M}(x \to 1) \propto (1-x)$, consistent with perturbative QCD. This endpoint behavior is linked to the short-range (high-momentum) wave function governed by color Coulomb interaction and relativistic kinematics. GEM accurately reproduces a power-law damping $\psi_0(k \to \infty) \propto 1/k_\perp^2$, aligning with perturbative QCD predictions. Furthermore, the electromagnetic form factors of pseudoscalar mesons in the low-$Q^2$ region fall off faster with GEM than with SGA. Overall, while both methods adequately describe static properties, GEM provides a more accurate description of structural properties, being more sensitive to details and asymptotic behaviors.

Often, a single Gaussian Ansatz (SGA) [9][10][11] or power-law Ansatz [12,13,34] is employed for the LFWFs, whose parameters are fitted to the decay constants without considering the Hamiltonian [7,11,13].Alternatively, the LFWFs can be derived from the Bethe-Salpeter amplitude on the light front [21], where the regulator parameter is fitted to the data.Nevertheless, once the model parameters are well-tuned, the predictions of such models can be sufficiently consistent with the data.
In another approach, LFWFs can be obtained by directly diagonalizing the light-front Hamiltonian, as exemplified in basis light-front quantization (BLFQ) [25,26].Here, basis functions are constructed as a product of longitudinal and transverse components, resulting in WFs with cylindrical symmetry [35], rather than spherical symmetry.However, it should be noted that LFWFs have a nontrivial (x, k ⊥ ) dependence and in general cannot be separated [36].Due to this construction, the spherical symmetry of the WFs is not fully realized.This contrasts with LFWFs within the SGA, which inherently exhibit spherical symmetry by construction [37][38][39].
In a different approach, the effective Hamiltonian in the instant form is constructed, and the parameters of the Ansatz WFs are determined using the variational principle.Subsequently, the WFs are mapped into LFWFs [8,40].While this approach has proven successful in describing mass spectra and various observables [14,41,42], the SGA has limitations in describing some data and asymptotic behavior.For instance, within the SGA, the distribution amplitudes (DAs) near the endpoints are suppressed compared to lattice QCD data [19], and the fall-off of the calculated transition form factor is slower compared to BaBar data [43].It is, therefore, crucial to assess the limitations of the Ansatz by contrasting it with an approach that aims to closely resemble the eigenstate of the Hamiltonian.
One method to achieve that is the Gaussian Expansion Method (GEM) [44][45][46], which has proven its flexibility in many systems from atomic and nuclear physics [47].This method has also been applied to the non-relativistic quark model to obtain mass spectra and other observables, not only for mesons [48][49][50] but also for the baryons [51] and multiquark systems [52][53][54][55][56][57][58].The GEM relies on the construction of realistic WFs by utilizing Gaussian basis functions with multiple range parameters.This allows to approximate any shape of the WF and can be used to find the eigenstates of a given Hamiltonian.
In this article, we investigate the structure of groundstate heavy mesons within the LFQM, utilizing LFWFs obtained through both SGA and GEM.We focus on heavy mesons due to their suitability for probing the nonrelativistic limit of the quark model.To accomplish this, we conduct a χ 2 fitting procedure on static properties such as the mass spectra and the decay constants for each method, determining the parameters associated with the effective Hamiltonian.We then analyze the LFWFs and other related quantities such as DAs and electromagnetic (EM) form factors to understand the distinctions between the methods and provide comparisons with the experimental and lattice QCD data.By contrasting the results for both methods, we shed light on the structure of the heavy mesons.
Our investigation reveals that both GEM and SGA yield comparable accuracy in reproducing static properties.However, we identify differences between the two methods for the DAs and EM form factors, which are both more sensitive to the details of the WFs.In particular, we observe that the DAs near the endpoints for the GEM are pronounced compared to those in SGA, exhibiting the behavior ϕ(x → 1) ∝ (1 − x) [36].Also, the S-wave WF in the high-momentum region shows a power-law damping, ψ(k → ∞) ∝ 1/k 2 ⊥ , in line with predictions from perturbative QCD [59].We emphasize that the endpoint behaviors are linked to the short-range region, influenced by relativistic kinematics and Coulomb interaction.In the low-Q 2 region, the EM form factors of pseudoscalar mesons fall off faster for the GEM than those for the SGA.
The article is structured as follows.In Section II, we explain the basic components of LFQM and distinguish between SGA and GEM methods.Additionally, we outline the procedure for obtaining LFWFs and other properties.In Section III, we discuss the numerical results obtained through both methods.Finally, the article concludes with Section IV, summarizing our findings.

II. METHOD
In this section, we first describe the model Hamiltonian in the instant form.We then explain the differences between the two methods: SGA and GEM.The asymptotic behaviors of the WFs are also discussed.After that, we outline the procedure for constructing the LFWFs from the instant from WFs calculated by the two methods.Additionally, we present the observables considered in this work within the LFQM and the fitting procedures used to determine the model parameters.

A. Effective Hamiltonian
First of all, let us consider the relativistic Schrodinger equation where M q q and Ψ q q are the eigenvalue and eigenfunction for mesons made of quark and antiquark.The Hamiltonian in the instant form is given by with the usual nonrelativistic kinetic energy replaced by the relativistic one as with the quark (antiquark) mass m q (m q ) and the relative momentum k = (k z , k ⊥ ).However, in this case, we cannot factorize the c.m. motion and we can only work in the rest frame of mesons (P c.m. = 0).This approach is often referred to as the relativized quark model.In this work, we focus on the ground state spin 0 (pseudoscalar, P) and spin 1 (vector, V) heavy mesons, containing one or two heavy quarks, c or b.We adopt the QCD-motivated interquark potential V q q , which consists of the sum of the confining, color Coulomb, and hyperfine potentials, as given by where the term ⟨S q • S q ⟩ yields the values of 1/4 and −3/4 for the vector and pseudoscalar mesons, respectively.Although the tensor potential may contribute via the S-and D-wave mixing, its contribution is known to be weak in the quark model [60], and therefore we neglect it in the present work.
Overall the model has eight parameters, four of which are the quark masses (m q , m s , m c , m b ).The a and b are parameters for the confining potential, and α s is the strong running coupling, taken as a constant parameter.Here, we smear the spin-spin interaction with a Gaussian function as where Λ determines the strength of the smearing effect, and we introduce a phenomenological quark mass dependence as Λ = Λµ 1/2 q with the reduced mass 1/µ q = 1/m q + 1/m q .This accommodates the fact that the size gets smaller for heavier mesons.Not only the hyperfine splitting but also the decay constants are affected due to such a dependence [14].In this work, the model parameters in the Hamiltonian are determined through a χ 2 fit, which will be explained in Section II D.

B. SGA and GEM
In this study, we consider two approaches to solve Eq. ( 1) with the Hamiltonian in the instant form: (i) single Gaussian Ansatz (SGA) and (ii) Gaussian expansion method (GEM), whose resulting WFs are both mapped subsequently into the LFWFs.Here, we use the same form of the model Hamiltonian for both SGA and GEM but allow a different set of parameters for each.In the following, we explain the two methods in more detail.

Single Gaussian Ansatz
The single Gaussian Ansatz (SGA) relies on making an Ansatz for the meson WFs in the form of a single Gaussian function [8,40].The trial WF in position space is given by and the WF in the momentum space, obtained through the Fourier transformation, is given by where the S-wave spherical harmonic Y 00 = 1/ √ 4π is already included.It is important to note that r represents the relative coordinate between the quark and antiquark.While the Hamiltonian parameters are adjusted by fitting the mass spectra, the Gaussian parameter ν for each meson is determined through the variational principle ∂M q q /∂ν = 0.
Once the model parameters of the Hamiltonian are well-tuned, the predictions from this approach can reasonably agree with experimental data for a wide range of observables.However, it is important to note that a single Gaussian is an eigenfunction of the harmonic oscillator (HO) potential, not a general Hamiltonian.Because of this, the SGA does not accurately reflect the correct shape and asymptotic behavior of the WF for the Hamiltonian in Eq. (2).Therefore, the eigenstate for a given Hamiltonian is expected to deviate from a single Gaussian shape, even though the size of the WF after fitting can be similar.

Gaussian Expansion Method
To overcome the limitations of the SGA, we employ the Gaussian expansion method (GEM) [44][45][46] to solve the relativistic Schrödinger equation in Eq. (1).This method can be understood as a generalization of the SGA, as we increase the number of Gaussian functions until the necessary accuracy for approximating the solution of Eq. ( 1) is achieved.Despite its similarity, it is important to note that GEM is conceptually different from the use of an Ansatz, which does not aim to find the eigenstate of the Hamiltonian.
In this method, we expand the WF in terms of a set of Gaussian basis functions, ϕ G n , each with a different Gaussian parameter ν n as where c n represents the expansion coefficient.Following the notation in Ref. [44], the Gaussian basis function in position space is given by and the basis function in the momentum space, obtained through the Fourier transformation, is given by In the GEM, it is necessary to determine two sets of parameters: the coefficients c n and the Gaussian parameters ν n .To obtain the c n expansion coefficients satisfying ∂M q q /∂c n = 0, we solve the generalized eigenvalue problem where the elements of the Hamiltonian matrix are H q q,nm = ϕ G n Ĥ ϕ G m , and the elements of the overlap matrix are To obtain the Gaussian parameters ν n , we consider a geometric progression [44] (12) which reduces them to only two (ν 1 , ν nmax ) while keeping a high accuracy of the calculation.Those two parameters are optimized using the Optim.jlpackage [61] to satisfy We note that the basis functions are non-orthogonal S nm ̸ = δ nm and the states are normalized as C. LFWFs and Observables

Mass Spectra
In both SGA and GEM, we determine the mass spectra of the ground state heavy mesons by solving Eq. (1) in their rest frame.In this study, we calculate the mass spectra using instant form dynamics, which in line with the approach in the NRQM but incorporates relativistic kinetic energy [48][49][50].This method contrasts with BLFQ [25], where mass spectra are directly derived from LFWFs.

LFWFs
After diagonalizing the Hamiltonian, we need to map the obtained WFs into the LFWFs [14,27].For that we proceed in the following way 1.The position-space WF ψ(r) is first transformed into the momentum-space WF ψ(k).
2. The radial part of ψ(k) can be supplied to LFWFs Φ(x, k ⊥ ) via the k z → x variable transformation.Note that the Jacobian factor ∂k z /∂x is necessary to maintain the rotational symmetry.
3. The spin and orbital part of the LFWFs are obtained via the interaction-independent Melosh transformation [62].
Before proceeding further, we make some remarks about this mapping procedure.In our relativized quark model, we replace the kinetic energy term by the relativistic one.
For fully relativistic formulation based on the Dirac equation, the effects of the small components induce (spindependent) relativistic corrections that are important in the mass spectrum.In the non-relativistic formulation, instead, we include the spin-spin, spin-orbit and tensor interactions explicitly in the Hamiltonian.Although the relativistic effects are not fully included in the WF, our spectrum reproduces the observed one by adjusting the Hamiltonian parameters as we will see later.Then we map the WF to LFWF by the Melosh transformation [62], which is independent of the interaction and consistent with the Bakamjian-Thomas construction [63].It is not easy to quantify the ambiguity coming from this approximation in the final results, but our approach can be justified by comparing the results with data.Thus, this approach provides a well-controlled connection between the WF and LFWF.
In the following, we provide a more detailed demonstration of the mapping.The LFWFs are expressed in terms of the Lorentz invariant internal variables where P µ = (P + , P − , P ⊥ ) and p µ i denote the fourmomentum of the meson and the i-th constituent quark, respectively.Here we define the longitudinal momentum fraction x ≡ x q and the transverse momentum k ⊥ ≡ k ⊥q .
We perform the variable transformation (k z , k ⊥ ) → (x, k ⊥ ).In this case, x can be related with k z as where Therefore, k z can be written as with the so-called invariant meson mass The radial part of LFWFs is given by where for ψ(k) we take the WF from the SGA and GEM.The Jacobian factor is expressed as which takes into account the variable transformation.
The spin and orbital angular momentum part, R JJz λqλq , of LFWFs is obtained via the interaction-independent Melosh transformation [62] from the spin and orbital angular momentum part of the relativistic WF in the instant form assigned to the quantum number J P C .The R JJz λqλq has covariant forms as with M0 ≡ M 2 0 − (m q − m q ) 2 and the Dirac spinors of quark u(p q ) and antiquark v(p q ).The vertex Γ M for the pseudoscalar and vector meson (with M = V or P) is given by where the polarization vectors ϵ µ (J z ) = (ϵ + , ϵ − , ϵ ⊥ ) are defined by with ϵ ⊥ (±1) = (1, ±i) / √ 2. We note that the vertex Γ P can contain not only the pseudoscalar coupling, but also the pseudovector coupling.Although this pseudovector is usually not important in the low-energy regime, it can affect the asymptotic behavior of form factors as Q 2 → ∞ [64].
The explicit forms of the spin and orbital WFs for the pseudoscalar and vector mesons are given by and respectively, where k The LFWF of the ground state heavy meson in momentum space is therefore given by and normalized as λq,λq Furthermore, the LFWFs are often constructed as products of Ψ L (x) and Ψ T (k ⊥ ).However, this construction often lead to a breaking in spherical symmetry of the WF [25].In our case, the LFWFs are transformed from Gaussian basis functions, thereby maintaining the spherical symmetry of the LFWFs.One evidence is that would yield the same results when using the LFWFs in Eq. ( 29).This consistency can be broken if LFWFs are constructed as Ψ = Ψ L (x)Ψ T (k ⊥ ).Interested readers may refer to the previous study [37] for more details.

Decay Constants
Now we provide formulae for the decay constants of the pseudoscalar meson f P and vector meson f with longitudinal and transverse polarizations.They are defined by where ϵ µ (J z ) and M V represent the polarization vector and the mass of the vector meson, respectively.
The explicit form of the decay constants computed in the LFQM using the plus current (µ = +) is given by [14] where Φ(x, k ⊥ ) is radial part of LFWFs and the operators O M read The decay constants can also be calculated using the minus (µ = −) and the transverse (µ =⊥) components of currents.The equivalence of the decay constants with various current components and polarizations has been demonstrated within this LFQM in a previous work [37].

Twist-2 Distribution Amplitudes
Next, we focus on the twist-2 DAs which give dominant contributions in the hard exclusive processes [65].The DAs are derived from matrix elements that connect free space and meson states with light-like separation, i.e., z 2 = 0. To establish a link between the DAs and the LFWFs, we apply the condition of the equal light-front time to the light-like vector z µ with z + = z ⊥ = 0.The twist-2 DAs for pseudoscalar mesons ϕ 2,P computed by choosing the plus current are given by [66] ⟨0| q(z)γ + γ 5 q(−z) where ξ = 2x − 1.The twist-2 DAs for the vector mesons with longitudinal ϕ ∥ 2,V and transverse polarizations ϕ ⊥ 2,V are computed as [67] ⟨0| q(z)γ + q(−z) ⟨0| q(z)σ ⊥+ q(−z) respectively.
In the LFQM, the ϕ M (x) can be obtained by the transverse momentum integration of the LFWF as [9] which are normalized as In the model calculation, it is difficult to precisely determine the scale µ, but it is associated with the modeling interaction.The scale dependence of DAs can be obtained by using the QCD evolution equation [65].where the parameters are fixed from using the SGA (n = 1).The figures also show the results when we expand the basis function up to n = 10 while keeping the parameters the same.Evidently, using the same parameters for n = 10 would result in a poor prediction, especially in the decay constants, which deviate significantly up to 30%.Thus, it is necessary to fit the parameters independently when using the GEM (n = 10).
Moreover, we provide the six-lowest ξ-moments of the DAs, defined as which can be compared with other models quantitatively.These ξ moments can be related to the Gegenbauer moments a n (µ) [9].Although features of DAs are reflected in the first few ξ moments, accurate parameterization to DAs of heavy-light mesons may require more moments due to a pronounced asymmetry [68].

Electromagnetic Form Factors
Furthermore, we want to test the two methods on the transfer momentum (Q 2 )-dependent quantities.For that reason, we include the EM form factor of pseudoscalar mesons where This EM form factor is sensitive to the quark masses and provides further tests on the details of the LFWFs.A more comprehensive study involving elastic and transition form factors is left for future work.
In LFQM, the EM form factor of pseudoscalar mesons is computed by using the Drell-Yan-West frame (q ⊥ and obtained by using the plus component of current J + .The explicit expression is given by [40] F (Q 2 ) = e q I + (Q 2 , m q , m q ) + e q I + (Q 2 , m q , m q ), (44) where the e q (e q ) is the electric charge of the quark (antiquark).The contribution of the quark and antiquark is calculated by where The form factor computed from the transverse current would also give the same results in this LFQM [69].Note that the EM form factor is normalized as F (Q 2 = 0) = e q + e q , and the corresponding charge radius is computed as In the case of heavy quarkonia such as η c and η b , we only consider the contribution from one of the quarks, as otherwise the EM form factor vanishes because the contributions from the quark and antiquark cancel each other out.However, the EM form factors for neutral mesons such as D 0 , B 0 , and B 0 s do not vanish due to different quark flavors.

D. Fitting Procedures
To determine the model parameters of the Hamiltonian in Eq. ( 1), we perform a χ 2 fit, as defined by where the minimization is performed using the Optim.jlpackage [61].The data set O obs i in the fit includes experimental data of both mass spectra and decay constants.
Fitted mass spectra (upper panel) and decay constants (lower panel) of ground-state heavy mesons using SGA (red) and GEM (blue).In addition, we display data from experiments (light green) [70] and lattice QCD calculations (light blue) [13,[71][72][73][74].We incorporated both observables in our fitting process.TABLE I. Model parameters for SGA and GEM which are independently fitted to the data by minimizing χ 2 .The parameters include the constituent quark masses (mq, ms, mc, m b ) in unit of GeV, confinement potential parameters a [GeV] and b [GeV 2 ], the strong coupling αs which is dimensionless and the smearing parameter Λ [GeV 1/2 ] with a quark mass dependence Λ = Λµ Here we have N data = 32 data points and Npar = 8 free parameters resulting in N dof = N data − Npar = 24.We obtain χ 2 /N dof is 0.95 and 0.46 for GEM and SGA, respectively.Additionally, we incorporate data from lattice QCD simulations for the decay constants [71][72][73][74].When minimizing χ 2 to fit the data, it is also necessary to optimize the Gaussian parameter for each iteration in the χ 2 fit to ensure the minimum of the energy.This is essential to satisfy the variational principle, as otherwise, the results are not valid due to not reaching the energy minimum.

Model
Here, we adjust the model error σ mod i = σ err O obs i , where σ err is the percentage error of the O obs i , to get χ 2 /N dof ≈ 1.The N dof represents the number of degrees of freedom obtained by subtracting the number of free parameters (N par = 8) from the number of data points (N data = 32).We set σ c err = 2%, σ b err = 1%, and σ f err = 5%, for the model error of the mass of charmed mesons (qc, sc, cc), bottom mesons (q b, s b, c b, b b), and decay constant for all of the heavy mesons, respectively.The inclusion of this additional σ mod in the χ 2 ensures an unbiased fit, especially when the data precision varies significantly [76].Note that σ mod becomes dominant when σ expt is negligibly small.Furthermore, we use a relative σ mod because mass spectra and decay constants differ by up to an order of magnitude.
It should be noted that the fitted parameter sets for SGA and GEM are significantly different.Indeed, Fig. 1 shows that considering a fixed parameter set for both methods, the decay constants may differ by up to 30%.This shows again that the two approaches are conceptually different.Consequently, we need to perform a fit to the data independently for both methods.In GEM, using 10 basis functions yields an accuracy of about 4 digits.

III. RESULTS AND DISCUSSION
In this section, we first present the obtained model parameters for SGA and GEM from fitting to static properties such mass spectra and decay constants.Subsequently, we continue to discuss the difference between the two methods for LFWFs and structural observables.

A. Model Parameters
The results for the fitted model parameters for both SGA and GEM are presented in Table I.Although the form of the Hamiltonian is the same for both methods, the model parameters obtained from the fit show differences.The obtained mass spectra and decay constants for ground-state heavy mesons, both included in our fitting process, are shown in Fig. 2.Although χ 2 /N dof for the SGA is smaller, it is less then one for both methods with σ mod ranging from 1-5% as explained previously.This shows that the SGA is simple yet powerful while the GEM can provide a more accurate WF but requires more effort to obtain a more refined model Hamiltonian to improve agreement with the data.
Initially, we attempted to fit with unconstrained parameters for SGA and GEM.In this setup, we found reasonable agreement with the data, but the obtained m q was rather small, around 60 MeV, which is smaller than the typical constituent quark mass m q = 200-350 MeV.Because of that, we tested other observables such as the EM radius with this parameter set, and the obtained radius for D + was around 0.480 fm 2 , much larger than the lattice QCD data of around 0.150 fm 2 discussed in Sec.III F. Furthermore, m q was more tightly constrained with the observables for light mesons, which were not included in the present work.Due to these observations, we set the lower bound for m q = 220 MeV to the commonly used value in Refs.[9,60] for our analysis.
The spatial dependence of the potentials for the SGA and GEM are shown in Fig. 3 and compared with the Godfrey-Isgur (GI) model [60] and Isgur-Scora-Grinstein-Wise (ISGW) with α s = 0.3 [75].While both potentials in this work are comparable with those of the literature, one can see that the potential for the GEM is more enhanced at a short distance and suppressed at a long distance as compared to those for the SGA.It is worth noting that the outcome for GEM is notably influenced by 0.0 0.5 1.0 Spatial dependence of the confinement and color Coulomb potentials for the SGA and GEM using the model parameters in Table I.The potentials are comparable with those of the GI [60] and ISGW potentials [75].
variations in the model Hamiltonian, leading to distinct shape, as illustrated in Fig. 3.In this case, the enhancement of the WF at the origin for the GEM as shown in Fig. 1, which is plausibly due to the use of relativistic kinetic energy, leads to the best fit with different potential parameters.

B. Wave functions and their asymptotic behaviors
In Fig. 4 we show as an example the WF (upper panel) and density (lower panel) of η c for both SGA and GEM.Unlike those obtained with the same parameters as shown in Fig. 1, the WFs of both methods are now much closer to each other.Nevertheless, we can identify some distinct differences, in particular, the WF for the GEM is more enhanced at the origin and extends to larger distances.If we replace the kinetic energy with the nonrelativistic one while keeping the parameters the same as those for the GEM, we obtain that the WF denoted as the GEM-NR is rather suppressed at the origin as shown in Fig. 4.This is mainly because the relativistic kinetic energy provides weaker repulsion near the origin than the nonrelativistic one.
Moreover, we plot the LFWFs (Ψ 00 ↑↓ − Ψ 00 ↓↑ )/2 of the η c The GEM and GEM-NR employ a relativistic and nonrelativistic kinetic energy, respectively.As compared to the one for the SGA, the WF for the GEM is more enhanced at short and long distances, but it is more suppressed at an intermediate distance.The WF near the origin for the GEM-NR is more suppressed than the one for the GEM. and D meson in the upper and lower panels of Fig. 5, respectively, exemplifying the cases of equal and unequal mass of constituents.The peak for the η c LFWF, which corresponds to (k z , k ⊥ ) = 0, is at x = 1/2 but the peak for the D meson LFWF is at x = m q /(m q + m q ) with x carried by the light quark.Note that the endpoints x = 0 and x = 1 correspond to k z = −∞ and k z = ∞, respectively.Evidently, the LFWF for the GEM extends more to both endpoints of x and the larger k ⊥ region, which we discuss in more detail in Sec.III.D.
Since the GEM accurately approximates the eigenstate of the Hamiltonian, while the SGA only captures the size of the WF, it is important to check the asymptotic behavior of the WF in two different limits: long distance (r → ∞) and short distance (r → 0).Of special interest is the short-distance behavior, which corresponds to the high-momentum behavior (k → ∞), which is dictated by perturbative QCD.
In the long-distance limit r → ∞, the Schrödinger equation for the S-wave with the power-law confinement reduces to with the reduced quark mass µ, the confinement parameter b, and exponent p.If we assume a decaying WF ψ(r) ∝ exp(−Br n ), we can obtain the asymptotic parameters n = (p + 2)/2 and B = 2µb/n 2 .While a HO confinement (br 2 ) results in a Gaussian asymptote, the linear confinement (br) considered here leads to which decays more slowly compared to the HO case, as expected.Note that the Airy function is the eigenfunction of the linear potential [77], and we reproduced here its asymptotic limit.
In the short distance, the dominant interaction is the color Coulomb potential, stemming from the one-gluon exchange and the relativistic effect becomes important as the relative momentum becomes large.As r → 0, the Schrödinger equation is reduced to where we consider the equal-mass case m = m q = m q for simplicity.Near the origin, the WF has a power-law behavior ψ(r → 0) → r n , which is also observed when using the Dirac Hamiltonian [78].Such behavior is well approximated by the GEM, as shown in Fig. 4, but differs significantly in the SGA.Additionally, this asymptotic WF is divergent, but the truncation of the number of basis functions in the GEM regulates it.
To explicitly demonstrate a power-law behavior in the model calculation, it is more effective to examine the WF in momentum space, as this reveals the behavior more clearly.The momentum-space WF is obtained via Fourier transformation as given by where we have used the plane-wave expansion for the S wave.We obtain Figure 6 shows the WF integrated over k z , for η c and as k → ∞.From the inset of the upper panel, it is evident that the GEM gives a linear dependence up to about k ⊥ = 50 GeV, which is eventually broken at higher momenta.This is because the GEM uses Gaussian basis functions, and the accuracy is limited by the narrowest (widest) basis function in position (momentum) space, and the number of basis functions.Our fit of the asymptotic WF between k ⊥ = 5-20 GeV on a loglog scale reveals a linear relationship with a slope of approximately −2.This corresponds to a damping factor of 1/k 2 ⊥ , which also implies that the ψ(r → 0) ∝ 1/r.This finding is consistent with the predictions of perturbative QCD [59] and the Bethe-Salpeter method [36].

C. Mass Spectra
The mass spectra for the ground state of heavy mesons obtained by the SGA and GEM are presented in Table II.Both methods yield similar results which are consistent with the experimental data of the Particle Data Group [70].Since the mass spectra originate from the average of short-and long-distance effects of the WFs, a similar result from both methods is expected.
Here we use different σ mod i for mesons containing the charm (σ c err = 2%) and bottom quark (σ b err = 1%) to obtain a better fit with the data for all heavy mesons.If we employ the same σ err = 2%, the results for heavier mesons, such as bottomonia, will be less accurate.On the other hand, if we use absolute values such as σ mod = 50 MeV, the results for the heavy-light meson become less accurate.The use for the relative error in σ mod becomes important if we include light mesons in the fit since their masses are much smaller.Furthermore, in this work, we analyze not only the mass spectra but also other observables.In this case, the results should have reasonable agreement with the data for other observables as well.If we fit only the mass spectra, the prediction for decay constants can be bad.Therefore, the mass spectra provide a necessary, but not sufficient condition to constrain the WF.

D. Decay Constants
The decay constants obtained for both SGA and GEM are tabulated in Table III, and compared with experimental data, extracted from the leptonic and weak decays [70], and lattice QCD data [71][72][73][74].Note that there are some discrepancies between both lattice QCD and experimental data, because of which we employ σ f err = 5% in the fit.Since there is no experimental data for f ⊥ V , we only compare our results with other theoretical model [13].We find that the results for both methods have reasonable agreement with the data as shown in Table III.Although in the nonrelativistic limit the decay constant is related to the WF at the origin via the famous Van Royen-Weisskopf formula [79], the results for both methods are comparable once the range parameters of the WFs are fitted to the data.Therefore, it is of great interest to analyze momentum-dependent quantities, instead of constant observables, to unveil the difference between the two methods.
While f P < f ∥ V is always held if the spin-spin interaction is treated perturbatively [41], f P can be larger than f ∥ V when the spin-spin interaction is treated nonperturbatively and the smearing parameter Λ plays a crucial role in determining the hierarchy as discussed previously [14].For instance, is observed despite the quark mass dependence added  [70] as well as lattice QCD data [71][72][73][74].We also compare the results for f ⊥ V by Ref. [13]. to the Λ parameter.At the moment, lattice QCD and experimental data for f ∥ Υ have some discrepancy.More precise data is therefore desirable to resolve the hierarchy which is useful to further constrain the WF.

SGA
Furthermore, from Eq. ( 34), we see that if the radial WFs of pseudoscalar and vector mesons are the same [41].Here we find that f P > 2f ⊥ V −f ∥ V when the spin-spin interaction is treated nonperturbatively and it applies to both SGA and GEM.

E. Twist-2 Distribution Amplitudes
In Fig. 7, we show the leading-twist DAs of pseudoscalar mesons for both SGA and GEM.Overall, the results are comparable, but they exhibit different behaviors near both endpoints.For heavy-light mesons, the difference is more evident near x = 1.For instance in the case of the B meson, the DAs in the two methods look quite different, and these features are commonly observed in other mesons.
The ϕ 2;B (x) for the GEM is more suppressed in the region of 0.1 < x < 0.3 and more enhanced in the region of x > 0.3.While the ϕ 2;B (x) in GEM is more extended to the x > 0.5 region, the ϕ 2;B (x) for the SGA is concentrated in the x < 0.5 region.The endpoint behavior at x = 1 is related to the high-momentum part of the WF and accordingly a short-distance part of the WF.For the SGA, this suppression can be directly inferred from the WF which is rather suppressed at a short distance as shown in Fig. 4. On the other hand, the WFs for the GEM are enhanced at a short distance.
We also plot the asymptotic behavior of the ϕ 2;ηc (x) in Fig. 8.It is evident that the DA near x = 1 for the GEM shows a linear dependence as ϕ 2;ηc (x → 1) = 2.13(0.99− x), while it behaves differently in the case of SGA.The linear behavior is needed to yield a meson PDF of (1 − x) 2 near x = 1, as predicted by perturbative QCD [80].It is worth noting that we fit the DA between 0.9 ≤ x ≤ 0.95.as they eventually start to deviate from the linear form when it is much closer to x = 1, similar to Fig. 6.In the BLFQ approach [25], this asymptotic WF is used as a basis function such as Ψ L (x) ∝ x α (1−x) β .This produces the asymptotic behavior of DAs by construction, but it may also break the spherical symmetry.Furthermore, a previous study [19] has shown that the DAs using SGA could not reproduce the enhancement near the endpoints seen in lattice QCD data [81].This clearly shows the limitations of the SGA due to the fixed form of the WF.These observations suggest the use of the eigenfunction of Hamiltonian instead of a single Gaussian Ansatz.In contrast, the superposition of Gaussian basis functions with different range parameters in GEM can produce a general WF with richer features.Furthermore, it is worth noting that the enhancement near the endpoints of DAs can alternatively be obtained using the power-law Ansatz [13], which is an asymptotic WF in the high momentum region.However, this power-law Ansatz has also some limitations in its shape and has other problems such as convergence issues.
In Table IV, we tabulate the six lowest ξ-moments for the pseudoscalar DAs computed within the SGA and GEM.We find that the ξ moments for the GEM are generally a bit larger than those from the SGA.Moreover, the ξ moments for mesons containing a bottom quark have larger deviations between both models.Such a large difference is reflected in the DAs as shown in Fig. 7.We also compare the ξ moments with the SGA and power-law ansatz computed in Ref. [13].For example, our computed odd ξ moments for the B meson, ⟨ξ 1 ⟩ = −0.629[−0.508]with SGA [GEM], are comparable to those in Ref. [13], ⟨ξ 1 ⟩ = −0.617[−0.531]with SGA [Power-law].This shows that our GEM results for the ξ moment are more in line with those for the power-law ansatz.This can be understood as they have more enhancements near x = 1, which make them less asymmetrical with respect to x = 1/2 and yield smaller absolute values of odd ξ moments.
Since we now know the qualitative difference in the predictions of both methods, it is also interesting to compare the DAs for pseudoscalar and vector mesons obtained by the GEM.The comparisons are provided in Fig. 9 where we plot the ϕ 2;P (x), ϕ ∥ 2;V (x), and ϕ ⊥ 2;V (x).The ϕ 2;P (x)  are similar to the ϕ ∥,⊥ 2;V (x) for bottom mesons, partly due to the heavy-quark symmetry, but there is some notice-able difference for the case of charm mesons.One can see that ϕ ∥,⊥ 2;V (x) are generally comparable with each other, but have higher peaks as compared to ϕ 2;P (x).When we use Gaussian basis functions with different range parameters (GEM), the resulting DAs always have a single peak.This is rather different than those obtained by expanding into the HO basis function where the mixture of higher nS can lead to oscillatory DAs [42] for the ground states such that the mixture is restricted to be very small [41].
For completeness, we also provide the six lowest ξ moments for vector meson DAs in Table V.The ξ moments for the bottom vector and pseudoscalar mesons are rather similar, following their DAs as shown in Fig. 9.For the charm mesons, the odd (even) ξ moments for ϕ ∥,⊥ 2,V are smaller (larger) than those of ϕ 2,P where the biggest difference is from the lowest odd ξ 1 and even ξ 2 moments.While the even ξ moments of ϕ ⊥ 2,V are generally smaller than those of ϕ ∥ 2,V .For the charmonia and bottomonia, the even (odd) ξ moments of ϕ ⊥ 2,V are larger (smaller) than those of ϕ ∥ 2,V for the heavy-light mesons.FIG. 9. Twist-2 DAs of pseudoscalar ϕ2;P (x) and vector heavy mesons with longitudinal ϕ ∥ 2;V (x) and transverse ϕ ⊥ 2;V (x) polarization.Although we include the spin-spin term nonperturbatively, the differences in the DAs of pseudoscalar and vector mesons are rather small.Only the mesons with charm quark show some visible difference in the DAs.TABLE V. Six lowest ξ moments of ϕ ∥ 2;V (x) and ϕ ⊥ 2;V (x) obtained with the GEM, where ξ = 2x − 1.The results of the EM form factors of the pseudoscalar heavy mesons for both SGA and GEM are presented in Fig. 10, together with the lattice QCD [82][83][84][85].We find that the results for D + and D + s in both methods are comparable and consistent with the lattice QCD [82,84].For the form factor of the η c , both methods seem consistent with Ref. [83], but only the result for the GEM can reproduce the lattice QCD (orange dashed line) of Dudek et al. [85].Note that we consider only the quark contribution for the case of η c and η b , otherwise, the form factor vanishes.
We see that the fall-off of EM form factors in low-Q 2 region for the GEM are generally faster than those for the SGA.Even so, the difference between them may depend on the fit as the form factor depends on the quark masses and the potential parameters, where typically the smaller quark masses produce a faster fall-off of the form factor.
The mean squared of the charge radii from the form factors are shown in Table VI, indicating that the obtained radii for the GEM are generally larger compared to those for the SGA.This can be understood from the density of the WF plotted in Fig. 4. Since the WF for the GEM extends more to long distances, the expected radius of r 2 EM is larger than those for the SGA.Nevertheless, the obtained radii in both models are consistent with current lattice QCD [82][83][84][85].In particular, for the B + meson, the results from the two methods show a large deviation, i.e., r 2 EM = 0.297(0.459)fm 2 for the SGA (GEM), respectively.Therefore, more lattice QCD data on this observable is necessary to further constrain the models.
It is also crucial to examine the fall-off of the form factor in the high-Q 2 region, as it is dictated by perturbative QCD.As a demonstration, in Fig. 11, our calculations show that Q 2 F D + (Q 2 ) decreases with increasing Q 2 rather than remaining constant.Although there are logarithmic corrections ln Q 2 due to the running of the strong coupling constant α s (Q 2 ) [87] that affect the Q 2 dependence, there could be other contributing factors.Then, we fitted the form factor Q 2 F D + (Q 2 ) using GEM for Q 2 values between 60 and 100 GeV 2 on a log-log scale as illustrated in Fig. 11 B 0 s FIG. 10.Computed EM form factor for pseudoscalar heavy mesons compared with available lattice QCD data [82][83][84][85].Note that only the quark contribution is considered for ηc and η b , as otherwise the form factors vanish.
show that F (Q 2 ) ∝ 1/(Q 2 ) n with n = 1.63, as indicated in the inset.When the form factor is adjusted by multiplying (dividing) by ln Q 2 , the exponent becomes n = 1.40 (n = 1.86), deviating from the perturbative QCD prediction of n = 1.In a previous study of the pion form factor [64], it was suggested that including the pseudovector component in the meson vertex using the Bethe-Salpeter approach, which is not considered in this work, could yield the correct asymptotic behavior.Furthermore, some nonperturbative effects may also give different Q 2 dependence in the form factor at high-Q 2 region [88].Further investigation into this matter is essential to fully understand the underlying mechanisms.

IV. CONCLUSION AND OUTLOOK
We have investigated the structure of heavy mesons using the single Gaussian Ansatz (SGA) and Gaussian expansion method (GEM) within the Light-Front Quark Model (LFQM).To accomplish this, we have concentrated our efforts on the ground state of heavy mesons and investigate not only their static properties, but also their structural properties.To determine the model parameters, we have performed a simultaneous χ 2 fit to the static properties such as mass spectra and decay constants and examined the difference in the predictions of both methods especially in the structural properties such as DAs and EM form factor.We found that both methods yield similar static properties such as mass spectra and decay constants, and given the model uncertainty, they exhibit reasonable agreement with experimental and lattice QCD data.However, they show differences in the LFWFs and structural properties.In particular, the asymptotic behaviors of the WFs ψ 0 (k → ∞) ∝ 1/k 2 ⊥ and DAs ϕ(x → 1) ∝ (1 − x) are correctly reproduced by the GEM, while they are not in the case of the SGA.These behaviors in the high-momentum region are governed by relativistic kine-matics and Coulombic one-gluon exchange, which produce a power-law fall-off of the WF ψ(r → 0) ∝ 1/r.Furthermore, the fall-off of the EM form factor in the low-Q 2 is faster for the GEM, giving better agreement when compared to the lattice QCD [85], especially for the η c meson.
For future work, several directions can be explored: first, it is crucial to investigate the form of the model Hamiltonian since the WF obtained in the GEM model is sensitive to the Hamiltonian, unlike in the SGA.Furthermore, expanding our calculations to include excited states and light mesons is of great importance in testing the applicability of GEM.Finally, our model can be tested on other form factors and hadron distributions.

FIG. 1 .
FIG. 1.(a) WF ϕη c (r), (b) density of WF r 2 |ϕη c (r)| 2 , (c) mass, and (d) decay constant of ηc(1S)where the parameters are fixed from using the SGA (n = 1).The figures also show the results when we expand the basis function up to n = 10 while keeping the parameters the same.Evidently, using the same parameters for n = 10 would result in a poor prediction, especially in the decay constants, which deviate significantly up to 30%.Thus, it is necessary to fit the parameters independently when using the GEM (n = 10).
D s D * s c J/ B B * B s B * s B c B *

2 FIG. 4 .
FIG.4.WF (upper panel) and the density (lower panel) of the WFs of ηc(1S) for both SGA, GEM, and GEM-NR.The GEM and GEM-NR employ a relativistic and nonrelativistic kinetic energy, respectively.As compared to the one for the SGA, the WF for the GEM is more enhanced at short and long distances, but it is more suppressed at an intermediate distance.The WF near the origin for the GEM-NR is more suppressed than the one for the GEM.

FIG. 5 .
FIG. 5. LFWF (Ψ 00↑↓ − Ψ 00 ↓↑ )/2 of the ηc (upper panel) and of the D meson (lower panel) for both SGA (left half) and GEM (right half).The LFWFs for the GEM are more extended to both endpoints of x but suppressed at the intermediate x as compared to those for the SGA.The dashed line represents kz = 0 obtained from Eq. (18).

FIG. 6 .
FIG. 6. (Upper panel) The k ⊥ dependence of the WF ψ(k ⊥ ) where the kz component has been integrated.The inset shows a k ⊥ dependence up to 100 GeV in a logarithmic scale.(Lower panel) a linear fit to a WF ψ(k ⊥ ) on a log-log scale.It implies that ψ(k → ∞) ∝ 1/k 2 ⊥ .

FIG. 7 .
FIG.7.Twist-2 DAs ϕ2;P (x) of pseudoscalar heavy mesons with various quark flavor contents for both SGA (red) and GEM (blue).Compared to the results for the SGA, we find the DAs for the GEM are more enhanced near both endpoints.

FIG. 11 .
FIG.11.EM form factor of the D + meson in the wide range of transfer momentum Q 2 .We perform a linear fit to the form factor F (Q 2 ) between Q 2 = 60-100 GeV 2 on a log-log scale.

TABLE II .
[70]rical results of mass spectra [MeV] of the ground state of heavy mesons for both SGA and GEM compared with the experimental data in Ref.[70].

TABLE III .
Numerical results of decay constants [MeV] of ground-state heavy mesons for both SGA and GEM, and compared with experimental

TABLE IV .
Six lowest ξ-moments of ϕ2;P (x) for both SGA and GEM, where ξ = 2x − 1.The ξ moments for the GEM are generally larger than those for the SGA, except for the even ξ moments of B and Bs mesons.