Light-front wavefunctions of mesons by design

We develop a mechanism to build the light-front wavefunctions (LFWFs) of meson bound states on a small-sized basis function representation. Unlike in a standard Hamiltonian formalism, the Hamiltonian in this method is implicit, and the information of the system is carried directly by the functional form and adjustable parameters of the LFWFs. In this work, we model the LFWFs for four charmonium states, $\eta_c$, $J/\psi$, $\psi'$, and $\psi(3770)$ as superpositions of orthonormal basis functions. We choose the basis functions as eigenfunctions of an effective Hamiltonian, which has a longitudinal confining potential in addition to the transverse confining potential from light-front holographic QCD. We determine the basis function parameters and superposition coefficients by employing both guidance from the nonrelativistic description of the meson states and the experimental measurements of the meson decay widths. With the obtained wavefunctions, we study the features of those meson states, including charge radii and parton distribution functions. We use the $J/\psi$ LFWF to calculate the meson production in diffractive deep inelastic scattering and ultra-peripheral heavy-ion collisions, and the $\eta_c$ LFWF to calculate its diphoton transition form factor. Both results show good agreement with experiments. The obtained LFWFs have simple-functional forms and can be readily used to predict additional experimental observables.


I. INTRODUCTION
Understanding and describing hadrons, the bound states of quantum chromodynamics (QCD), is crucial to increase our comprehension of the strong interaction and the constitution of matter. The charmonium sector has attracted extensive experimental investigations, including the mass spectrum, transitions between excited and low-lying states, and photoproduction of the vector mesons in heavy-ion collisions [1][2][3][4][5][6]. Theoretical efforts also contribute from an array of complementary perspectives, both Euclidean and Minkowskian formalisms. Euclidean formulations of quantum field theories such as Dyson-Schwinger equations [7,8] and lattice gauge theory [9,10] offer methods of performing a first-principles computation of the charmonium spectrum and other observables. On the other hand, Hamiltonian methods formulated in Minkowski spacetime also provide a detailed description of the meson's internal structure and dynamics through wavefunctions. The wavefunctions play a central role in describing the bound states and computing the physical observables in Hamiltonian formalism.
In a standard Hamiltonian formalism, the wavefunctions are solved from the Schrödinger(-like) equations where the Hamiltonian governs the physics of the system. Since the discovery of the J/ψ resonance in 1974 [11], various potential models were developed to describe the heavy quarkonium system, including the Buchmüller-Tye potential [12], power-law potential [13], logarithmic potential [14], and the Cornell potential [15]. These phenomenological models from the early years were inspired by various aspects of QCD, and they successfully described the spectrum, especially the ψ family and the Υ family.
Later on, nonrelativistic QCD (NRQCD) was developed as an effective field theory by incorporating standard quantum field theory techniques such as dimensional regularization [16][17][18]. It captures the nonrelativistic nature of the heavy system, and the relativistic corrections can be incorporated systematically, although calculations show that the relativistic corrections may be large for selected charmonium observables [19].
One shared feature of these studies is that the meson wavefunctions are given in their rest frames. Hadrons in high-energy processes, however, are correctly described by wavefunctions on the light front. A Lorentz transformation is required to apply rest-frame wavefunctions to modern deep inelastic scattering (DIS) experiments. In principle, the light-front wavefunctions can be generated directly from the light-front Hamiltonian approach [20][21][22], which combines light-front quantization and affords boost-invariant light-front wavefunctions (LFWFs). Light-front holography (LFH) exploits the AdS/CFT correspondence between string states in anti-de Sitter (AdS) space and conformal field theories (CFT) in physical space-time to obtain a semiclassical first approximation to QCD [23,24]. It generates effective potentials for valence quarks of hadron bound states, and the resulting LFWFs are relativistic and analytically tractable. Basis lightfront quantization (BLFQ) [25] has been used to improve LFH by incorporating a longitudinal confinement and a realistic one-gluon exchange interaction [26][27][28].
A different path to obtain the meson LFWF is to model it directly or determine it from other formalisms. In such an approach, the LFWF is not solved from an eigenvalue equation and does not require one to assume a specific form for the Hamiltonian. However, the functional form can be inspired by a phenomenological Hamiltonian when modeling the LFWF directly. The vector meson wavefunction is modeled as predominantly a quark-antiquark state. The Dosch, Gousset, Kulzinger, and Pirner (DGKP) model [29], and the widely used boosted Gaussian [3,30,31] are in this category. In these parametrizations, the helicity and polarization structure is the same as in the photon perturbatively calculated in QCD. One main advantage of such modeled wavefunctions is their simplicity. They have become an important element for calculating meson production cross sections, e.g., exclusive processes at the electron-ion collider (EIC) [32]. There are also works in determining the LFWF from other formalisms. LFWFs determined from the Dyson-Schwinger and Bethe-Salpeter approach embed information from higher Fock states, which is achieved by projecting the covariant Bethe-Salpeter wavefunctions onto the light front [33][34][35][36]. Insights from NRQCD could also be carried through LFWFs by boosting the meson wavefunctions from the rest frame to the light-front frame and employing Melosh rotations on the spin structure [37,38].
As a complementary study to the existing modeled LFWFs, we propose a method of designing the LFWFs of meson bound states with a simple-functional form. A unique feature of the formalism in this work is its close relation to the light-front Hamiltonians through the choice of basis functions. More specifically, for designing the charmonium LFWFs, we choose eigenfunctions of a generalized light-front holographic confining potential as the basis functions, which were first introduced in a study on heavy quarkonia by the BLFQ approach [26,27]. Consequently, this work is closely related to the light-front Hamiltonians through the choice of basis functions, and shares several advantages of the LFH framework and the BLFQ framework. Additionally, this choice of functional forms provides us with some guidance on interpreting the meson's internal structure. While the works on and based on LFH give connections and predictions across the entire mass spectrum of hadrons [24,[26][27][28][39][40][41][42][43][44][45][46][47][48], our proposed approach has more flexibility in adjusting the wavefunctions and facilitates a better agreement with targeted experimental observables such as the dilepton and diphoton decay widths studied in this work.
The implicit Hamiltonian that describes the designed system is thereby understood as an effective Hamiltonian which is an extension of the LFH/BLFQ Hamiltonian, and its information is carried directly by the functional form and the adjustable parameters in the wavefunctions. We determine the parameters and basis coefficients by adopting guidance from the nonrelativistic description of the meson states and the experimental measurements of the meson's decay widths.
In designing the charmonia LFWFs, our major considerations and motivations are the following: • Approximation to QCD. The designed LFWFs inherit an approximation to QCD from LFH through the basis functions.
• Symmetries. The designed states are invariant under kinematical symmetries, including mirror paritym P and charge conjugation C, and approximately under rotational symmetries, including the total spin J and parity P. In addition, when calculating Lorentz invariants, rotational symmetry is respected by using different current components and different polarized states.
• Matching the NR limit. In the NR limit, the LFWFs reduce to the solutions of the Schrödinger equation with a spherically symmetric harmonic oscillator potential, and the three-dimensional rotational symmetry is restored.
• Decay widths. The LFWFs give the correct diphoton decay width for the pseudoscalar and the dilepton decay width for the vector meson. These constraints are most sensitive to the behavior of the wavefunctions at short distances, which is also the most important region for perturbative scattering processes.
We call our approach "by design", because we are choosing by hand to apply exactly the constraints that we consider to be the most important for phenomenological applications of the wavefunctions. It is an alternative phenomenological approach to obtain LFWFs that is not limited by a particular choice of Hamiltonian. This involves a certain amount of judgement as to how many basis functions to include and which constraints to impose. We make these choices in such a way that we can exactly satisfy all the constraints that we are using. Alternative schemes are also possible, for example, having the number of observables exceed the number of parameters, and one then needs to choose weights for different observables. One primary advantage of our approach is that the resulting LFWFs are analytically tractable and can be used to calculate a wide variety of physical observables.
In this paper, we calculate the charge radii and parton distribution functions of those meson states with the obtained wavefunctions. We estimate the masses of those designed states by evaluating their expectation with existing Hamiltonians, and they are in a reasonable range compared with experimental values. We use the wavefunction of J/ψ to calculate meson production in heavy-ion collisions and compare it with other model calculations. We use the wavefunction of η c to calculate the diphoton transition form factor. We found that all observables calculated are in reasonable agreement with experimental measurements.
The layout of this paper is as follows. We first introduce the formalism of the basis function representation in Sec. II. We then construct the LFWFs for J/ψ, η c , ψ , and ψ(3770) in Sec. III. With the obtained wavefunctions, we study selected key features of those meson states in Sec. IV and calculate the J/ψ production in high-energy scattering and the η c diphoton decay. We conclude the work in Sec. V.

II. BASIS FUNCTION REPRESENTATION
In this section, we introduce our formalism of designing the meson LFWFs in a basis function representation and the basis functions we use in this work.
A. Light-front wavefunctions in a basis space Consider a meson state h consisting of a quark and an antiquark, with momentum (P + , P ⊥ ), and expand its wavefunction on an orthonormal basis {β 1 , β 2 , . . . , β N β }, where C h,i are the basis coefficients for h and N β is the number of basis states. Here we are writing the wavefunction in a relative coordinate, where x = p + q /P + is the longitudinal momentum fraction of the quark and k ⊥ = k q⊥ − x P ⊥ is the relative transverse momentum.
The wavefunctions should satisfy the orthonormalization relation Physical quantities and observables (O) such as decay widths and charge radius are functions ( f O ) of the basis coefficients, The constraints Eqs. (2) and (3) form a system of equations, and the unknowns are the basis coefficients C h,i and could also include parameters in the basis functions. The procedure of designing LFWFs is, in essence, solving such a system of equations.
In the continuum limit of N β → ∞, both the number of constraints and the number of unknowns are infinite. For the purpose of designing LFWFs with a simple-functional form, it is favorable to use a small number of basis states. At the same time, the constraints are chosen from observables that are most relevant to the physics one wants the wavefunctions to describe. A solution can be determined uniquely when the number of equations is equal to the number of unknowns. It could also be true that the system is overdetermined if there are more equations than unknowns, in which case one could obtain a solution by using optimization algorithms to minimize the deviation from the constraint. If there are more unknowns than equations, such that the system is underdetermined, one might obtain multiple solutions, with additional criteria needed to choose the preferred one. The approach we choose in this paper is to work in a relatively small basis and to impose exactly enough constraints to obtain a unique solution that is a physically reasonable description of the lowest states J/ψ, ψ , ψ(3770) and η c . In the future, this work could be extended to a larger basis and to include additional constraints.

B. A generalized holographic basis representation
In this paper, we are working with a highly truncated basis space, and our goal is to describe the charmonium bound states in a simple-functional form. It is, therefore, crucial to choose a basis function that has been successfully adopted for solving light-front Hamiltonian problems. For this purpose, we take the basis introduced in BLFQ studies on heavy quarkonia [26,27]. These basis functions are the eigenfunctions of a generalized holographic confining potential, where m q (mq) is the mass of the quark (antiquark), and for charmonium we have m q = mq = m f . The first two terms are the light-front kinetic energy of the constituent quark and antiquark, the third term is a transverse holographic confining potential, and the last term is a longitudinal confining potential. The confining strength is characterized by the parameter κ. This Hamiltonian, although with independent strength parameters for the transverse and longitudinal confining terms, is referred to as BLFQ 0 in a recent work on light mesons [28].
This basis, obtained from the eigenfunctions of the Hamiltonian (4), consists of a two-dimensional (2D) harmonic oscillator (HO) function φ nm in the transverse direction and a modified Jacobi polynomial χ l in the longitudinal direction. The transverse basis function in the momentum space reads where k ⊥ = | k ⊥ |, θ k = arg(k x + ik y ), n = 0, 1, 2, . . . is the principal number, and m = 0, ±1, ±2, . . . is the orbital number. The transverse confining strength κ functions as the basis parameter here. The orthonormality relation is Note that the HO functions are a natural generalization of the Gaussian type wavefunctions widely adopted in the literature [3,29]. The longitudinal basis function is where P (α,β) l is the Jacobi polynomial with l = 0, 1, 2, . . .. The dimensionless basis parameters α and β are related to the longitudinal confining strength and the fermion mass in H 0 as α = 4m 2 q /κ and β = 4m 2 q /κ. For charmonia, α = β. The orthonormality relation is 1 4π Each basis state is characterized by five quantum numbers, {n, m, l, s,s}, where s (s) is the light-front helicity of the quark (antiquark). The basis is constructed to conserve the magnetic projection of the total angular momentum: m j = m + s +s, and m is interpreted as the orbital angular momentum projection.
The LFWF of a meson state h is written as an expansion on this basis function representation: h (n, m, l, s,s)β n,m,l,s,s ( k ⊥ , x) , (9) where ψ h (n, m, l, s,s) is the basis coefficient, and the basis function is defined as β n,m,l,s,s ( k ⊥ , x) ≡ ψ nml ( k ⊥ , x)σ ss . The spatial dependence is in ψ nml ( k ⊥ , x), which combines the transverse and the longitudinal basis functions as The parameters in the transverse basis function and those in the longitudinal basis function are connected by α = β = 4m 2 f /κ, so we have two adjustable parameters for the basis function, κ and m f . Figure 1 presents some special cases of the basis function ψ nml ( k ⊥ , x) as functions of k ⊥ and x.
The charge conjugation on the basis reads The basis states ψ nml σ ss are categorized by charge conjugation in Table II. We use those eigenstates of mirror parity and charge conjugation as building blocks to construct meson states.

Light-front spectroscopic state
We have seen that the spatial part of the basis function, ψ nml , is the eigenfunction of H 0 in Eq (4). In the nonrelativistic limit -10 -10 m f κ, the potential in H 0 reduces to the 3D harmonic oscillator potential. Consequently, the basis functions ψ nml are also related to the 3D harmonic oscillators (HOs). Such a relation is helpful in designing the LFWFs of charmonium, which has been studied extensively in the NRQCD framework. We shall first use this nonrelativistic limit to construct approximate orbital angular momentum states (denoted by ψ LF−W ) as combinations of our basis functions. We will then combine these spatial wavefunctions with the helicity structure to construct what we call light-front spectroscopic states ψ LF−n , jPC . These states are exact eigenstates of mirror parity and charge conjugation, which are good symmetries on the light front. Although they do not exactly correspond to orbital angular momentum eigenstates, they are close enough to be clearly identified with states of specific n (the principle number), j (the total angular momentum), and s (the parton spin).
We first construct the light-front spatial state ψ LF−W that has an approximate orbital angular momentum, using the basis function ψ nml . The "W" part contains the information of the orbital angular momentum and its projection m(= − , − + 1, . . . , ), and the principal number n. We take the spectroscopic notation, in which = 0, 1, 2 corresponds to the S, P, and D, respectively, and n = 1, 2, . . . labels the energy level in the ascending order. The notation LF − W distinguishes the light-front spatial state from the spectroscopic W state of the 3D HO. In the nonrelativistic limit, the former should reduce to the latter with the corresponding n, , m numbers. Building ψ LF−W , therefore, suggests a transformation between the 3D HO and the light-front functional basis ψ nml . Such a transformation is non-trivial and not exact since angular momentum is dynamical on the light front. However, considering the similarity of our basis function and the 3D HO in the cylindrical coordinates, we could use the transformation of 3D HO from the spherical coordinates to the cylindrical coordinates as a guidance. The transformation coefficients are calculated in Appendix. B.
The lowest basis state in the chosen basis function representation is the ground state ψ 0,0,0 . We expect it to provide a close approximation to the 1S state, The 1P state has three components differing by magnetic projections, and their light-front basis correspondences are In the 3D HOs, the radial excited state and the angular excited state are distinct from each other, whereas in the light-front basis, the radial and angular excitations mix automatically. The first radial excited state, the 2S state, can be approximated as a linear superposition of ψ 1,0,0 and ψ 0,0,2 . We take their relative coefficients from the cylindrical representation, which leads to The 1D state has five components differentiated by magnetic projections, and their light-front basis correspondences are These light-front spatial states ψ LF−W are orthogonal to each other. In the nonrelativistic limit, they reduce to the spherical harmonics with their corresponding angular numbers (see the derivation in Appendix C). We assume that these ψ LF−W states are the spatial components of the physical heavy quarkonia LFWFs. We now combine the spatially dependent part with the quark-antiquark helicity structure into the form of ψ LF−W σ ss . Each of these states inherits the principal quantum number n, the approximate orbital number , and its projection m from the spatial part, and the parton spin projection m s (= s+s) from the helicity part. Consequently, it acquires the m j (= m + m s ) number as well. We know from Sec. II B that we can further combine these components into eigenstates of mirror paritym P and charge conjugation C (see the eigenstates listed in Tables I and II). Then by specifying the total spin j, we can construct the light-front spectroscopic state ψ LF−n , jPC , which is identified by j PC values. The total spin is a vector summation of the orbital angular momentum and the parton spin, j = + s. Note that for the light-front basis state, the value of orbital angular momentum is approximate, and the parton spin number s is 0 (1) for the singlet σ − (triplet {σ ↑↑ , σ + , σ ↓↓ }) configuration. In the following, we find the light-front spectroscopic states for meson states with j PC = 0 −+ and 1 −− .
A pseudoscalar state with quantum number j PC = 0 −+ (recall that j is only approximated in the valence sector) could be a LF-nS state or LF-nP state. The LF-nS component is The LF-nP component is The vector state with j PC = 1 −− could be a LF-S state or LF-D state. Each LF-nS m j state has one component: Each LF-1D m j state has multiple components with different helicity structures, and the relative magnitude and phase among those components are not determined. Here we take the Clebsch-Gordan (CG) coefficients such that in the nonrelativistic limit, the light-front wavefunction components reduce to the spectroscopic state with j = 1, = 2, s = 1, We will use these light-front spectroscopic states as components in designing the charmonium LFWFs.

The Hamiltonian is implicit
Unlike the standard Hamiltonian approach, the heavy quarkonia LFWFs in this work are constructed directly on a chosen basis representation, without solving the eigenvalue equation of an explicit Hamiltonian. The advantage of doing so is that we are free from interpretating any modeled Hamiltonian as the physical Hamitonian responsible for the mass spectroscopy. Instead, we use quantities that can be directly calculated from the wavefunction. The basis space we use to design the wavefunctions is much smaller than in the computational method BLFQ. Instead, our approach here is reminiscent of a mean-field theory where a many-body problem could be reduced into an effective one-body problem, and there is a simplified reference Hamiltonian with an effective interaction correction. The basis states β n,m,l,s,s ( k ⊥ , x) are a chosen subset of the eigenfunctions of a generalized holographic Hamiltonian H 0 , as in Eq. (4). We could take H 0 as a reference Hamiltonian, and interpret the full implicit Hamiltonian as H 0 plus an effective interaction ∆H, The effective interaction ∆H should contain all the remaining interaction effects from truncated states, including interaction effects from higher Fock sectors. We could estimate the invariant mass for a designed state by calculating its expectation value of H, which gives a range instead of a definite value since ∆H is unknown. The expectation value of H 0 for a basis state β n,m,l,s,s is The expectation value of H for a state by design is therefore n,m,l δ n,n δ m,m δ l,l δ s,s δs ,s + β n ,m ,l ,s ,s | ∆H |β n,m,l,s,s .
Here we give a tentative interpretation on ∆H by approximating it as an effective one-gluon exchange interaction V OGE from Ref. [26]. The one-gluon exchange interaction includes contributions from the |qqg sector. This corresponds to approximately interpreting the implicit Hamiltonian as the BLFQ Hamiltonian, which is H BLFQ = H 0 + V OGE [26]. The effective one-gluon-exchange interaction, written explicitly, is where C F = (N 2 c − 1)/(2N c ) = 4/3 is the color factor, and the energy denominator is the average 4-momentum squared carried by the exchanged gluon q 2 = [−(k − k) 2 − (k −k) 2 ]/2 (see the details of the form and parameters in Ref. [27]).

III. WAVEFUNCTION AND THE DECAY WIDTH
In this section, we design the LFWFs for the four charmonium states, η c , J/ψ, ψ , and ψ(3770) in the basis function representation developed in Sec. II B. We determine the basis parameters and remaining coefficients in the wavefunctions by matching the calculated decay widths to the experimental values, in particular, the diphoton decay width for the pseudoscalar meson, and the dilepton decay width for the vector meson.

A. The vector meson dilepton decay and the decay constant
The vector meson can decay into a dilepton pair via a virtual photon, see Fig. 2 for the process. The transition amplitude is factorized into the decay constant, which is defined with the electromagnetic current via the local vacuum-to-hadron matrix element Here m V is the mass of the vector meson V and f V the decay constant. It is related to the experimental decay width in the particle rest frame as [49,50] Here, Q f is the dimensionless fractional charge of the constituent quarks, Q c = +2/3 for the charm quark and Q b = −1/3 for the bottom quark. The decay width sums over the contributions from all polarized states of the vector meson and the virtual photon, and the size of the individual contribution depends on the reference frame. On the contrary, the decay constant defined from Eq. (31) is Lorentz invariant and should not depend on the polarization of the vector meson or the current component in the calculation. As a result, for the LFWFs obtained from the Hamiltonian formalism or model calculations, the decay constants extracted from different polarized states could provide a measure of the violations of rotational symmetry in the system, for example, in Ref. [51]. For the formalism we are developing here, this provides a powerful constraint in designing the LFWF; that is, we require the decay constant calculated from different polarized states of the vector meson to be the same as the experimental value. In this way, we impose rotational invariance in the LFWF at the level of electromagnetic decay width.
In the LFWF representation, the decay constant extracted from Eq. (31) is an integral of the meson LFWF. The calculations using different current components in combination with different polarized states have been analyzed and discussed in Chapter. 2 of Ref. [52]. In brief, for the m j = 0 state, using the J + and the J ⊥ current components are equivalent, and the result from the J − current only agrees with the two on the condition that the meson mass equals the invariant mass of the constituent quark and antiquark. The latter is not guaranteed for the LFWF in the valence sector since the bound state mass also contains interactions besides the kinetic energy of the valence particles. We therefore adopt the J + or equivalently the J ⊥ current to calculate the decay constant from the m j = 0 state, Recall that we adopt the notations for spin configurations as For the m j = ±1 states, the J + channel is not available since e + (P, m j = ±1) = 0, and using the J ⊥ and J − currents give the same result, Here, we use the notations k R = k x + ik y = k ⊥ e iθ k and in Appendix.A 1 for the definitions of these quantities). Note that Eqs. (34a) and (34b) are equivalent up to an overall minus sign, which can be seen directly by noting that the two polarized states m j = ±1 are related by the mirror parity, as in Eq. (17). But what matters in the decay width is the absolute value of the decay constant, so the two equations give the same constraint on the transverse polarized LFWF of the vector meson.
In the basis function representation introduced in Section. II, the integrals in Eqs. (33) and (34) can be solved analytically, and those equations reduce to a summation over weighted basis coefficients. The expressions of the decay constant in the basis function representation can be found in Appendix. D.
In the nonrelativistic (NR) limit of The two choices of extracting the decay constants, Eqs. (33) and (34), reduce to the same form, In designing the LFWF of the vector mesons, we impose the constraint that the decay constant calculated from different polarized states of the vector meson to be the same as the experimental value, The significance of this condition is twofold. First, it matches the structure of the designed LFWF to the physical state, especially at short distances. Second, it enforces rotational symmetry in the vector meson state examined by an electromagnetic probe. The simplest decay channel of a neutral pseudoscalar meson for our purposes is the decay into two photons, P(P) → γ(q 1 ) + γ(q 2 ) (see Fig. 3 for the process). It involves only one transition form factor, from which the decay width Γ P→γγ can be extracted exactly at Q 2 1 = Q 2 2 = 0, i.e., taking both photons on-shell. The width can also be approximately described by a decay constant. Similar to the dilepton decay for the vector meson, the diphoton decay provides a test of the structure of the pseudoscalar meson.
In designing the pseudoscalar meson LFWF, we use the decay width calculated through the transition form factor as a constraint. However, it is also helpful to check the decay constant, especially since this approximate route has a simpler form in the basis function representation.

The diphoton decay width in terms of the pseudoscalar-photon transition form factor
The transition form factor F(Q 2 1 , Q 2 2 ) for the P(P) → γ * (q 1 ) + γ * (q 2 ) transition is defined from the Lorentz covariant decomposition of the electromagnetic transition matrix element [53,54], where P µ and q 1 µ are the four momenta of the incoming pseudoscalar meson and the outgoing photon, and q 2 µ = P µ − q 1 µ is the momentum of the other photon. Here, m P is the meson mass, and e σ,λ 1 the polarization vector of the final photon γ * (q 1 , λ 1 ), with λ 1 = 0, ±1 the magnetic projection. The transition amplitude P(P) → γ * (q 1 ) + γ * (q 2 ) is obtained by contracting I µ with the polarization vector of the other photon, where µ,λ 2 is the polarization vector of the final-state virtual photon γ * (q 2 , λ 2 ) with its spin projection λ 2 = 0, ±1. The twophoton decay width is calculated by averaging over the initial particle polarization (which is just 1 for the pseudoscalar) and summing over the final polarization (λ 1 and λ 2 ), The transition form factor defined by Eq. (37)  Here we take both choices and require the two derived diphoton decay widths to be the same as the experimental value. In this way, we impose rotational invariance in the LFWF at the level of the electromagnetic decay width, in the same spirit as what we do for the vector meson. Though the pseudoscalar meson has only one polarization state m j = 0, extracting the transition form factor using different current components and polarized states of the photon probes the pseudoscalar LFWF in different ways.
Taking the impulse approximation, in which the interaction of the external current with the meson is the summation of its coupling to the quark and the antiquark, the transition matrix elements in the LFWF representation read as an overlap of the pseudoscalar meson and the final photon γ * (q 1 , λ 1 ). Here we take the photon wavefunction calculated from light-cone perturbation theory to the lowest order (see the explicit expression in Appendix E). The transition form factor at Q 2 1 = Q 2 2 = 0 using the transverse current J + and the transverse polarized photon (λ 1 = ±1) is For simplicity, we have introduced two scalar functions that are related to the singlet and triplet helicity components as φ 0/P (k ⊥ , Using the transverse current J ⊥ and the longitudinally polarized photon (λ 1 = 0), we get According to the formalism of the covariant light-front dynamics (CLFD), the general helicity structure of a two-body pseudoscalar bound state takes the form ψ ss/P, [56,57]. It follows that the two spin components are related as With this condition, the transition form factors extracted using the two different currents, Eqs. (40) and (41), would be equivalent.
In the nonrelativistic (NR) limit of The two Equations. (40) and (41), reduce to the same form, In designing the LFWF of the pseudoscalar meson, we impose the constraint that the transition form factors calculated from different currents are the same as the experimental value converted from the diphoton decay width according to Eq. (39), In this constraint, the first equation enforces rotational symmetry in the pseudoscalar meson decay width (i.e. the integral over k ⊥ , x), but not at the level of the wavefunction [as in the CLFD condition for the integrand in Eq. (42)]. The second equation then fixes coefficients in the LFWF based on the experimental measurement.

The diphoton decay width in terms of the pseudoscalar decay constant
The decay constant of the pseudoscalar f P is defined with the axial current via the local vacuum-to-hadron matrix element, In the LFWF representation, the pseudoscalar decay constant is an integral of the meson LFWF. The calculation using the J + and the J ⊥ current components are equivalent, whereas the result extracted with the J − current has some nontrivial dependence on the meson's momentum (see Chapter. 2 of Ref. [52] for more details). Here, we take the J + (or equivalently the J ⊥ ) current to calculate the decay constant of the pseudoscalar, The expressions of the decay constant in the basis function representation can be found in Appendix. D.
The decay constant of the pseudoscalar meson is related to the diphoton decay width through a single pole fit to the transition form factor [49,50], f P ≈ F Pγ (0, 0)m 2 P /(4Q 2 f ) in the leading order approximation, thus We could also see this approximation explicitly in the LFWF representation by comparing Eq. (46) to Eq. (41). The lefthand side of Eq. (47) would reduce to the right-hand side in the nonrelativistic limit of m P = 2m f and m f >> k ⊥ .
The decay constant f η c could also be extracted from the B meson decay channels by assuming factorization, As noted above, to design the pseudoscalar LFWF that gives the correct diphoton decay width, we use the more accurate equations from the transition form factor rather than the decay constant here. Still, calculating the decay constant as a simple analytical form provides an insightful approximation and a cross-check, especially for heavy systems.
We list the experimental values of the meson masses and the dilepton (diphoton) decay width for the vector (pseudoscalar) meson, of selected charmonium states, according to the Particle Data Group (PDG) [64] in Table. III. The decay constant of the vector meson is converted from the associated dilepton decay width according to Eqs. (32) with α em = 1/134. Note that we take the effective electromagnetic coupling at the meson mass scale. The decay constant and transition form factor at zero four-momentum transfer of the pseudoscalar meson are converted from the corresponding diphoton decay width according to Eqs. (47), and (39) with α em = 1/137. Note that for this process, the coupling is taken at zero momentum transfer.
C. The LFWF of J/ψ as a 1 −− state In the charmonia system, J/ψ is the ground vector state with J PC = 1 −− . The dominant contribution to the J/ψ state should be the ground light-front state, 1S state. For instance, the BLFQ calculations predict that the probability of finding the J/ψ in a LF-1S state is larger than 99% [26,27]. In consideration of that, we set the LFWF of J/ψ to be a pure LF-1S state, ψ LF−1S ,1−− , such that all three polarization states of J/ψ  (4) 0.262 (18) 0.0977 (34) have the same spatial dependence:   Note that these are not the mass eigenvalues of the designed J/ψ state, since we are not working with some specific Hamil-tonian. These estimated masses are reasonably close to the experimental value in Table. III. Based on the nonrelativistic limit, we would expect the η c to be mostly a ψ LF−1S ,0−+ state with partons in the spin singlet state, i.e. ψ 0,0,0 σ − . We would also expect the ψ to be predominantly a ψ LF−2S ,1−− state, and thus have a large ψ 1,0,0 component [see Eq. (21)] with partons in the spin triplet state, i.e. ψ 1,0,0 {σ ↑↑ , σ + , σ ↓↓ } for m j = 1, 0, and −1 states respectively. It turns out that out that the pseudoscalar meson decay constant calculated from a ψ 0,0,0 wavefunction, and the vector meson decay constant calculated from ψ 1,0,0 have exactly the same expressions as a function of κ and m f as we have for the J/ψ [also see Eqs. (D1), (D2), and (D3) and the associated discussions]. Since these LFWF's are expected to be close to η c and ψ , it is interesting to plot also the curves corresponding to the experimental values f ψ and f η c in the same κ, m f plane. These curves are also shown in Fig. 4. One immediately observes that, having fixed κ and m f using the J/ψ decay, the decay widths of the ψ 0,0,0 σ − and ψ 1,0,0 states do not match the experimental values for η c and ψ . Therefore, we must allow both of these mesons to have some contribution from other basis states to fit their decay widths. In the previous section, we have constructed the LFWF for J/ψ as 1 −− state in Eq. (48). We find the values of κ and m f by fitting the calculated J/ψ decay constant to the experimental value, and we will use the same values in constructing other states. The pseudoscalar meson η c with quantum number 0 −+ is the ground state in the charmonium system. It is a 1S state in the nonrelativistic limit. In this work, we build the η c LFWF as a relativistic bound state, admitting a predominant LF-1S component with admixtures of LF-2S and LF-1P components, The LF-1P component has a pure relativistic origin. It is forbidden in the NR limit by the nonrelativistic parity relation P = (−1) +1 = 1 ( = 1 for the P wave), but allowed in our case as long as the charge conjugation and mirror parity symmetries are satisfied. We determine the values of the basis coefficients by matching the diphoton decay width Γ η c →γγ to the experimental value. Using the framework developed in Sec. III B, we impose the constraint that the diphoton decay width calculated through the transition form factor F η c γ (0, 0)| J + and F η c γ (0, 0)| J ⊥ match the PDG value simultaneously, as illustrated in Eq. (44). For comparison, we also use the decay constant formula, Eq. (46), to determine the basis coefficients by fitting to the diphoton decay width. In this case, only the spin singlet component is involved, and the LF-1P component is absent, i.e., C η c ,1P = 0. The results are presented in Table. IV. The obtained η c LFWF, for which the basis coefficients are listed in the third column, has a large LF-1S component and small LF-2S and LF-1P components. This set of basis coefficients are very close to the one obtained via the decay constant fit. However, there is an essential difference between the two; as we have discussed in Sec. III B, the LF-1P component is necessary for preserving the rotational symmetry detected via the diphoton decay. In comparison, the BLFQ LFWF has a smaller percentage for the dominant LF-1S component, leaving probability available for higher excited modes. We present a plot of the spatial wavefunction of η c in Fig. 6. In addition, we compare the relation of the two helicity components to the condition from CLFD in Eq. (42) in Fig. 7. We see that the actual σ ↑↑ component in our LFWF and the one obtained from the σ − component via the CLFD condition have the same overall form and roughly the same magnitude, but are not exactly equal.
E. The LFWF of ψ as a 1 −− state As the first excited vector state, ψ should receive a dominant contribution from the LF-2S state. However, mixing with other states is also required to obtain the measured decay constant. Thus, we allow admixture of the the LF-1S state. Hence, we construct ψ as a combination of LF-1S and LF-2S  (49). The values of basis coefficients are obtained by fitting the diphoton decay width Γ ηc→γγ to the PDG value (see Table. III) at κ =κ, m f =m f . Uncertainties come from those of the parameters and the PDG decay width. In the second column, the values of basis coefficients are obtained by fitting the diphoton decay width via the decay constant formula according to Eq. (46). There, the LF-1P component is taken away. In the third column, the values of basis coefficients are obtained by fitting the diphoton decay width via the transition form factor F ηcγ (0, 0) according to Eqs. (40) and (41). The corresponding LFWF is what we provide from this work. The percentage of each basis state is calculated by taking the square of the corresponding basis coefficient. The BLFQ values are extracted from the LFWF that is solved from the Hamiltonian formalism [27]. The meson mass is estimated by evaluating the expectation value of H according to Eq. (27).
value by f ηc value by F ηcγ (0, 0) percentage (%) BLFQ percentage (%) [27], normalized C ηc,1S 0.987 (7) (2) 3.02 (5) (GeV -1 ) states,  Table. V. The LFWFs of ψ are presented in Fig. 8. In addition, we compare different m j states of ψ in the longitudinal and the transverse dimension separately in Fig. 9. The ψ (m j =0) state and the ψ (m j =±1) state largely resemble each other in both the x and k ⊥ directions. This is expected by looking at their spatial decomposition, each being a predominantly (> 90%, see Table. V) LF-2S state. The ψ (m j =±1) state has a slightly higher peak at {x = 1/2, k ⊥ = 0 ⊥ } compared to the ψ (m j =0) state, due to the fact that the ψ (m j =±1) state has a somewhat larger LF-1S component (9.2% compared to 5.7%, see Table. V). The vector meson ψ(3770) is recognized as primarily a 1D wave, with admixtures such as 1S and 2S waves, in reference to potential models and BLFQ calculation. Here, we design the ψ(3770) state as a linear combination of the LF-1D, LF-1S, and LF-2S states, We impose two constraints: (1) the LF-1D component is dominant, (2) ψ(3770) is orthogonal to ψ . With these considerations, we solve the basis coefficients of ψ(3770) by fitting its decay constant to the PDG value [i.e., imposing Eq. (36)], at κ =κ, m f =m f . We also estimate the meson mass by calculating the expectation value of H defined in Eq. (27). The results are presented in Table. VI. We present the different spin components of the ψ(3770) state in Fig. 10. Those plots display the features one would expect from 1D and 2S states.

IV. MESON STATES AND OBSERVABLES
In this section, we calculate several observables for the charmonium states using the LFWFs we have constructed. For the η c observables, we adopt the amplitudes presented in the third column of Table IV. We first estimate the masses of those states, and then calculate their charge radii and parton distribution functions. We use the J/ψ wavefunction to calculate exclusive meson production in DIS and ultraperipheral heavy-ion collisions and compare with other model calculations. We use the η c wavefunction to calculate the diphoton transition form factor and compare it with the experimental measurement.
A. An estimated mass spectrum When designing the LFWFs in Sec. III, we estimate the masses for the constructed states by evaluating their result-   (51). The values of basis coefficients are obtained by fitting the decay constant to the PDG value, at κ =κ, m f =m f . Uncertainties come from those of the parameters and those of the PDG decay constants. The percentage of each basis state is calculated by taking the square of the corresponding basis coefficient. The BLFQ values are extracted from the LFWF that is solved from the Hamiltonian formalism [27]. The meson mass is estimated by evaluating the expectation value of H according to Eq. (27). x ing expectation values for the operators H 0 and H BLFQ . We now use those values to construct an estimated mass spectrum. The results are shown in Fig. 11, with a comparison to the PDG values in Ref. [64] and the BLFQ values in Ref. [27]. For each meson state, we take its values for √ H 0 and H BLFQ calculated from the m j = 0 state, and draw them in dashed and dotted lines respectively. For all the four charmonium states, √ H 0 > H BLFQ , showing that the one-gluon-exchange operator reduces the mass expectation values. We draw vertical solid lines that extend from

B. Charge radii
The charge radius of the meson bound state is defined in terms of the slope of the charge form factor at zero momentum transfer. It provides important insight on the spatial structure of the system, With the constructed LFWFs, the form factors can be obtained from the Drell-Yan-West formula within the Drell-Yan frame where q = P − P, and Q 2 = −q 2 = q 2 ⊥ . For the pseudoscalar η c , it directly produces the charge form factor G 0 (Q 2 ) = I 0,0 (Q 2 ). For vector mesons, we adopt the prescription of Grach and Kondratyuk [65], where η = Q 2 /(4M 2 h ), M h is the mass of the hadron. Table VII lists the r.m.s. radii of the states studied in this work. The values calculated from the BLFQ wavefunctions [27] are also listed for comparison. From our results, the radius of η c is about 46% larger than the radius of the J/ψ. This is because the J/ψ is a LF-1S state by design, whereas η c admits small LF-2S and LF-1P components. Both the LF-2S and LF-1P states have a larger radius than the LF-1S state.

C. Parton distributions
The quark Parton Distribution Functions (PDFs) represent the probability of finding a quark carrying momentum fraction x in the hadron state. They are essential ingredients in describing the hadron structure for hard scattering processes such as deep inelastic scattering from a hadron target and the hadron-hadron Drell-Yan process [66,67]. In the light-front formalism, the PDFs of a hadron state can be evaluated by integrating out the transverse momentum of the wavefunction overlaps. For the pseudoscalar meson η c , the quark's PDF is defined as This PDF is normalized to unity in our model in the sense that there is one valence quark, For vector mesons, the PDFs are defined in terms of the quark-quark correlation function probed in the spin-one target lepton-hadron scattering [66][67][68][69][70][71]. The PDFs appear as the , indicating an open range for the estimated masses of the constructed states. The PDG values in the red dots are experimental measurements from Ref. [64]. The BLFQ values in the yellow triangles are solved from the Hamiltonian approach in Ref. [27]. parametrization coefficients in front of the Dirac γ-matrices, and there are four time-reversal even distributions at the lead-ing twist. We write those PDFs in terms of the light-front helicity matrix elements, defined as (57) where the initial (final) state helicity of the hadron is m j (m j ) and that of the quark is s (s ). The unpolarized PDF f 1 (x) is expressed as It represents the unpolarized quark distributions in the unpolarized spin-one hadron. Similar to the pseudoscalar PDF q(x) in Eq. (55), f 1 (x) is also normalized to unity in our model, The tensor polarized PDF f 1LL (x) represents the difference of unpolarized quark distributions in the transversely polarized spin-one hadron with spin projection m j = 0 and m j = ±1, and is sensitive to the quark's orbital angular momentum [72][73][74]. It is expressed as The longitudinally polarized PDF h 1 (x) is expressed as It describes the distribution of the longitudinally polarized quark in the longitudinally polarized meson. The transversely polarized PDF g 1 (x) is expressed as It describes the distribution of the transversely polarized quark in the transversely polarized meson. Note that there are different conventions in naming the polarized PDFs h 1 (x) and g 1 (x). For example, in Ref. [75], the former is referred to as the transversely polarized PDF whereas the latter as the longitudinally polarized PDF, which is linked with calling the λ = 1 nucleon a longitudinally polarized state. Figure 12 shows the PDFs of η c , J/ψ, ψ and ψ(3770). The PDF of η c is peaked at x = 1/2, reflecting its structure as a predominantly LF-1S wave. For J/ψ, the three PDFs h 1 (x), g 1 (x), and f 1 (x) are identical, and the tensor polarized PDF f 1LL (x) is 0. This is because the designed J/ψ is a pure LF-1S wave with spin-triplet configuration, so its spatial dependence of each polarized state is the same. In the case of ψ , the three PDFs h 1 (x), g 1 (x), and f 1 (x) are only slightly different, and they reflect excitations in the longitudinal direction from the large LF-2S component. The ψ 's tensor polarized PDF f 1LL (x) deviates slightly from 0, indicating the resemblance between the transversely and longitudinally polarized states. The PDFs of ψ(3770) admit several interesting features. The unpolarized PDF f 1 (x) of ψ(3770) has an extensive flat region, indicating its large angular excitation. A recent study on ρ meson PDF also reveal a large flat region using the light-front holographic wavefunction [71]. Compared to J/ψ and ψ , the polarized PDF h 1 (x) and g 1 (x) of ψ(3770) are each different from f 1 (x). In addition, the tensor polarized PDF f 1LL (x) of ψ(3770) has a much larger amplitude, indicating the difference between the spatial dependences of different polarized states.

D. Vector meson production
In this section, we study exclusive charmonium production in diffractive deep inelastic scattering and ultra-peripheral heavy ion collisions within the dipole picture. We employ the LFWF of J/ψ obtained in Sec. III. We also make calculations of charmonium production using the BLFQ LFWF [26] and boosted-Gaussian LFWF for comparisons.
For exclusive heavy quarkonium production in DIS, the am-plitude in the dipole model can be calculated as [3] A γ * p→hp where T and L specify the transverse and longitudinal polarization of the virtual photon (with virtuality Q 2 ) and the produced quarkonium, and t = −| ∆ ⊥ | 2 is the momentum transfer squared. On the right-hand side, the transverse size of the color dipole is denoted by r ⊥ , the LF longitudinal momentum fraction of the quark is denoted by x, the impact parameter of the dipole relative to the proton is denoted by b ⊥ and x B is the Bjorken variable. Here, ψ γ and ψ h are the LFWFs of the virtual photon and the exclusively produced quarkonium, respectively (see the explicit expression of the photon LFWF in Appendix E). The cross section is related to the amplitude via Furthermore, we implement two phenomenological corrections in the calculation of the cross section: the contribution from the real part of the scattering amplitude [3], and the skewedness correction [76], which takes into account the fact that two gluons interacting with the dipole are carrying slightly different momentum fractions (consult Ref. [77] for the details of the implementations). We take the leading order perturbative calculations to obtain the photon's wavefunction. With the J/ψ LFWF from Eq. (48), the overlap functions are, written explicitly, with φ(r ⊥ , x) =κ 4(2α + 1) where e = √ 4πα em , Q f = Q c = 2/3 is the charge fraction carried by the quark, K 0 the modified Bessel function of the second kind, 2 ≡ x(1 − x)Q 2 + m 2 f and N c = 3 the number of colours. In Fig. 13, we show the overlap between the photon and the J/ψ wavefunctions integrated over x at different photon virtualities. To be precise, we plot the quantity We take three different versions of the J/ψ LFWF, the one designed in this work (denoted as "1-basis" BLFQ in the figure), the BLFQ LFWF [26], and the boosted Gaussian with x PDFs ψ' m c = 1.27 GeV [5]. The LFWF in this work is built from the same basis functions used in BLFQ. However, the BLFQ LFWF is solved by diagonalizing the light-front Hamiltonian in a much larger basis space and fitting to the meson mass spectrum, whereas our LFWF and the boosted Gaussian are obtained by fitting to the decay widths. Our LFWF exhibits similarity with the boosted Gaussian by having a simple analytical form, but the structures are different. From Fig. 13, we see that the J/ψ-photon overlaps calculated with the three J/ψ LFWFs are of the same magnitude and similar overall shape at each Q 2 . In the case of the longitudinally polarized states, as in Fig. 13(a), the overlaps from this work and the boosted Gaussian are very close and slightly higher than BLFQ. For the transversely polarized states, as in Fig. 13(b), the result from this work is roughly in between that from BLFQ and boosted Gaussian. We use the bCGC dipole model for the dipole cross section, dσ qq Here, Q s ≡ Q s (x B ) = (x B,0 /x B ) λ s /2 Q 0 and Q 0 = 1 GeV; γ s , κ s , λ s are parameters to be determined by inclusive DIS data [78]; A and B should be evaluated by continuity conditions at r ⊥ Q s = 2. We use one of the parametrizations in Ref. [79] for this investigation, which we provide in Table VIII.
We then calculate the J/ψ production in the kinematic range of the HERA experiment [80,81]. Various cross sections obtained as a function of the kinematic variables Q 2 , W, and t reasonably agree with experimental data. As an illustration, we present some representative results in Fig. 14, to-  [79]. gether with calculations using BLFQ and boosted Gaussian wavefunctions for comparison. In all three panels, the solid curves are calculated with the J/ψ's LFWF designed in this work, the dotted curves are calculated with BLFQ vector meson LFWF, and the dot-dashed curves are calculated with the boosted Gaussian LFWF of Ref. [5], respectively. The bCGC parametrization listed in Table VIII for dipole cross section was used for all wavefunctions. Figure 14(a) shows the total J/ψ cross section as function of (Q 2 + m 2 J/ψ ) for photon-proton c.m. energy W = 90 GeV. In Fig. 14(b), we show the total J/ψ cross section as function of W at various values of Q 2 . The differential cross section dσ/ dt is shown in Fig. 14(c) as function of the momentum transfer t. Qualitatively, all three wavefunctions provide reasonable descriptions of the J/ψ cross section data at HERA. The calculations using the small basis LFWF give very similar results to those using the boosted Gaussian LFWF. The BLFQ LFWF calculation generally underestimates the J/ψ production at HERA, especially in the small Q 2 regime. However, the J/ψ cross section at small Q 2 may have a stronger dependence on the dipole cross section model and the photon wavefunction [77].
We also apply the proposed J/ψ LFWF to calculate the coherent production of J/ψ at LHC at mid-rapidity in Fig. 15, using the same procedure we adopted in Refs. [77,85]. Here, the (r Solid lines are obtained using the small-basis LFWF (this work), dotted lines the BLFQ LFWF [26], and dot-dashed lines the boosted Gaussian LFWF [5].
solid curve, the dotted curve, and the dot-dashed curve show the predictions of the small-basis LFWF, the BLFQ LFWF [26], and the boosted Gaussian LFWF, respectively, for the coherent production of J/ψ in Pb-Pb ultra-peripheral collision at √ s NN = 2.76 TeV, compared to the measurements of the ALICE [82,83] and CMS collaborations [84] at the LHC. The bCGC parametrization listed in Table VIII for the dipole cross section was used for all wavefunctions. The prediction of the small-basis LFWF is within the statistical uncertainty of the experimental data. The prediction of the boosted Gaussian LFWF slightly overshoots the data, and that of the BLFQ LFWF underestimates the experimental data. Based on the above discussion, we arrive at the conclusion that the small-basis LFWF for J/ψ designed in this work can make quantitatively reasonable predictions for diffractive charmonium production in both ep collisions and ultraperipheral collisions.
The ψ cross section calculated with our designed LFWF for ψ is, on the other hand, far below the experimental data. In contrast to the ground state J/ψ, the ψ wavefunction has a node in the r ⊥ direction, so in calculating the scattering amplitude in the dipole model, there is a cancellation between the negative and the positive regions. The value of the ψ cross section is therefore very sensitive to the location of the node relative to the typical transverse separation of the virtual cc pair in the dipole model. This cancellation turns out to be very dramatic with our designed ψ LFWF, resulting in a greatly suppressed cross section. For example, due to the sensitivity to the location of this node, a 5% increase (decrease) in κ (with no other changes) results in a factor of ∼ 3 increase (a factor of ∼ 5 decrease) in the ψ production in Pb-Pb ultra-peripheral collision at √ s NN = 5.02 TeV. On the contrary, the J/Ψ production is insensitive to such a change of κ. We do not present the ψ production in this paper, and we hope to return to this aspect of the ψ LFWF in future work.
E. The γ * γ → η c transition form factor The η c meson is produced at e + e − colliders in the process e + e − → e + e − η c via the diphoton production mechanism [86]. In this section, we study the γγ * → η c transition with the LFWF of η c obtained in Sec. III. We have defined the diphoton transition form factor F Pγ (Q 2 1 , Q 2 2 ) for the process P(P) → γ * (q 1 ) + γ * (q 2 ) in Sec.III B, and used its value at Q 2 1 = Q 2 2 = 0 to determine the basis coefficients in the η c LFWF. We now use the obtained η c LFWF to calculate the transition form factor for the case of one photon on-shell and the other being spacelike.
The transition form factor can be extracted with either the J + or the J ⊥ current, and the results should be the same by Lorentz invariance. Therefore, using both currents for the calculation would help check the rotational symmetry embedded in the LFWF. We take the Drell-Yan frame, such that the vertex photon has zero longitudinal momentum, i.e., q + 2 = 0, which is the preferred frame for LFWF in the valence sector [87]. In the LFWF representation, the transition form factor extracted from the J + current reads
y FIG. 15. Predictions of the small-basis LFWF (solid curves), the BLFQ LFWF [26] (dotted curves) and the boosted Gaussian LFWF [5] (dot-dashed curves) for the coherent production of J/ψ production in Pb-Pb ultra-peripheral collision at √ s NN = 2.76 TeV, compared with the measurements by the ALICE collaboration [82,83] and CMS collaboration [84] at LHC. Error bars show statistical uncertainties only. and the expression from the J ⊥ current In both expressions, By Bose symmetry, the transition form factor should also be symmetric under the exchange of Q 2 1 and Q 2 2 . This means that in the limit of one on-shell photon, we should have F η c γ (Q 2 ) ≡ F η c γ (Q 2 , 0) = F η c γ (0, Q 2 ). However, such a symmetry is not explicit in the expressions (69) and (70). We examine this symmetry by taking both the two limits of Q 2 1 = 0 and Q 2 2 = 0. We present the results for the η c diphoton transition form factor in the format of the normalized transition form factor |F η c γ (Q 2 )/F η c γ (0)| in Fig. 16. We found that the results calculated with the two different current components, and by taking the two limits of Q 2 i = 0(i = 1, 2) agree. This indicates that the designed η c LFWF quite closely preserves both the rotational symmetry and the Bose symmetry for this observable. We also compare our results with the experimental data from BaBar, finding a reasonable agreement.  Table. IV). The solid line is obtained by taking the limit of 0)], and the dashed by taking Q 2 The BaBar experimental data is taken from Ref. [86].

V. SUMMARY
We proposed a method to build the LFWFs of meson bound states on a small-sized basis function representation. In this work, the basis functions are the eigenfunctions of an effective Hamiltonian developed from light-front holography. However, the basis coefficients and parameters of the basis functions are obtained using constraints on the wavefunction directly, not through the Hamiltonian. We use physical constraints, including the orthonormalization relation, symmetries, insights from nonrelativistic state identification, and decay widths from experimental measurements, to determine the parameters within the basis function and the basis coefficients. The resulting LFWFs inherit the physical interpretation of relativistic bound states from the phenomenological Hamiltonians of LFH and BLFQ, while admitting simple-functional forms that are feasible in calculating observables. We have adopted a "by design" approach where we choose by hand a set of sufficiently many phenomenologically most important constraints to achieve a unique determination of the parameters.
With this formalism, we designed the LFWFs for η c , J/ψ, ψ , and ψ(3770). First, we make the assumption that the J/ψ state is the ground vector state of the charmonium system and its spatial wavefunction is a LF-1S state. The two adjustable parameters in the basis function are determined by the J/ψ decay constant, which is established by its experimental dilepton decay width. We then construct the pseudoscalar state η c and the two excited vector states ψ and ψ(3770) as superpositions of the ground and excited light-front basis states with proper spin structure assignments. The basis coefficients are determined by their diphoton/dilepton decay widths, as well as other theoretical considerations, such as spatial symmetries and the orthogonality conditions. This step-by-step approach could be extended by including further constraints and could be developed towards a simultaneous global analysis in the future.
Using the charmonium LFWFs by design, we calculate several physical quantities that are accessible by experimental measurements. For instance, we calculate the charge radii and parton distribution functions; we also estimated the masses of the charmonium states we designed by evaluating the expectation value of an approximated Hamiltonian guided by the full BLFQ formalism. Our predictions for masses, charge radii, and parton distribution functions using the constructed LFWFs are in reasonable agreement with the experimental measurements and are quantitatively consistent with other established methods such as Dyson-Schwinger Equations and Lattice calculations. Furthermore, we calculate the J/ψ production in DIS at HERA and UPC at LHC using the constructed LFWF in the dipole model, the theoretical prediction is consistent with experiment data within uncertainties, and are comparable to calculations using the boosted Gaussian and BLFQ LFWFs. We calculate the η c diphoton transition form factor using the obtained LFWF and find a reasonable agreement with experimental data.
With this work, we provided light-front wavefunctions of mesons in a simple-functional form while retaining physical interpretation and matching to a variety of selected experimental observables. We anticipate these analytical LFWFs can be used for making predictions of various physical processes that involve the meson states, e.g., exclusive processes at the EIC [32].

Light-Front coordinates
The light-front coordinates are defined as x µ = (x + , x − , x 1 , x 2 ), where x + = x 0 + x 3 is the light-front time, x − = x 0 − x 3 is the longitudinal coordinate, and x ⊥ = (x 1 , x 2 ) are the transverse coordinates. We also write the transverse components with subscript x (y) in place of 1 (2), for example, r ⊥ = (r x , r y ). The covariant vectors are obtained by x µ = g µν x ν , with the metric tensors g µν and g µν . The nonzero components of the metric tensors are, The inner product of two 4-vectors is therefore For a transverse vector k ⊥ = (k x , k y ), we will write its complex forms as where k ⊥ = | k ⊥ | and θ k = arg k R .

Polarization vectors
The polarization vector for a vector boson with momentum k µ , mass m, and helicity λ is where ⊥ ± = (1, ±i)/ √ 2. The polarization vector for a photon with momentum q µ , virtuality Q 2 = −q 2 > 0, and helicity λ is µ λ=0 (q) =( In this Appendix, we calculate the transformation coefficients of the three-dimensional harmonic oscillators (3D HOs) from the spherical coordinate to the cylindrical coordinate. The meson states, especially heavy quarkonia, are sometimes identified as 3D-HO states in the spherical coordinate, e.g., 1S wave, 1P wave, etc. On the other hand, the light-front basis functions we take, as in Eq. (9), are very similar to the 3D HOs in the cylindrical coordinate. Therefore, we take the transformation between the two coordinates as a guidance to help us construct and identify 3D-HO states in the chosen light-front basis functions.
A vector in the Cartesian coordinates is given by η = (η x , η y , η z ) = η x e x + η y e y + η z e z . Then in the cylindrical coordinates, η = (η ρ , η z , θ), and in the spherical coordinates, η = (η, θ, φ), where The 3D-HO state in a spherical representation is written as |n, , m , where n is the radial quantum number, and m are the orbital angular momentum and its z component, and the total energy of the state is N = 2n + , Φ n, ,m (η) ≡ n, , m|η, θ, φ =(−1) n 4π2 (n + )! n!(2n Note that the quantum numbers (n, , m) here are different from the (n, m, l) in the light-front basis functions as in Eq. (9). Note especially that the "l" in the latter set is the light-front longitudinal quantum number. The spherical harmonics read with the phase convention where P m are the associated Legendre polynomials without the Condon-Shortley phase (−1) m (to avoid counting the phase twice). The 3D-HO state in a cylindrical representation is written as |n ρ , m, n z , where n ρ is the quantum number for the ρ coor-dinate, and the total energy of the state is N = 2n ρ + n z + |m|, Φ n ρ ,m,n z (η) ≡ n ρ , m, n z |η ρ , η z , θ = (−1) n ρ 2π 2 2n ρ +m n ρ !(n ρ + m)! 1 n z ! η n z z η 2n ρ +m ρ 1 √ 2π e imφ . (B4) The relation of the two representations is |n, , m = n ρ ,n z δ 2n+ ,2n ρ +n z +|m| n ρ , m, n z |n, , m |n ρ , m, n z , (B5) and the transformation coefficient is derived in Ref. [88]. Written explicitly, n ρ , m, n z |n, , m =δ 2n+ ,2n ρ +n z +|m| (−1) n+m+n ρ where and The transformation coefficients for selected states are calculated according to Eqs. (B5) and (B6) and listed in Tables. IX, X and XI.  In this Appendix, we write out the nonrelativistic (NR) limit of the spatial part of the light-front spectroscopic states, ψ LF−W , and compare them to the spherical harmonic oscillators. In the NR limit of | k| m f and x → 1/2 + k z /(2m f ), the light-front basis function ψ nml reduces to where C nml is a positive constant depending on n, m, and l, C nml =κ −1 4πn! (n + |m|)! 4π(2l + 2α + 1) Γ(l + 1)Γ(l + 2α + 1) Γ(l + α + 1) 2 2 −α , and R(k) is the Gaussian part, Here the spherical polar angles θ and φ represent the colatitude and azimuthal angle of k, respectively. To compare the NR limit of the basis function with the spherical harmonic oscillators, we expand them in the small momentum region. The solutions of the Schrödinger equation with a spherically symmetric harmonic oscillator potential are in the form of Ψ n m (r, θ r , φ r ) = R n (r)Y m (θ r , φ r ), where R n (r) is the radial function. In the small momentum limit of | k| → 0, the spherical harmonic oscillator in the momentum space has the following dependence on k, lim k→0Ψ n m (k, θ, φ) = f (n, )k Y m (θ, φ) , where f (n, ) is some function of the n and inherited from the radial wavefunction. In the spectroscopic notation, = 0, 1, 2 correspond to the S, P, and D wave, respectively, and n = 1, 2, . . . labels the energy level in the ascending order. For each LF spectroscopic state in the NR limit, we identify its quantum number by comparing its dependence on k to the right hand side of Eq. (C4), and n by its energy level.
The NR limit of the LF-1S state is that of the ground basis state ψ 000 ψ LF−1S ,NR = lim k→0 ψ 000,NR ( k) = C 000 = Ck 0 Y 0,0 (θ, φ) . (C5) Here and in the following, we use C (and C , C , . . .) to indicate a positive constant. Thus the NR limit of the LF-1S state is indeed a 1S state, by seeing the radial number n = 0 and Y 0,0 as the S wave.

(C15)
Note that there is a '-' sign on the NR limit of the LF-1P-1 state. It follows that the NR limit of the LF-1P1 state is the m = 1 component of the 1P state, and that of the LF-1P-1 state is negative of the m = −1 component of the 1P state. The above states are listed in Table. XII.