Mixing effects on 1S and 2S state heavy mesons in the light-front quark model

The mass spectra and wave functions of both $1S$ and $2S$ state heavy pseudoscalar ($P$) and vector ($V$) mesons are analyzed within the light-front quark model. Important empirical constraints employed in our analysis of the mass spectra and wave functions are the experimental mass-gap relation, $\Delta M_P>\Delta M_V$, where $\Delta M_{P(V)}=M^{2S}_{P(V)}-M^{1S}_{P(V)}$ and the hierarchy of the decay constants, $f_{1S}>f_{2S}$, between $1S$ and $2S$ meson states. We maintain the orthogonality of the trial wave functions of the $1S$ and $2S$ states in our variational calculation of the Hamiltonian with the Coulomb plus confining potentials and treat the hyperfine interaction perturbatively for the heavy-heavy and heavy-light $P$ and $V$ mesons due to the nature of the heavy quark symmetry. Realizing that the empirical constraints cannot be satisfied without mixing of the $1S$ and $2S$ states, we find the lower bound of the mixing angle $\theta$ between $1S$ and $2S$ states as $\theta_c = {\rm cot}^{-1}(2\sqrt{6})/2\simeq 6^{\circ}$ and obtain the optimum value of the mixing angle around $12^\circ$ to cover both the charm and bottom flavors of the heavy quark. The mixing effects are found to be more significant to the $2S$ state mesons than to the $1S$ state mesons. The properties of $1S$ and $2S$ state mesons including the mass spectra, decay constants, twist-2 distribution amplitudes, and electromagnetic form factors are computed. Our results are found to be in a good agreement with the available data and lattice simulations. In particular, the $2S$ state pseudoscalar $D_s$ meson is predicted to have a mass of $2600$ MeV, which is very close to the mass of the newly discovered $D_{s0}(2590)^+$ meson by the LHCb Collaboration. This supports the interpretation of the observed state as a radial excitation of the $D^+_s$ meson.


I. INTRODUCTION
Quantum chromodynamics (QCD) is a unique theory of strong interactions with its non-perturbative nature in the lowenergy regime and its asymptotically free nature in the highenergy regime. Building the effective degrees of freedom that describe the strongly interacting system in the low-energy regime is one of the crucial issues in understanding the link between the first-principle QCD and the constituent quark model (CQM) that has proven to provide successful and intuitive descriptions of hadrons. In particular, the light-front dynamics (LFD) is found to provide an effective way to handle the relativistic effects thanks to its distinguished features of the rational energy-momentum dispersion relation. It carries the maximum number (seven) of the kinetic (or interactionindependent) generators rendering the less effort in dynamics to get the QCD solutions that reflect the full Poincaré symmetries [1,2]. Effectively, the light-front quark model (LFQM) based on the LFD turns out to be one of the most successful hadronic models in describing various properties of hadrons.
While the LFQM analyses have been quite successful in describing the properties of ground state mesons [3][4][5][6][7][8][9][10][11][12][13][14][15], the structures and properties of the excited hadron states are yet to be understood in LFQM more extensively as their nature is still veiled and not well explored compared to the ground states. One of the most challenging problems in the quest of the excited states is to clarify whether the observed state belongs to the standard quark-antiquark excitation or an exotic state. For instance, the newly observed 0 (2590) + meson with a mass of 2591 ± 6 ± 7 MeV by the LHCb Collaboration [16] has been proposed as a radial excitation of the + meson. However, the observed mass is quite smaller than the available CQM predictions. For example, the relativized quark model of Ref. [17] predicts 2680 MeV and the relativistic quark model based on the quasipotential approach predicts 2688 MeV [18,19]. In recent works, therefore, some nonstandard quark-antiquark behaviors attributed to this resonance have been discussed [20,21].
In particular, the radially excited states of hadrons are important in understanding the strong interactions as they give information complementary to the orbitally excited states. They have been observed in light and heavy quark sectors of hadrons although some of them are yet to be confirmed according to the Particle Data Group (PDG) [22]. The most well-known example in baryon spectrum would be the Roper resonance [23]. Such states with multistrangeness were recently discussed in Ref. [24], and the excited states in meson spectrum were discussed and summarized, for example, in Ref. [25]. Most notable empirical hierarchy appears in the radially excited 2 state as well as the ground 1 state of heavy pseudoscalar ( ) and vector ( ) mesons. Namely, the two important constraints that we notice from the empirical hierarchy are (i) the experimental mass gap relation Δ > Δ , where Δ ( ) = 2 ( ) − 1 ( ) , and (ii) the hierarchy of the decay constants 1 > 2 . While these constraints on the mass spectra and the decay constants apply both for the light and arXiv:2205.04075v1 [hep-ph] 9 May 2022 heavy meson sectors, the difference between and becomes much larger in the light meson sectors. The reason for the smaller difference between and in the heavy-light system may well be attributed to the heavy quark symmetry [26][27][28][29] as well as the perturbative nature of the hyperfine interaction in heavy quark systems. This may be contrasted with the chiral symmetry reflected in the light meson sectors, which deserves a separate analysis and discussion. Due to these significant differences in the underlying symmetry between the light and heavy meson sectors, we apply the two constraints here only for the heavy meson sectors and discuss the heavy-heavy and heavy-light and mesons in the present work. Effectively, we analyze the mass spectra and wave functions of the radially excited 2 state and the ground 1 state of heavy and meson sectors within the framework of the LFQM and discuss various properties of these mesons.
In the previous LFQM analyses of the mass spectra and decay constants of the 1 state mesons performed by two of us with the QCD-motivated effective Hamiltonian [10][11][12][13][14][15], the trial wave functions were chosen as either the pure harmonic oscillator (HO) wave function 1 [10][11][12][13] or an expansion in the HO basis functions, i.e., Φ = max =1 with max = 2 [14] or 3 [15]. Through the analyses of the 1 state mesons, it was shown that the physical observables are not much sensitive to the number of HO bases used in the trial wave functions, Φ = max =1 , once the optimum values of the model parameters are fitted. For the combined analysis of (1 , 2 ) state heavy mesons in the present work, we take into account the two important constraints in obtaining the optimum values of our model parameters. We note that tends to be smaller as gets larger since the decay constant of a hadron is proportional to its radial wave function at the origin, ( = 0). The available experimental data have also confirmed this tendency. In the literature, however, some difficulties in the combined analysis of the ground and radially excited states have been observed. For instance, in the LFQM analysis of (1 , 2 ) state heavy Υ(¯) systems, it was found that using the HO wave functions ( = 1, 2) leads to the reverse order problem of 1 < 2 [7]. In order to resolve this problem and to obtain the correct hierarchy of 1 > 2 for the heavy quarkonium system using two pure ( 1 , 2 ) wave functions, the authors of Refs. [30,31] had to choose different HO model parameters for different states, which breaks the orthogonality condition between the two wave functions 1 and 2 . A similar problem, namely, the breakdown of orthogonality condition between the resultant 1 and 2 state light-front (LF) wave functions, appears in the analysis of heavy quarkonium system performed with the basis light-front quantization (BLFQ) approach in a holographic basis [32].
The main purpose of the present work is thus to extend the previous LFQM analyses including both 1 and 2 state and heavy meson sectors to remedy the difficulties in the combined analysis using the trial wave functions for 1 and 2 states as mixtures of the two HO wave functions 1 and 2 . In particular, our trial wave functions for 1 and 2 states satisfy naturally the orthogonality condition. One of the key findings in this work is the criterion of the mixing angle between 1 and 2 for reproducing the correct order of the mass gap (i.e., Δ > Δ ) and decay constants (i.e., 1 > 2 ) between 1 and 2 states. Various properties of heavy (1 , 2 ) state mesons such as mass spectra, decay constants, distribution amplitudes (DAs), and electromagnetic form factors are also scrutinized. Moreover, we obtain the mass of the radially excited (2 ) state as ≈ 2600 MeV, which leaves the possibility that the 0 (2590) + observed recently by the LHCb Collaboration [16] can be interpreted as the standard quark-antiquark radial excitation. This paper is organized as follows. In Sec. II, we briefly introduce the effective Hamiltonian and the trial wave functions adopted in the present approach. We also describe how to determine the model parameters via the variational analysis. Subsequently, we describe the mass spectra of the 1 and 2 state heavy mesons and the role of the mixing is discussed as well. In Sec. III, we summarize various properties of heavy mesons including decay constants, DAs, and electromagnetic form factors obtained in our LFQM formalism. Section IV presents our numerical results for those quantities of (1 , 2 ) state heavy pseudoscalar and vector mesons. Finally, we summarize and conclude in Sec. V.

II. MODEL DESCRIPTION
The key idea in our LFQM [10][11][12][13][14] for the 1 state mesons is to treat the radial wave function as a trial function for the variational principle to the QCD-motivated effective Hamiltonian saturating the Fock state expansion by the constituent quark and antiquark. In this section, we briefly summarize our LFQM and discuss some distinguished features for the trial wave functions and their mixing angles that emerge from the inclusion of the radially excited 2 state in addition to the 1 state.

A. Effective Hamiltonian
The meson system at rest is described as an interacting bound system of effectively dressed valence quark and antiquark satisfying the eigenvalue equation of the QCD-motivated effective Hamiltonian, where¯and Ψ¯are the mass eigenvalue and eigenfunction of the¯meson state, respectively. We take the Hamiltonian in the quark-antiquark center of mass frame as where 0 is the kinetic energy part of the quark and antiquark with three-momentum k = (k ⊥ , ). The effective potential is given by [10][11][12][13] where Conf is the linear confining potential, with and being parameters to be determined later. The Coulomb potential and hyperfine interaction potential stemming from the effective one-gluon exchanges for the -wave mesons are written as We take the strong coupling constant as a parameter and S · S¯ is 1/4 and −3/4 for vector and pseudoscalar mesons, respectively. As we consider the heavy meson sector in this work, we handle Hyp perturbatively employing the contact hyperfine interaction, i.e., ∇ 2 Coul = (16 /3) 3 (r), which is a fairly good approximation for the analysis of heavy meson mass spectroscopy.
The LF wave function is represented by the Lorentz invariant internal variables = + / + , k ⊥ = p ⊥ − P ⊥ , and helicity , where = ( + , − , P ⊥ ) is the four-momentum of the meson and is the four-momentum of the th ( = 1, 2) constituent quark, which leads to the constraints 2 =1 = 1 and 2 =1 k ⊥ = 0. We assign = 1 to the quark and = 2 to the antiquark, and define ≡ 1 with k ⊥ ≡ k ⊥1 . Then the three-momentum k = ( , k ⊥ ) can be written as k = ( , k ⊥ ) via the relation, is the boost-invariant meson mass squared. Therefore, the variable transformation { , k ⊥ } → { , k ⊥ } accompanies the Jacobian factor, which we take into account for the normalization of the radial part of the wave function. The LF wave function, Ψ¯= Ψ of the state pseudoscalar and vector mesons, in momentum space is then given by where Φ ( , k ⊥ ) is the radial wave function and R¯is the spin-orbit wave function that is obtained by the interactionindependent Melosh transformation from the ordinary spinorbit wave function assigned by the quantum number . The covariant forms of R¯for pseudoscalar and vector mesons are given by [3] where˜0 ≡ √︃ 2 0 − ( −¯) 2 . The polarization vectors ( ) = ( + , − , ⊥ ) of the vector meson are given by [3] (±1) = 0, 2 + ⊥ (±) · P ⊥ , ⊥ (±) , where so that the spin-orbit wave functions R¯satisfy the unitary condition automatically, i.e., R¯ R¯ = 1.
For the 1 and 2 state radial wave functions Φ of Eq. (9), we allow the mixing between the two lowest order HO wave functions ( 1 , 2 ) by writing where 1 ( , k ⊥ ) = and is the parameter which is inversely proportional to the range of the wave function and can be used as the variational parameter in our mass spectroscopic analysis. It should be noted that the wave functions include the Jacobian factor / so that the HO bases satisfy the following normalization: From the orthonormality of Φ ( = 1, 2) defined in Eq. (13) and the unitarity of R¯, one can easily see that Φ and Ψ of Eq. (9) satisfy the same normalization as . We denote (Φ 1 , Φ 2 ) for ≠ 0 and (Φ 1 , Φ 2 ) = ( 1 , 2 ) for = 0 as "mixed" and "pure" (1 , 2 ) states, respectively. As we shall discuss below, the mixing scheme turns out to be crucial to reproduce the experimental data for both masses and decay constants of heavy mesons.

B. Variational method to effective Hamiltonian
The present LFQM for the combined analysis of the 1 and 2 state heavy mesons has several parameters, namely, the constituent quark masses ( , , , ) with being the light and quark mass, the potential parameters ( , , ), the HO parameter for each (¯) content, and the mixing angle . We first determine the values of these parameters by reproducing the mass spectra based on the variational principle. Then we compute other observables of heavy mesons such as decay constants, distribution amplitudes (DAs), and electromagnetic form factors.
Here we follow the procedure adopted in Refs. [10][11][12][13], namely, we consider the central potential 0 = Conf + Coul as well as the kinetic energy 0 in variational calculation via Then the remaining Ψ¯ hyp Ψ¯ is treated as a perturbation so that we have values common for both pseudoscalar and vector mesons of the same (¯) content. This constrains the model parameters. Since the spin-orbit wave function satisfies the exact unitarity, we have the mass eigenvalue of the meson as¯= Ψ¯ ¯ Ψ¯ = Φ ¯ Φ . The analytic forms of the mass eigenvalues ( 1¯, 2¯) for the mixed (1 , 2 ) state mesons are then obtained as where ( 1 , 2 ) = (cos , sin ), = 2 / 2 , ( ) is the modified Bessel function of the second kind of order , and ( , , ) is the Tricomi's (confluent hypergeometric) function. The mass eigenvalues for the pure (1 , 2 ) states can be read by setting = 0, i.e., ( 1 = 1, 2 = 0) in Eq. (17).
In order to explore the mixing effects and to determine the optimal value of the mixing angle , we utilize the empirical constraint on the mass gap Δ ( ) = 2 ( ) − 1 ( ) between the 1 and 2 state heavy pseudoscalar and vector mesons. The mass gap Δ ( ) for pseudoscalar (vector) mesons in our LFQM is decomposed as where we separate the four different contributions, i.e., 0 , Conf , Coul , and Hyp , to the total mass gap for the taxonomical analysis in our numerical calculations. From the available experimental data for the 1 and 2 state heavy meson pairs, ( , * ), ( , /Ψ), and ( , Υ) [22], we observe that the mass gaps between pseudoscalar mesons (Δ ) are greater than the corresponding mass gaps between vector mesons (Δ ), i.e., Δ > Δ . In our LFQM calculation, Δ Kin+Conf+Coul = Δ Kin+Conf+Coul due to the usage of common parameters for both pseudoscalar and vector mesons of the same quark flavor contents as shown in Eq. (17), and thus the mass gap is exclusively governed by the hyperfine interaction Hyp and can be readily obtained as Equation (19) combined with the relation Δ > Δ provides a very important constraint on the mixing angle . It is evident that the pure ( 1 , 2 ) states with = 0 • always leads to Δ < Δ , which shows that the introduction of the mixing is inevitable. Furthermore, one can find that the condition of Δ > Δ gives the constraint, This concludes that the lower bound of the physical mixing angle, , is determined as > = cot −1 (2 √ 6)/2 6 • .

C. Model parameters
As we have discussed in the previous subsection, the parameters of heavy mesons in the present model for (1 , 2 ) state mesons include four quark masses ( , , , ) with ( = , ), seven variational HO parameters ( ,  ,  , , , , ), three potential parameters ( , , ), and the mixing angle . The variational principle in Eq. (16) leads to a constraint in the parameter space, which relates the strong coupling constant and the other parameters i.e., = ( , , , ,¯,¯). This indicates that the variational parameters¯are automatically determined once other model parameters such as the quark masses, the strong coupling constant, the string tension, and the mixing angle are fixed.
In this study of heavy mesons, we take = 0.22 GeV, = 0.45 GeV, and the widely-used string tension = 0.18 GeV 2 [17,33,34] as inputs, which were adopted in our  previous LFQM analysis [10][11][12][13] for 1 state mesons. This leaves five parameters, i.e., ( , , , , ), to be determined. In order to determine those five unknowns, we use two masses of the 1 state heavy mesons as inputs. Among many possible choices of two input masses, we find that the use of the ( , * ) pair masses as inputs produces other meson masses well enough compared to the data. Since we have only two equations ( , * ) with five unknowns to be determined, we first try to find the best fit parameters for the pure (1 , 2 ) state case without mixing ( = 0 • ). In this case, we need to choose two input parameters from ( , , , ). Through our analyses with various combinations, we found that = 1.68 GeV and = 5.10 GeV give satisfactory results. We then obtain the remaining potential parameters, = −0.538 GeV and = 0.425, by solving Eq. (17) for ( 1 , 1 * ) using their measured values. We also note that Conf and Coul are flavor-and scale-independent so that the confining potential constant and the strong coupling are the same for all heavy mesons considered in this work. Therefore, once and are determined, the values of seven parameters are automatically computed and all the other meson masses are our predictions.
Using the same quark masses ( , , , ) and the string tension as in the = 0 • case but taking into account of the two experimental constraints, Δ Hyp > Δ Hyp and 1 > 2 , we obtain the optimum value = 12 • of the mixing angle as well as other model parameters to cover both charm and bottom flavors of the heavy quark. This would be enough for verifying the mixing effects on the physical quantities of heavy mesons.
We summarize our best fits for the model parameters obtained for the mixed state ( = 12 • ) case in Table I. For the comparison purpose of mixing effects, we also include the best fits for the model parameters obtained for the pure state ( = 0 • ) case. This shows that the values of and are not significantly different in both cases, but the values of the parameters become smaller with mixing, which results in different meson properties.
With the model parameters determined as above, we can The mixing angle in general depends on the quark flavor contents of mesons. We found = 9.8 • , 17.6 • , and 13.9 • for ( , * ), ( , /Ψ), and ( , Υ), using the measured masses in PDG [22], i.e., Δ − Δ = 62 MeV, 64 MeV, and 37 MeV, respectively. However, the paucity of data does not allow us to estimate the mixing angles for the other mesons and, therefore, we analyze the quantities assuming = 12 • , which was found to be fair in explaining the available data. compare the central potential 0 with other model calculations. In Fig. 1, we present the central potentials 0 ( ) up to 2 fm for the pure and mixed configurations ( = 0 • and 12 • , respectively) and compare them with the well-known GI model [17] and ISGW2 model [34]. We also plot the potential of the previous works of two of us in Refs. [10][11][12][13] as the CJ model. As one may expect from the similarities of the model parameters ( , , ), the central potentials obtained from the two different mixing scenarios are almost the same and they are also quite comparable with the results from the GI and ISGW2 models as well as the CJ model.
The mixing effects on the 1 and 2 state radial wave functions are shown in Fig. 2 for the bottomonium (¯) case. In Figs. 2(a) and 2(b), we compare the two radial wave functions of the pure states (dashed lines) and the mixed Φ states (solid lines) for = 1, 2. Shown in Fig. 2(c) are the ratios . This shows that the small mixing ( = 12 • ) significantly modifies the wave function of the 2 state, while the 1 radial wave function is barely modified. The mixed state radial wave functions Φ and their ratios for various heavy-heavy (¯, ) and heavy-light (¯,¯) with ( = , ) quark states are also given in Figs. 2(d)-2(f). Since the range of the radial wave function is inversely proportional to the value of the parameter, it is quite natural to observe that the wave function of bottomonium state is narrower than the other states. We also note that the value of the radial wave function at the origin ( = 0) is proportional to the parameter. As shown in Figs. 2(d) and 2(e), the bottomonium wave function at the origin is the highest among those of meson wave functions for both the 1 or 2 state. However, the ratio Φ at the origin takes the same value independent of the quark flavor contents once the mixing angle is fixed as can be seen in Fig. 2(f).

D. Mass spectra
In this subsection, we compute the mass spectra of the 1 and 2 state heavy mesons. With the parameters given in Table I, the mass formulas in Eq. (17) are used to obtain the mass spectra. The meson masses obtained for = 0 • and = 12 • cases are summarized in Table II. As discussed before, all the meson masses apart from the two inputs ( 1 , 1 * ) are our predictions. For comparison, we also list the experimental data of Ref. [22] and the predictions of the GI model [17] and those of the relativistic quark model (RQM) in Refs. [18,19]. Our predictions obtained from both mixing scenarios are found to be overall in a good agreement with the experimental data [22]. However, there are some delicate but important mixing effects in mass spectra. While there are no big differences in the predicted masses for the 1 state mesons between the pure and mixed scenarios, the predictions for the 2 state mesons obtained from the mixed case agree better with the experimental data. This also can be seen by performing 2 analysis, which gives 2 = 0.009 for the mixed scenario and 2 = 0.024 for the pure one.
In particular, our result, (2 ) = 2600 MeV for the 2 The 2 is computed as 2 = ( − ) 2 / 2 , where and are the experimental data and theoretical prediction, respectively. state of the meson with = 12 • is very close to the observed mass of the recently discovered 0 (2590) + with = 0 − [16]. This supports the interpretation of the observed 0 (2590) + state as a standard quark-antiquark radial excitation of the + meson as claimed by the LHCb Collaboration [16]. For making a definite conclusion on the structure of the 0 (2590) + , however, we still need more detailed and precise experimental studies on various properties of the 0 (2590) + .
The LHCb Collaboration [35] also reported traces of (5840) and (5960), and confirmed the observation of the (5960) by the CDF Collaboration [36]. Although their existence as resonances awaits confirmation and their quantum numbers are yet to be identified, these states are suggested as the 2 states of and * mesons in Ref. [35]. However, if these are the (2 ) and * (2 ) states, then it violates the mass hierarchy by giving Δ < Δ as Δ ≈ 561 MeV and Δ ≈ 635 MeV. This is in contradiction with our predictions on the mass gaps, Δ ≈ 612 MeV and Δ ≈ 561 MeV, which observe the relation Δ > Δ . Therefore, verifying the (2 ) and * (2 ) states is crucial to understand the structure of the radially excited heavy meson states. Experimental searches for these states in , * , and * are thus highly anticipated.
The top panel of Fig. 3 shows the mass spectra of the 1 and 2 state heavy mesons. The middle panel of Fig. 3 represents the mass gaps between the 1 and 2 state mesons and the four different contributions to this mass gap (Δ Kin , Δ Conf , Δ Coul , Δ Hyp ) are depicted in the bottom panel of Fig. 3. The dashed and solid lines in the upper and middle panels represent our results obtained with the pure ( = 0 • ) and mixed ( = 12 • ) cases, respectively. The decomposition of Δ shown in the bottom panel of Fig. 3 is for the mixed case. We also note in this taxonomical analysis that the contributions of the heavier and lighter quarks in the kinetic energy part are further separated and denoted as 'Kin 1' and 'Kin 2,' respectively, when the quark contents are different. As one can see from the mass gap Δ , the observed mass gap relation, Δ > Δ , cannot be realized without introducing the mixing angle. It is also interesting to see from the available data that the mass gaps between the 1 and 2 states are around 600 MeV, and the values are almost flavor-independent. Similar mass gap is also observed for the radially excited states of baryons with various flavors [37].
While the mass gap Δ seems almost flavor-independent as shown in the middle panel of Fig. 3, the four different contributions, (Δ Kin , Δ Conf , Δ Coul , Δ Hyp ), which make up Δ , are flavor-dependent as one can see from the bottom panel of Fig. 3. For instance, comparing the Coulomb and confinement interactions, one can easily find from Eq. (17) that Δ Conf ∝ −1 while Δ Coul ∝ . This relation provides an intuitive explanation for the observation that the confinement is dominant at large distances, while the Coulomb interaction arising from the one-gluon exchange dominates at short distances. This tendency can be clearly seen in the bottom panel of Fig. 3. By comparing the mass gap components of ( , ) or ( / , Υ), one can see that the green area for Δ Coul becomes larger for bottom quark systems. We also find that the kinetic energy is one of the important components in the mass gap. In particular, for mesons with different quark and antiquark masses such as and mesons, we find that the light quark (Kin2) gains more kinetic energy than the heavy quark (Kin1) as one may expect intuitively.

III. APPLICATIONS
With the model parameters fixed by mass spectrum, we predict various properties of the 1 and 2 state heavy mesons in this section. The standard LFQM adopting the spin-orbit wave functions and the polarization vectors of a vector meson defined in Eqs. (10) and (11) is based on the requirement of all constituents being on their respective mass shell (i.e., → 0 ). This on-mass-shell condition of quark and antiquark is completely different from the manifestly covariant models, which allow the quark and antiquark to be off-massshell allowing ≠ 0 . For instance, the invariant mass 0 included in the polarization vector (0) of Eq. (11) in the standard LFQM needs to be replaced by the physical mass in the manifestly covariant model.
The complications coming from the binding energy issue do not appear in the analysis of meson mass spectra since only the radial wave functions are needed. However, it does matter for the calculations of other physical observables such as the decay constants and form factors, which we will discuss below. In the previous works of Refs. [38][39][40][41][42] for the decay constants, DAs for pseudoscalar and vector mesons, and weak transition form factors between two pseudoscalar mesons, it was shown that the self-consistent LFQM description of those physical observables can be achieved if and only if every physical mass appeared in the matrix elements is replaced by the invariant mass 0 . In other words, the replacement of the physical mass in the integrand of the amplitude by the invariant mass 0 (denoted by CJ-scheme for convenience) results in the physical observables that are independent of the current components and polarizations used in the analysis. This → 0 mapping is indeed proven to be an effective way of including the treacherous points such as the light-front zero modes and the instantaneous contributions. As the comprehensive and rigorous analysis of decay constants and DAs for pseudoscalar and vector mesons can be found in Refs. [38][39][40][41], here we just summarize the final theoretical results for those physical quantities for completeness.

A. Decay constants
The decay constants of a pseudoscalar meson and a vector meson with a four-momentum and a mass are defined by as and , respectively, where ( ) is the polarization vector of a vector meson given by Eq. (11). For the case of pseudoscalar mesons, it has been explicitly shown that the decay constants of the pure 1 state mesons obtained from the plus ( = +) and minus ( = −) components of the currents are exactly the same [41]. In the present work, we extend it to the cases of mixed 1 and 2 states. Denoting (±) obtained from the plus and minus components of the currents, the explicit forms of (±) are given by [41] (±) = √ 6 where A = (1 − ) +¯and with A = A ( ↔¯). Here, Φ( , k ⊥ ) denotes the wave functions (Φ 1 , Φ 2 ) defined in Eq. (13) for (1 , 2 ) state decay constants.
For the vector meson case, it was also explicitly shown that the decay constants for the pure 1 state mesons obtained from the plus ( = +) component of the currents with the longitudinal polarization (0) and the perpendicular ( =⊥) components of the currents with the transverse polarizations (±) are exactly the same to each other [38]. Denoting obtained from the plus and perpendicular components of the currents as (+) and (⊥) , respectively, their explicit forms read [38] (+,⊥) = √ 6 where and LF = 0 + +¯. In Ref. [43], both (+) and (⊥) were computed in the BLFQ approach, but the two results are found to depend on the adopted component of the current, which was ascribed to the measure of rotational symmetry violation of the model. In our LFQM calculation, by using the CJ-scheme, however, we could confirm numerically that ≡ (+) = (−) and ≡ (+) = (⊥) for both (1 , 2 ) state mesons.

B. Twist-2 distribution amplitudes
The twist-2 quark DAs, tw−2 ( ) ( ), for pseudoscalar and vector mesons are related with the decay constants obtained from the plus component of the currents through [44] where tw−2 ( ) ( , ) is obtained by the k ⊥ integration of the LF wave function up to the transverse momentum scale . Here, (≥ |k ⊥ |) can be regarded as the energy scale that separates the perturbative and non-perturbative regimes. The twist-2 DA then describes the probability amplitudes to find the hadron in a state with a minimum number of Fock constituents and small transverse momentum separation. While the typical transverse momentum cutoff for the light meson sectors [38][39][40] was estimated as ≈ 1 GeV, the values of the scale for the heavy meson sectors appear shifted to the larger values as we discuss in the numerical results of Sec.IV B.
The normalized quark DA is defined as˜t w−2 The quark DAs can be usually expanded in Gegenbaur poly- which are closely related with the Gegenbaur moments ( ). The explicit relations between and ( ) can be found, for example, in Ref. [44].

C. Electromagnetic form factors and charge radii
We also compute the electromagnetic form factors of heavy pseudoscalar mesons as well as their charge radii. Our calculation is carried out by using the Drell-Yan-West frame ( + = 0) with q 2 ⊥ = 2 = − 2 . The electromagnetic form factor of the pseudoscalar meson can be expressed for the "+" component of the current as [10] ( where (¯) is the electric charge of quark (antiquark), and where k ⊥ = k ⊥ + (1 − )q ⊥ . The electromagnetic form factor is normalized as (0) = +¯, and the charge radius of the meson is calculated by  [22] and lattice simulations [45][46][47][48]. For the 2 state, we include the QCD sum rule result of Ref. [49]. (Lower panel) The ratio = 2 / 1 compared to the experimental data and lattice simulations.

A. Decay constants
In Fig. 4, we present our numerical results for the decay constants of the 1 (middle panel) and 2 (upper panel) state heavy mesons. The results obtained with the pure ( = 0 • ) and mixed ( = 12 • ) wave functions are represented by triangles and squares, respectively. For comparison, we show the available experimental data [22] and other theoretical predictions from lattice simulations [45][46][47][48]50] and the QCD sum rules [49].
While the decay constants of the 1 state mesons are rather robust against the change of the mixing angle, those of the 2 state mesons are shown to be quite sensitive to the mixing angle, and clearly a better description is achieved thanks to the mixing effects. For instance, our predictions for the heavy quarkonia obtained with the mixed wave functions are ( /Ψ(1 ) , /Ψ(2 ) ) = (390, 274) MeV and ( Υ(1 ) , Υ(2 ) ) = (666, 498) MeV. These values not only satisfy the hierarchy of 1 > 2 but also consistent with the experimental data [22]: ( Υ(2 ) ) = (689(5), 497(5)) MeV. The full numerical results of our LFQM compared with the available experimental data [22] as well as other theoretical model predictions [12,14,15,32,[51][52][53][54][55][56][57] are summarized in Tables III   and IV. Displayed in the third panel of Fig. 4 are our results on the ratio = 2 / 1 obtained with = 0 • , 6 • , and (12 ± 1) • cases, which are compared with the available experimental data [22] and the lattice simulations [48]. In particular, we present the results with = (12 ± 1) • as a band to check the sensitivity of the ratio on the variation of the mixing angle around = 12 • . As one can see from the experimental data for ( / , Υ), should be less than unity for heavy meson systems. The same constraint, < 1, for light mesons was also discussed in Ref. [9]. In our case with the pure state ( = 0 • ), most heavy mesons except ( , , ) mesons violate the constraint < 1. For the critical mixing ( = = 6 • ) case, one can see that the constraint < 1 is satisfied for all heavy mesons. However, the predictions with this critical mixing angle are still not comparable with the available experimental data. For the mixing angles = (12 ± 1) • , our results for ( / , Υ) are now quite close to the data. Further experimental measurements on the 2 decay constants are, therefore, highly desirable for testing our mixing angle effects. Fig. 5 are the normalized twist-2 DAs, tw−2 ( ) ( , ), for the 1 and 2 state heavy pseudoscalar and  --410 ± 20 ---Sum rules [56] 387 ± 7 418 ± 9 ----BS [52] 292 ± 25 459 ± 28 ---496 ± 20 BS [53] 385 -519(1) -709 -LFQM (CJ) [12] 326 360 349 369 507 529 LFQM (CJ2) [14] 353 361 389 391 605 611 vector mesons with = 12 • . Since the qualitative behaviors obtained with the pure ( = 0 • ) states are not much different from the mixed case, we do not give the results for the pure case in Fig. 5. In this figure, the convention is chosen so that the heavier quark in (¯) configuration carries the longitudinal momentum fraction and the lighter quark carries the fraction of 1 − . The DAs for the 1 state heavy pseudoscalar and vector mesons are given by the black solid and red dashed lines, respectively. We also compare our results for 1 pseudoscalar heavy mesons with the results of Ref. [59] (brown dot-dashed lines) obtained by employing a continuum approach to the hadron bound-state problem. The DAs for 1 state heavy pseudoscalar mesons are not much different from those for the corresponding vector mesons within our LFQM mainly because the parameters are same for both pseudoscalar and vector mesons. Although some quantitative differences can be found, in particular, for the and heavy-light mesons, the qualitative behaviors of our results for pseudoscalar mesons are similar to those of Ref. [59]. For the 1 state heavy quarkonia (¯,¯), although both DAs are symmetric under → 1 − , the shape of the DA is narrower for the bottomonium than the charmonium. However, the DAs for heavy-light systems become more asymmetric and more sharply peaked as the mass difference between the two constituents grows. In particular, one can see that the peak of the DA for heavy-heavy system such as the is more attracted to the center compared to the heavy-light system.

Shown in
The DAs of the 2 state pseudoscalar and vector mesons are shown by the blue solid and green dashed lines, respectively, in Fig. 5. This shows that the differences between pseudoscalar and vector mesons are more pronounced for the 2 states than for the 1 states. This tendency is opposite to that of BLFQ results [60] where the differences are more pronounced for the 1 states. We also found that the locations of the two extrema for the quarkonia systems move towards the end points as the quark mass decreases. In addition, the valley of the DAs is found to be much lower than those in Ref. [60]. Finally, we mention that the DAs for 2 state light meson sector such as ( , ) reported in Ref. [61] show similar qualitative behaviors found in the present work for the 2 state heavy mesons.
The normalized twist-2 pseudoscalar and vector meson DAs are rewritten as where the LF wave function corresponding to˜t w−2 ( ) ( , ) is denoted as Ψ tw−2 ( ) ( , k ⊥ ). Shown in Fig. 6 are the 3- dimensional (3D) plots of Ψ tw−2 ( , k ⊥ ) and Ψ tw−2 ( , k ⊥ ) for the 2 states of and mesons, respectively. Equation (32) implies that the normalized twist-2 DAs˜t w−2 ( ) ( , ) shown in Fig. 5 are obtained by the k ⊥ -integration of Ψ tw−2 ( ) ( , k ⊥ ). In our LFQM calculation with the Gaussian wave functions, we observe that |k ⊥ | → ∞ corresponds to the ultra violet (UV) cutoffs ( max ⊥ ) or energy scale around 2 GeV for ( ( ) , ( ) , ), 3 GeV for , and 4 GeV for , respectively. In other words, the wave functions for heavy-heavy systems have longer transverse momentum tails than for the heavy-light systems. On the other hand, the wave functions for heavy-light systems have deeper negative valleys than for heavy-heavy systems. This property in the twist-2 light-front wave function explains why only the twist-2 DAs for the 2 state heavy-light system have negative regions. Of course, the small UV cutoff such as | max ⊥ | < 2 GeV for heavy-heavy system may cause˜t w−2 ( ) to have negative regions as well. Similar observations are made for the 2 state vector meson DAs.
In Table V, we provide the -moments defined in Eq. (28) up to = 6 for both 1 and 2 state heavy mesons. For heavy quarkonia, the odd-moments vanish because of the symmetric shape in their DAs. For heavy-light mesons, the odd-moments reflect the asymmetry of the DAs coming from the mass difference between the quark and antiquark. As one can see, the first moment decreases as ( −¯) gets smaller. For example, we have 1 (1 ) = 0.644, 1 (1 ) = 0.614, and 1 (1 ) = 0.390. The DAs for 1 state mesons can be well reproduced with the first few moments up to = 4. However, those for 2 state mesons require much higher moments beyond = 6 to reproduce the full results. For the 2 state mesons, the first and third moments (2 ) ( = 1, 3) have negative values. This originates from the fact that the DA of (2 ) occupies more in the 0 < < 0.5 domain compared to the other DAs. In particular, our predictions of 2 , 4 , 6 = (0.179, 0.048, 0.019) for the 2 state of are comparable to (0.16, 0.046, 0.016) of Ref. [62], which are obtained from the leading twist light-front wave functions with the Cornell potential.

C. Electromagnetic form factors and radii
In Fig. 7, we present the electromagnetic form factors of 1 (solid lines) and 2 (dashed lines) state heavy pseudoscalar mesons obtained with = 12 • . For comparison, we show the available lattice simulation data of Refs. [63][64][65]. Since the form factors of heavy quarkonia ( , ) obtained from both quark and antiquark contributions vanish, we show only the contribution from the quark part for the comparison with the available lattice simulation results. This shows that our results for the 1 state ( + , + , ) mesons match well with the lattice simulation results. We find that the form factor of 2 state mesons are in general steeper than those for the corresponding 1 state mesons. In a heavy-light system such as ( + ( ) , + ), the main contribution to the form factor for the region of 2 > 6 GeV 2 comes from the heavy quark and the light quark contribution is negligible at high 2 regions. On the other hand, for the form factor of the meson, both and quark contributions are almost equally important for the intermediate 2 region. Figure 8 presents the calculated charge radii 2 of the (1 , 2 ) state heavy mesons obtained with = 0 • and = 12 • cases by triangles and boxes, respectively. Since there is no experimental data, our results are compared with the lattice simulation results of Refs. [63][64][65][66]. The full results are summarized in Table VI with other model predictions from Refs. [6,67,68]. We find that our predictions are overall in a good agreement with the lattice results. Since the heavy quark is sitting near the center of a meson while the light quark is moving actively, the radius of a heavy-light meson would be mostly governed by the motion of the light quark. This can be noticed by comparison with the results of the LFQM of Ref. [6].
We find that the charge radii for 2 state mesons are more sensitive to the mixing angle than for the 1 state cases. As shown in Fig. 8, the charge radii for 2 state mesons are larger in the mixed case than those in the pure one. The larger radii can be understood from the 2 wave functions given in Fig. 2. Namely, the reduction of the wave function at the origin results in a larger radius since the wave function is more spreading to a larger distance. This is the opposite behavior of the decay  [63], and green circles and magenta circles for B1 and C1 ensembles of Refs. [64,65], respectively. Only one quark contribution is considered for quarkonia.   Fig. 4, which is proportional to the wave function at the origin. Therefore, these observables reflect the structure of the wave functions from different points of view.

V. SUMMARY
In the present work, we have investigated 1 and 2 state heavy mesons employing the pure and mixed harmonic oscillator wave functions. We invoked the variational principle adopting the linear plus Coulomb potential, and treated the hyperfine interaction perturbatively as a contact term to distinguish vector and pseudoscalar mesons. The variational principle allowed us to obtain a constraint for model parameters.
With the fixed quark masses, all model parameters are determined by two meson masses and this led us to predict and test other physical quantities.
We have analyzed the mass spectra, decay constants, distribution amplitudes, electromagnetic form factors, and charge radii of the 1 and 2 state heavy mesons. As for the mass spectra, our predictions are in a good agreement with the available experimental data. Although no apparent difference is found for the masses of the 1 heavy mesons in the pure and mixed cases, the predicted masses of the 2 heavy mesons are appreciably modified and have a better agreement with the available data when the mixing is introduced. Our results support the speculation that the newly observed (2590) by the LHCb Collaboration [16] can be interpreted as the radially excited (2 ) state. We also observe that the mass gaps between the 1 and 2 state mesons are around 600 MeV and they are not sensitive to the flavor contents of mesons. In LFQM, we found that mass gaps of pseudoscalar mesons can be made larger than those of vector mesons regardless of the quark flavor contents only if we use the mixing angle ≥ = 6 • . Such behavior can be explained by Eq. (19), which shows that the hierarchy appears in the opposite direction without mixing.
The mixing effects are crucial to understand the properties of 2 state mesons. As for the decay constants, we could obtain a good agreement with the experimental and lattice simulation data for 1 state mesons even without the mixing effects. However, for the 2S states, the mixing effects are essential to get the correct order of decay constants. By introducing a small mixing, we noticed that the ratio Φ = Φ 2 2 /Φ 2 1 becomes smaller than unity. The optimum value of the mixing angle is obtained as 12 • to cover both the charm and bottom flavors of the heavy quark.
For the DAs of the 1 states, our prediction is found to be similar to those reported in Ref. [53]. For the 2S states, some DAs have the nodal structure arising from the structure of the wave functions. We note that the difference between DAs for vector and pseudoscalar mesons are more pronounced for the 2 states. In addition, we find that the DAs are saturated up to several GeV for the transverse momentum and it has longer tails for mesons with heavier quarks. For completeness, the corresponding -moments up to = 6 are computed in the present work.
The electromagnetic form factors and charge radii for , , and mesons are also computed and found to be comparable with the available lattice simulation data of Refs. [63][64][65]. The mixing effects lead to larger radii of the 2 states since their wave functions are more spread in space, which results in the reduction of the wave functions near the origin. This is opposite to the behavior of the decay constants that are reduced by the mixing.
In this work, we have focused on the heavy meson sector in LFQM. However, a combined analysis for both light and heavy meson sectors is also of great importance for scrutinizing the dependence of physical quantities on quark masses. In the future work, we would consider smearing the hyperfine interaction to treat it nonperturbatively as a part of the entire Hamiltonian for the application of the variational principle. A global analysis would be also required for more rigorous inves-tigations to discuss the uncertainties of the model parameters. While the investigation along this directions is under progress, more precise measurements on the physical properties of heavy mesons as well as observations of undiscovered heavy meson states are essential to test phenomenological models on the structure of heavy mesons.