# Pseudoscalar meson decay constants and distribution amplitudes up to the
twist-4
in the light-front quark model

Ahmad Jafar Arifi<EMAIL_ADDRESS>Few-Body Systems in Physics Laboratory,
RIKEN Nishina Center, Wako 351-0198, Japan  Research Center for Nuclear
Physics (RCNP), Osaka University, Ibaraki, Osaka 567-0047, Japan  Asia Pacific
Center for Theoretical Physics (APCTP), Pohang, Gyeongbuk 37673, Korea  Ho-
Meoyng Choi<EMAIL_ADDRESS>Department of Physics Education, Teachers
College, Kyungpook National University, Daegu 41566, Korea  Chueng-Ryong Ji
<EMAIL_ADDRESS>Department of Physics, North Carolina State University, Raleigh,
NC 27695-8202, USA

###### Abstract

In the light-front quark model (LFQM) amenable to the simultaneous study of
both the mass spectroscopy and the wave function related observables, we
examine the decay constants and distribution amplitudes (DAs) up to the
twist-4. The analysis of the heavy pseudoscalar mesons is carried out both in
the $1S$ and $2S$ states. This investigation involves calculating the local
and nonlocal matrix elements $\langle 0|{\bar{q}}{\Gamma}q|P\rangle$ using
three distinct current operators
${\Gamma}=(\gamma^{\mu}\gamma_{5},i\gamma_{5},\sigma^{\mu\nu}\gamma_{5})$.
Considering a general reference frame where ${\bf P}_{\perp}\neq 0$ and
investigating all available current components, we examine not only the frame-
independence but also the component-independence of the decay constants. The
explicit findings from our study provide the evidence for the equality of the
three pseudoscalar meson decay constants obtained from the three distinct
current operators $\Gamma$. The notable agreement in decay constants is
attained by imposing the Bakamjian-Thomas construction of the LFQM, namely the
meson state is constructed by the noninteracting quark and antiquark
representations while the interaction is added to the mass operator, which
provides the self-consistency condition replacing the physical mass $M$ with
the invariant mass $M_{0}$ for the noninteracting quark-antiquark
representation of the meson state. In addition to obtaining the process-
independent pseudoscalar meson decay constant, regardless of the choice of
current operators $\Gamma$, we further demonstrate its explicit Lorentz and
rotation invariance. In particular, we present the analysis conducted on the
twist-4 DA derived from the minus component of the axial-vector current.
Finally, we discuss the various twist DAs and their $\xi$-moments associated
with the $1S$ and $2S$ heavy pseudoscalar mesons.

##  I Introduction

The distribution amplitudes (DAs) of mesons are important non-perturbative
ingredients in comprehending a range of the light-cone dominated processes
that can be treated via collinear factorization Lepage:1980fj ; ER80 ;
Chernyak:1983ej , as they offer valuable insights into the nonperturbative
makeup of hadrons and the distribution of partons in relation to their
longitudinal momentum fractions within these particles. The meson’s DA is
typically defined as a matrix element of a quark-antiquark bilocal light-cone
operator between the vacuum and the meson state in the light-front dynamics
(LFD) BPP  which provides a natural separation of the meson’s momentum into
its longitudinal and transverse components. Thus, the LFD appears to be a
practical and rigorous framework for computing the DAs of mesons categorizing
them according to their increasing twist.

While the leading-twist DA describes the longitudinal momentum distribution of
valence quarks inside the meson providing a straightforward interpretation of
the partonic structure of the meson, the higher-twist DAs are considerably
more abundant as they take into account of various effects including the
transverse motion of quarks or antiquarks, the higher Fock states involving
extra gluons and/or quark-antiquark pairs, etc. Ball:2006wn . In the light-
cone dominated hard processes, the leading-twist DAs give the dominant
contributions and the higher twist contributions are suppressed by a power of
the hard scale. As a result, the study of higher twist DAs has received less
attention in comparison to the leading twist DAs in the analyses of the hard
processes. However, with the higher statistical precision of experimental data
expected from KEKII, LHC, JLAB, and the forthcoming Electron-Ion-Collider
(EIC) Accardi:2012qut ; AbdulKhalek:2021gbh , the relevance of the higher
twist effects in hadron structure increases, accentuating the growing
importance of further exploring higher twist contributions, e.g. in the
formalism of TMD factorization Braun:2022gzl , in the wide-angle
photoproduction of pions PK18 , and in the nonleptonic $B$-meson decay BBNS
etc. Thus, the quest to obtain essential nonperturbative insights into QCD has
spurred numerous theoretical investigations aimed at calculating not only the
leading-twist DA but also the higher-twist DAs using various nonperturbative
techniques, such as the QCD sum rule Ball:2006wn ; Ball90 ; HWZ04 ; Agaev05 ;
BMS06 ; MPS10 , the chiral-quark model derived from the instanton vacuum PPRWG
; NK06 , the Nambu-Jona-Lasinio (NJL) model AB02 ; PR01 , the Dyson-Schwinger
equation (DSE) approach CRS13 ; SCCRSZ , and the light-front quark model
(LFQM) Hwang10 ; CJ07 ; C13 ; C15 ; Choi:2017uos .

In particular, the LFQM is the theoretical framework based on the LFD that has
been highly successful in explaining simultaneously both the mass spectra and
the wave function related observables including the electroweak properties of
mesons CJ99 ; CJ99B ; HJZR ; Dhiman ; Arifi22 ; Jaus90 ; Jaus91 ; CCH97 ;
Choi07 . In this model, mesons are treated as bound states of constituent
quarks and antiquarks. The LFQM typically places the constituent quark and
antiquark on their respective mass shells, and the spin-orbit (SO) wave
function is obtained through the Melosh transformation Mel1  which is
independent of the interaction and is uniquely determined from the ordinary
quantum state representation $J^{PC}$. For the construction of the more
phenomenologically accessible LFQM, we applied the variational principle with
the trial radial wave function, typically the Gaussian radial wave function,
with the Melosh transformed spin-orbit wave functions for the on-mass shell
constituent quark and antiquark to provide simultaneous analyses of both the
mass spectra and the wave function related observables such as the decay
constants, form factors, etc. CJ99 ; CJ99B ; HJZR ; Dhiman ; Arifi22 ; CCH97 ;
Choi07 . Our LFQM follows the Bakamjian-Thomas (BT) construction BT53 ; KP91 ,
where the meson state is constructed by the noninteracting quark and antiquark
representations while the interaction is added to the mass operator via
$M=M_{0}+V_{Q{\bar{Q}}}$ applying the variational principle with typical
Gaussian radial wave function as the trial wave function for the variational
calculation. Due to the absence of manifest covariance, however, it is
challenging to identify the LF zero-mode contributions in the phenomenological
LFQM by itself. In particular, as the number of DAs proliferates with the
higher twist, the computation of the higher twist DAs involves not only the
good component (e.g. $J^{+}=J^{0}+J^{3}$) of the current but also the bad
component (e.g. $J^{-}=J^{0}-J^{3}$) of the current to identify the
proliferated number of DAs with the more number of current components.
However, employing the bad component in the computation is often quite
challenging due to the involvement of the light-front (LF) zero modes and/or
the instantaneous contributions to restore the Lorentz and gauge invariance
Keister ; deMelo:1997cb ; CJ98 ; DSHwang ; Jaus99 ; deMelo:02 ; BCJ01 . Thus,
it is important to conduct a rigorous study of the higher twist DAs and
address the challenges involved in order to gain the better understanding of
the hadron structures.

To pin down the treacherous points involving the LF zero modes and the
instantaneous contributions, one may utilize the manifestly covariant Bethe-
Salpter (BS) field-theoretic model. Using the LF projection of the manifestly
covariant BS model, one can provide the corresponding LFQM with the multipole
ansatz for the meson-quark vertex function deMelo:1997cb ; deMelo:02 ; BCJ01 .
This type of LFQM is useful for providing the theoretical guidance on how to
analyze the LF zero modes although the obtained results are in general semi-
realistic. To account for the Lorentz structure of a hadronic matrix element,
the light-like four-vector $\omega$ was introduced and the corresponding
covariant approach was developed originally in Ref.Kar76 111In this
formulation, the state vector is defined on a plane characterized by the
invariant equation $\omega\cdot x=0$, where $\omega$ represents an arbitrary
light-like four vector $\omega=(\omega_{0},{\vec{\omega}})$ satisfying
$\omega^{2}=0$. The special choice $\omega=(1,0,0,-1)$ corresponds to the LF
or null plane $\omega\cdot x=x^{+}=x^{0}+x^{3}=0$.. Subsequently, the authors
of Ref. Carbonell98  developed a method for identifying and separating
spurious contributions, enabling them to determine the physically meaningful
contributions to the hadronic form factors and coupling constants that are
independent of the choice of $\omega$. By employing this covariant methodology
as described in Refs. Carbonell98 ; Kar76 , Jaus Jaus99  employed a manifestly
covariant BS model as a reference to devise a fundamentally distinct technique
for addressing this issue. In this approach, Jaus developed a way of
identifying the LF zero-mode contributions by removing spurious contributions
proportional to the light-like vector $\omega^{\mu}$ in the physical
observables. As the $\omega$-dependent contributions violate the covariance,
they may be eliminated by including the LF zero-mode contributions and at the
same time restoring the covariance of the current matrix elements in the
solvable BS model. Jaus identified the light-front zero-mode contribution
corresponding it to the removal of the spurious $\omega$-dependence and then
applied the LF zero-mode contributions identified in the BS model to the LFQM
simply by replacing the multipole type vertex function in the BS model with
the Gaussian radial wave function.

However, two of us C13  found that Jaus’s prescription to identify the zero
mode is only valid in the BS model with the multipole type vertex function but
not in the LFQM with the Gaussian radial wave function as it causes a serious
problem of impeding the self-consistency in the computation of physical
observables, e.g. the decay constant of a $\rho$ meson gives different results
for different polarization (longitudinal and transverse) states in the LFQM.
This finding has also been confirmed by others Chang:2018zjq ; Chang22 ,
indicating that the LF zero-mode contributions do depend on the model wave
functions. In Ref. C13 , we then identified a specific matching condition for
the first time between the manifestly covariant BS model and our LFQM, which
we called the “Type II” link (see Eq. (49) in C13 ). This unique matching
condition ensures the self-consistency of the LFQM analysis. For example, it
was demonstrated C13  by using the “Type II” link that the two $\rho$ meson
decay constants obtained from the longitudinal and transverse polarizations
exhibit the equivalence in the LFQM numerical results. One of the key
ingredients in the “Type II” link is the replacement of the physical meson
mass $M$ that appeared in the integrand for the matrix element calculation
with the invariant mass $M_{0}$, which is equivalent to imposing the on-mass
shell condition of the constituent quark and antiquark in the LFQM. Enforcing
the on-mass shell condition for the constituents is tantamount to ensuring
four-momentum conservation $P=p_{1}+p_{2}$ at the meson-quark vertex, where
the meson and quark (antiquark) momenta are denoted as $P$ and $p_{1(2)}$,
respectively. Such a replacement ($M\to M_{0}$) is indeed consistent with the
BT construction BT53 ; KP91 , namely, the meson state is constructed by the
noninteracting quark and antiquark representations while the interaction is
added to the mass operator via $M=M_{0}+V_{Q{\bar{Q}}}$. This replacement can
also be viewed as an effective way to include the LF zero-mode contributions
and restore the Lorentz symmetry of the model, in particular the rotational
symmetry, when compared to the covariant BS model C13 .

Subsequent works using the same “Type II” link C13  have been made for the
analyses of twist-2 and twist-3 DAs of the light pseudoscalar $(\pi,K)$ mesons
C15 ; Choi:2017uos  through the matrix elements $\langle
0|{\bar{q}}(z){\Gamma}q(-z)|P\rangle$ of the nonlocal operators
${\Gamma}=(\gamma^{\mu}\gamma_{5},i\gamma_{5},\sigma^{\mu\nu}\gamma_{5})$,
discussing the link between the chiral symmetry of QCD and the numerical
results of the LFQM. In the very recent work Arifi:2022qnd , the decay
constant for pseudoscalar meson with the axial-vector
$(\gamma^{\mu}\gamma_{5})$ current in the equal quark and antiquark mass case
was investigated in the LFQM and the self-consistent result independent of the
current components was obtained. The decay constant of the vector meson was
also investigated for both longitudinal and transverse polarizations,
obtaining a self-consistent result independent of all possible combinations of
the current components and the polarizations. In particular, in Ref
Arifi:2022qnd , it was explicitly demonstrated that the decay constants
obtained via the “Type II” link between the BS model and the LFQM are
precisely equivalent to those obtained directly in the LFQM, where the on-mass
shell condition of the constituents is enforced.

In this work, we extend our previous LFQM analyses C15 ; Choi:2017uos  for the
decay constant and the DAs of the ground $1S$ state light pseudoscalar mesons
through the matrix elements $\langle 0|{\bar{q}}(z){\Gamma}q(-z)|P\rangle$ of
the operators
${\Gamma}=(\gamma^{\mu}\gamma_{5},i\gamma_{5},\sigma^{\mu\nu}\gamma_{5})$ to
include both $(1S,2S)$ state heavy-light and heavy-heavy pseudoscalar mesons
Arifi:2022qnd . In particular, we shall explicitly show that the three
pseudoscalar meson decay constants defined through the three different
operators $\Gamma$ are all identical numerically in our LFQM constrained by
the on-mass condition of the constituents. Namely, we obtain the process-
independent pseudoscalar meson decay constant regardless of the current
operators $\Gamma$ used at the level of one-body current matrix element
computation, as the independence of the current operators $\Gamma$ means the
independence of the decay process. The new two-particle twist-4 DA is also
obtained from the minus component of the axial vector current
($\Gamma=\gamma^{\mu}\gamma_{5}$). We also investigate the different helicity
contributions to the decay constants defined through different operators
$\Gamma$ and perform a quantitative analysis of each helicity component for
different heavy-light and heavy-heavy pseudoscalar meson systems. For the
numerical calculations, we present the results both for the ground state
($1S$) and the radially excited state ($2S$) of heavy pseudoscalar mesons,
which were discussed in our recent work Arifi22 . We then scrutinize the shape
of the leading- and higher-twist DAs and their $\xi$-moments, where $\xi=2x-1$
with the LF longitudinal momentum fraction $x$ of the constituent. The “Type
II” link between the covariant BS model and the LFQM is further discussed for
the deeper understanding of the underlying physics involved.

The paper is organized as follows: In Section II, we describe the LFQM and the
light-front wave functions of $1S$ and $2S$ pseudoscalar heavy meson. In
Section III, we examine the pseudoscalar decay constants derived from the
three distinct current operators
$\Gamma=(\gamma^{\mu}\gamma_{5},i\gamma_{5},\sigma^{\mu\nu}\gamma_{5})$ and
establish the process independence and rotational invariance of the decay
constants. In Section IV, we discuss the DAs up to the twist-4 obtained from
the three local and nonlocal current operators $\Gamma$ and their
$\xi$-moments. Finally, we summarize our findings in Section V. In the
Appendix, the “Type II” link between the manifestly covariant BS model and the
LFQM is demonstrated for the completeness.

##  II Light-front quark model

When applied to meson states reflecting the feature of BT construction BT53 ;
KP91 , the LFQM employs a noninteracting $q{\bar{q}}$ representation to
describe the Fock state that is composed of the constituent quark ($q$) and
antiquark (${\bar{q}})$ while the interactions are incorporated into the mass
operator $M:=M_{0}+V_{q{\bar{q}}}$ to ensure compliance with the group
structure satisfying the Poincaré algebraic commutation relations. The
interactions are then encoded in the light-front wave function (LFWF)
$\Psi_{q{\bar{q}}}$, which satisfies the eigenvalue equation
$H_{q\bar{q}}|\Psi_{q{\bar{q}}}\rangle=(M_{0}+V_{q\bar{q}})|\Psi_{q{\bar{q}}}\rangle=M_{q\bar{q}}|\Psi\rangle$
of our QCD-motivated effective Hamiltonian CJ99 ; CJ99B ; HJZR ; Dhiman ;
Arifi22 .

Our LFQM for the $1S$ state CJ99 ; CJ99B ; HJZR ; Dhiman ; Arifi22  and $2S$
state Arifi22  pseudoscalar and vector mesons is based on the central concept
of using the radial wave function as a variational trial function for the QCD-
motivated effective Hamiltonian $H_{q{\bar{q}}}$, which results in the
determination of the mass eigenvalues $M_{q{\bar{q}}}$. Once the values of the
model parameters are determined by the variational analysis of the mass
spectra, those determined model parameters are used to describe different
observables including decay constants and electromagnetic and weak form
factors etc. CJ99 ; CJ99B ; HJZR ; Dhiman ; Arifi22 ; CJ07 ; C21

For the self-consistent analysis of the decay constants and the higher twist
DAs for the $(1S,2S)$ state heavy pseudoscalar mesons performed in this work,
we provide a brief overview of the LFWFs for $1S$ and $2S$ state heavy
pseudoscalar mesons presented in Ref. Arifi22  focusing on the important
aspects of LFWFs constrained by the on-mass shell condition of the
constituents.

The four-momentum $P$ of the meson in terms of the LF components is defined as
$P=(P^{+},P^{-},{\bf P}_{\perp})$, where $P^{+}=P^{0}+P^{3}$ and
$P^{-}=P^{0}-P^{3}$ are the LF longitudinal momentum and the LF energy,
respectively, and ${\bf P}_{\perp}=(P^{1},P^{2})$ are the transverse momenta.
Here, we take the metric convention as $P^{2}=P^{+}P^{-}-{\bf P}^{2}_{\perp}$.
The meson state $|{\rm M}(P,J,J_{z})\rangle$ of momentum $P$ and spin state
$(J,J_{z})$ can then be constructed as follows CCH97 ; Cheng04 ; Cheng:2004ew

$\displaystyle\ket{\rm M}$ $\displaystyle=$ $\displaystyle\int\left[{\rm
d}^{3}{\bar{p}}_{1}\right]\left[{\rm
d}^{3}{\bar{p}}_{2}\right]2(2\pi)^{3}\delta^{3}\left({\bar{P}}-{\bar{p}}_{1}-{\bar{p}}_{2}\right)$
(1)
$\displaystyle\times\mbox{}\sum_{\lambda_{1},\lambda_{2}}\Psi_{\lambda_{1}\lambda_{2}}^{JJ_{z}}(x,{\bf
k}_{\perp})\ket{q_{\lambda_{1}}(p_{1})\bar{q}_{\lambda_{2}}(p_{2})},\quad$

where $p^{\mu}_{i}$ and $\lambda_{i}$ are the on-mass shell
($p^{2}_{i}=m^{2}_{i}$) momenta and the helicities of the constituent quark
($i=1$) and antiquark $(i=2)$, respectively, with the LF three momentum
defined by ${\bar{p}}=(p^{+},{\bf p}_{\perp})$ and $\left[{\rm
d}^{3}{\bar{p}}\right]\equiv{\rm d}p^{+}{\rm
d}^{2}\mathbf{p}_{\perp}/(16\pi^{3})$. The LF internal relative variables
$(x,{\bf k}_{\perp})$ are defined by $x_{i}=p^{+}_{i}/P^{+}$ and ${\bf
k}_{i\perp}={\bf p}_{i\perp}-x_{i}{\bf P}_{\perp}$, where $\sum_{i}x_{i}=1$
and $\sum_{i}{\bf k}_{i\perp}=0$ and we set $x=x_{1}$ and ${\bf
k}_{\perp}={\bf k}_{1\perp}$. This meson state satisfies the following
normalization

$\displaystyle\langle{{\rm M}(P^{\prime},J^{\prime},J^{\prime}_{z})}|{\rm
M}(P,J,J_{z})\rangle$
$\displaystyle=2(2\pi)^{3}P^{+}\delta^{3}({\bar{P}}^{\prime}-{\bar{P}})\delta_{J^{\prime}J}\delta_{J^{\prime}_{z}J_{z}}.$
(2)

In momentum space, the LFWF $\Psi_{q{\bar{q}}}(x,{\bf k}_{\perp})$ of a meson
can be decomposed as

$\displaystyle\Psi^{JJ_{z}}_{\lambda_{1}\lambda_{2}}(x,\mathbf{k}_{\bot})=\Phi(x,\mathbf{k}_{\bot})\
\mathcal{R}^{JJ_{z}}_{\lambda_{1}\lambda_{2}}(x,\mathbf{k}_{\bot}),$ (3)

where $\Phi(x,\mathbf{k}_{\bot})$ is the radial wave function that was used as
our trial function for the mass spectroscopic analysis CJ99 ; CJ99B ; HJZR ;
Dhiman ; Arifi22  and $\mathcal{R}^{JJ_{z}}_{\lambda_{1}\lambda_{2}}$ is the
SO wave function obtained by the interaction-independent Melosh transformation
Mel1  for the corresponding meson quantum number $J^{PC}$.

We should note that one crucial aspect of the LF formulation for a bound
state, as depicted in Eq. (1), is the frame-independence of the LFWF BPP . In
other words, the hadron’s internal variables $(x,{\bf k}_{\perp})$ of the wave
function remain unaffected by boosts to any physical $(P^{+},{\bf P}_{\perp})$
frame, which is not the case in the instant formulation. We shall explicitly
show this boost invariance of the decay constant computed in general ${\bf
P}_{\perp}\neq 0$ frame.

The SO wave function for a pseudoscalar meson obtained from the interaction-
independent Melosh transformation can be written as the covariant form Jaus90
; Jaus91  consistent with the BT construction, which is given by

$\displaystyle\mathcal{R}^{00}_{\lambda_{1}\lambda_{2}}$ $\displaystyle=$
$\displaystyle\frac{\bar{u}_{\lambda_{1}}(p_{1})\gamma_{5}v_{\lambda_{2}}(p_{2})}{\sqrt{2}{\tilde{M}_{0}}},$
(4) $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{2}\sqrt{\mathcal{A}^{2}+\mathbf{k}_{\perp}^{2}}}\begin{pmatrix}-k^{L}&\mathcal{A}\\\
-\mathcal{A}&-k^{R}\\\ \end{pmatrix},$ (5)

where
${\tilde{M}_{0}}^{2}=M_{0}^{2}-(m_{1}-m_{2})^{2}=(\mathcal{A}^{2}+\mathbf{k}_{\perp}^{2})/x(1-x)$,
$k^{R(L)}=k^{1}\pm ik^{2}$ and $\mathcal{A}=(1-x)m_{1}+xm_{2}$. The boost-
invariant meson mass squared is given by

$\displaystyle
M_{0}^{2}=\frac{\mathbf{k}_{\bot}^{2}+m_{1}^{2}}{x}+\frac{\mathbf{k}_{\bot}^{2}+m_{2}^{2}}{1-x}.$
(6)

Note that the SO wave function satisfies the unitary condition,
$\sum_{\lambda^{\prime}s}{\cal R}^{\dagger}{\cal R}=1$. It is worth noting
that the LFWF $\Psi$ depends on the interaction-independent invariant mass
$M_{0}$ that follows the BT construction as the meson is constructed in the
noninteracting representation.

For the $1S$ and $2S$ state radial wave functions $\Phi_{ns}$ of Eq. (3), we
allow the mixing between the two lowest order harmonic oscillator (HO) wave
functions $(\phi_{1S},\phi_{2S})$ by writing Arifi22

$\displaystyle\begin{pmatrix}\Phi_{1S}\\\ \
\Phi_{2S}\end{pmatrix}=\begin{pmatrix}\cos\theta&\sin\theta\\\
-\sin\theta&\cos\theta\end{pmatrix}\begin{pmatrix}{\phi}_{1S}\\\
{\phi}_{2S}\end{pmatrix},$ (7)

where

$\displaystyle\phi_{1S}({\vec{k}})$ $\displaystyle=$
$\displaystyle\frac{4\pi^{3/4}}{\beta^{3/2}}e^{-{\vec{k}}^{2}/2\beta^{2}},$
$\displaystyle\phi_{2S}(\vec{k})$ $\displaystyle=$
$\displaystyle\frac{4\pi^{3/4}}{\sqrt{6}\beta^{7/2}}\left(2{\vec{k}}^{2}-3\beta^{2}\right)e^{-{\vec{k}}^{2}/2\beta^{2}}.$
(8)

Here, ${\vec{k}}=(k_{z},\mathbf{k}_{\perp})$ is the three momentum and $\beta$
represents a parameter that serves as the variational parameter in our mass
spectroscopic analysis Arifi22 . The rotationally invariant HO wave functions
$\phi_{nS}(\vec{k})$ in Eq. (II) satisfy

$\int\;\frac{{\rm
d}^{3}{\vec{k}}}{2(2\pi)^{3}}\;\absolutevalue{\phi_{nS}({\vec{k}})}^{2}=1.$
(9)

Sustaining this rotationally invariant property of the wave function, one can
transform the normalization of $\phi_{nS}({\vec{k}})$ to that of
$\phi_{nS}(x,{\bf k}_{\perp})$ via the variable transformation
$(k_{z},\mathbf{k}_{\perp})\to(x,\mathbf{k}_{\perp})$ as follows

$\int_{0}^{1}{\rm d}x\int\frac{{\rm
d}^{2}\mathbf{k}_{\bot}}{2(2\pi)^{3}}\absolutevalue{\phi_{nS}(x,\mathbf{k}_{\bot})}^{2}=1.$
(10)

We note that the wave functions $\phi_{nS}(x,\mathbf{k}_{\bot})$ include the
Jacobian factor $\partial k_{z}/\partial x$ as

$\displaystyle\phi_{nS}(x,\mathbf{k}_{\perp})=\sqrt{\frac{\partial
k_{z}}{\partial x}}\phi_{nS}(\vec{k})$ (11)

because of the variable transformation
$(k_{z},\mathbf{k}_{\perp})\to(x,\mathbf{k}_{\perp})$ and $k_{z}(=k^{3})$ and
$x$ are related by CCH97

$x=\frac{E_{1}-k_{z}}{E_{1}+E_{2}},\;1-x=\frac{E_{2}+k_{z}}{E_{1}+E_{2}},$
(12)

where $E_{i}=\sqrt{m^{2}_{i}+{\vec{k}}^{2}}$. We then have $M_{0}=E_{1}+E_{2}$
and

$\displaystyle
k_{z}=\left(x-\frac{1}{2}\right)M_{0}+\frac{m^{2}_{1}-m^{2}_{2}}{2M_{0}}.$
(13)

The Jacobian factor is then given by

$\frac{\partial k_{z}}{\partial x}=\frac{E_{1}E_{2}}{x(1-x)M_{0}},$ (14)

or in terms of $(x,{\bf k})$

$\frac{\partial k_{z}}{\partial
x}=\frac{M_{0}}{4x(1-x)}\left[1-\frac{(m_{1}^{2}-m_{2}^{2})^{2}}{M_{0}^{4}}\right].$
(15)

It should be noted that the total LFWF $\Psi$ given by Eq. (3) meets the same
normalization given by Eq. (10). This is due to the meson state $|{\rm
M}(P,J,J_{z})\rangle$ fulfilling the condition of Eq. (II), and the SO wave
function adhering to the unitary condition. Especially, the inclusion of the
Jacobian factor in defining $\phi_{nS}(x,\mathbf{k}_{\perp})$ is the key
aspect to retaining the rotational invariance of the model wave function and
obtaining the self-consistent, i.e. current-component and boost invariant,
results of physical observables. The quantitative effects of the Jacobian
factor on the decay constants, twist-2 DAs, and electromagnetic form factors
of mesons were also discussed in Ref. CJM0 .

##  III Decay constants

For the decay constants and the leading-and higher-twist DAs of pseudoscalar
mesons, one may obtain them from the matrix elements $\langle
0|{\bar{q}}(z){\Gamma}q(-z)|P\rangle$ of the the following three possible
nonlocal operators
${\Gamma}=(\gamma^{\mu}\gamma_{5},i\gamma_{5},\sigma^{\mu\nu}\gamma_{5})$,
where $z^{\mu}$ is the light-like vector ($z^{2}=0$). For the calculation of
the decay constant $f_{\rm P}$, while $f_{\rm P}$ can be defined by using a
local operator with axial-vector $(\Gamma_{\rm A}=\gamma^{\mu}\gamma_{5})$ and
pseudoscalar $(\Gamma_{\rm P}=i\gamma_{5})$ current as Ball:2006wn ; Ball90

$\displaystyle\bra{0}\bar{q}(0)\gamma^{\mu}\gamma_{5}q(0)\ket{{\rm P}(P)}$
$\displaystyle=$ $\displaystyle if_{\rm P}P^{\mu},$ (16)
$\displaystyle\bra{0}\bar{q}(0)i\gamma_{5}q(0)\ket{{\rm P}(P)}$
$\displaystyle=$ $\displaystyle f_{\rm P}\mu_{M},$ (17)

it can also be computed by utilizing the nonlocal matrix element in the case
of the pseudotensor current $(\Gamma_{\rm T}=\sigma^{\mu\nu}\gamma_{5})$ as
defined by the subsequent equation Ball90 :

$\displaystyle\bra{0}\bar{q}(z)\sigma^{\mu\nu}\gamma_{5}q(-z)\ket{{\rm
P}(P)}=-\frac{i}{3}f_{\rm P}\left(1-\rho_{+}\right)\mu_{M}$
$\displaystyle\times(P^{\mu}z^{\nu}-P^{\nu}z^{\mu})\int_{0}^{1}{\rm d}x\ {\rm
e}^{i\zeta P\cdot z}\psi_{3;{\rm P}}(x),\quad\quad$ (18)

where $\sigma^{\mu\nu}=\frac{i}{2}[\gamma^{\mu},\gamma^{\nu}]$,
$\mu_{M}=M^{2}/(m_{1}+m_{2})$, $\rho_{+}=(m_{1}+m_{2})^{2}/M^{2}$, and
$\zeta=2x-1$. In this definition of Eq. (III), the two-particle twist 3 DA
$\psi_{3;{\rm P}}(x)$ is normalized to unity $\int^{1}_{0}{\rm d}x\psi_{3;{\rm
P}}(x)=1$. As a reference, in Ref. Ball:2006wn , defining the matrix element
$\bra{0}\bar{q}(z)\sigma^{\mu\nu}\gamma_{5}q(-z)\ket{{\rm P}(P)}$, the authors
removed the term $(1-\rho_{+})$ on the right-hand side (RHS) of Eq. (III) by
normalizing $\psi_{3;{\rm P}}(x)$ in such a way that $\int^{1}_{0}{\rm
d}x\;\psi_{3;{\rm P}}(x)=1-\rho_{+}$. In the previous LFQM analysis
Choi:2017uos  for $\psi_{3;{\rm P}}(x)$, two of us used the definition of Ref.
Ball:2006wn  rather than Eq. (III). However, in this study, we opt to use Eq.
(III) as we observe that this definition yields the same decay constant as
those obtained from Eqs. (16) and (17) in which the same normalization for the
leading and higher twist DAs is used as $\psi_{3;{\rm P}}(x)$ defined in Eq.
(III).

###  III.1 Process-Independence

In this subsection, we shall first compute the decay constants defined by the
three different operators $\Gamma=(\Gamma_{\rm A},\Gamma_{\rm P},\Gamma_{\rm
T})$ and show their equivalence, i.e. process-independent decay constant in
our LFQM. As the different decay operators are used for the different decay
processes, the decay constant’s independence of the current operators $\Gamma$
means the independence of the decay process for the decay constant as a
physical observable. The leading-and higher-twist DAs obtained from the matrix
elements $\langle 0|{\bar{q}}(z){\Gamma}q(-z)|P\rangle$ will be analyzed
separately in the next section.

In the LFQM, the decay amplitudes for the operators
$\Gamma=(\Gamma_{A},\Gamma_{\rm P})$ given by Eqs. (16) and (17) can be
defined at the level of one-body local current matrix element as

$\displaystyle\bra{0}\bar{q}\Gamma q\ket{P}$ $\displaystyle=$
$\displaystyle\sqrt{N_{c}}\int_{0}^{1}{\rm d}x\int\frac{{\rm
d}^{2}\mathbf{k}_{\bot}}{16\pi^{3}}\ \Phi(x,\mathbf{k}_{\perp})$ (19)
$\displaystyle\times\mbox{}\sum_{\lambda_{1},\lambda_{2}}{\cal
R}_{\lambda_{1}\lambda_{2}}^{00}\left[\frac{\bar{v}_{\lambda_{2}}(p_{2})}{\sqrt{x_{2}}}\Gamma\frac{u_{\lambda_{1}}(p_{1})}{\sqrt{x_{1}}}\right],\quad$

where $N_{c}=3$ arises from the color factor implicit in the wave function
Jaus91 ; CCH97 . Denoting the decay constants $f_{\rm P}$ corresponding to the
current operators $(\Gamma_{\rm A},\Gamma_{\rm P},\Gamma_{\rm T})$ as $(f_{\rm
A},f_{\rm P},f_{\rm T})$, we may provide the generic form for the decay
constants $f_{\rm A(P)}$ obtained from the two local operators $\Gamma_{\rm
A(P)}$ as Arifi:2022qnd

$\displaystyle f_{\rm A(P)}$ $\displaystyle=$
$\displaystyle\sqrt{N_{c}}\int_{0}^{1}{\rm d}x\int\frac{{\rm
d}^{2}\mathbf{k}_{\bot}}{16\pi^{3}}\ \Phi(x,\mathbf{k}_{\perp})$
$\displaystyle\times\mbox{}\frac{1}{{\cal P}_{\rm
A(P)}}\sum_{\lambda_{1},\lambda_{2}}{\cal
R}_{\lambda_{1}\lambda_{2}}^{00}\left[\frac{\bar{v}_{\lambda_{2}}(p_{2})}{\sqrt{x_{2}}}\Gamma_{\rm
A(P)}\frac{u_{\lambda_{1}}(p_{1})}{\sqrt{x_{1}}}\right],$

where we incorporate the Lorentz structures ${\cal P}_{\rm
A(P)}=iP^{\mu}(\mu_{M})$ on the RHS of Eqs. (16) and (17) into the integral.
This incorporation of Lorentz structures into the integral is to assure the
consistent one-body current level of approximation in the computation of the
decay constant, which is the crucial aspect of our recently developed LFQM
analysis C13 ; C15 ; Choi:2017uos ; Arifi:2022qnd ; C21  to obtain the self-
consistent, i.e. current-component and boost invariant as well as the process
(e.g. $\Gamma_{\rm A(P)}$) independent physical observables by replacing all
physical mass $M$ appeared in Eq. (III.1) with the invariant mass $M_{0}$. So
far, most LF calculations of the decay constant, e.g. $f_{\rm A}$ from
$\Gamma_{\rm A}$, used $\mu=+$ or $\perp$ since $(P^{+},{\bf P}_{\perp})$ do
not involve any physical mass $M$. On the other hand, the minus component of
the axial-vector current involves $P^{-}=(M^{2}+{\bf P}^{2}_{\perp})/P^{+}$
and one fails to produce the same result as the one obtained from the currents
with $(\mu=+,\perp)$ if one uses the physical mass in the calculation. The
difference between $\mu=-$ and $\mu=(+,\perp)$ was identified as the LF
treacherous points such as the instantaneous and zero-contributions to the
minus current in the solvable covariant model C13 . However, in our LFQM
consistent with BT construction, we showed Arifi:2022qnd ; C21  that the
result from the minus component of axial-vector current gives the same result
as the one obtained from the currents with $(\mu=+,\perp)$ if we replace $M$
with $M_{0}$. In this work, we shall show that the result obtained from
$\Gamma_{\rm P}$ also gives the same result as the one obtained from
$\Gamma_{\rm A}$ as far as we use $M\to M_{0}$ prescription. This may be
regarded as the effective LF zero-mode inclusion at the level of the one-body
matrix element calculation in the LFQM is consistent with BT construction.

For the pseudotensor current ($\Gamma^{\mu\nu}_{\rm
T}=\sigma^{\mu\nu}\gamma_{5}$), since the decay constant can be computed only
in the nonlocal limit (i.e. $z^{\mu}\neq 0$), the calculation incorporating
$\psi_{3;{\rm P}}(x)$ is inevitably required in the process of deriving the
decay constant from the pseudotensor current. From the light-like vector
$z^{\mu}$ with $z^{+}={\bf z}_{\perp}=0$, there are two possible ways to
compute the nonlocal matrix element by choosing $\mu\nu=+-$ or $\perp-$.

Table 1: Various helicity contributions $H_{\lambda_{1}\lambda_{2}}$ to the current operators $\mathcal{O}$ for all possible components of the currents, where $\mathcal{A}=(1-x)m_{1}+xm_{2}$, $\mathcal{B}=(1-x)m_{1}-xm_{2}$, $\Delta_{1}=(m_{1}^{2}+\mathbf{k}_{\perp}^{2})/x$, $\Delta_{2}=(m_{2}^{2}+\mathbf{k}_{\perp}^{2})/(1-x)$, $\tilde{M}_{0}^{2}={(\mathbf{k}_{\perp}^{2}+\mathcal{A}^{2})}/{x(1-x)}$, $\mu^{0}_{M}=M_{0}^{2}/(m_{1}+m_{2})$, and $\rho^{0}_{+}=(m_{1}+m_{2})^{2}/M_{0}^{2}$. We should note that the primed symbols in pseudotensor results represent functions of $x^{\prime}$ rather than $x$. $\mathcal{O}$ | Current | $H_{\uparrow\uparrow}+H_{\downarrow\downarrow}$ | $H_{\uparrow\downarrow}+H_{\downarrow\uparrow}$ | $\mathcal{O}=\sum H_{\lambda_{1}\lambda_{2}}$
---|---|---|---|---
$\mathcal{O}_{\rm A}^{\mu}$ | $\Gamma^{+}_{\rm A},\Gamma^{\perp}_{\rm A}$ | $0$ | 2$\mathcal{A}$ | 2$\mathcal{A}$
$\Gamma^{-}_{\rm A}$ | $\dfrac{2(m_{1}+m_{2})\mathbf{k}_{\perp}^{2}}{x(1-x)(M_{0}^{2}+\mathbf{P}_{\perp}^{2})}$ | $\dfrac{2\mathcal{A}[x(1-x)\mathbf{P}_{\perp}^{2}-\mathbf{k}_{\perp}^{2}+m_{1}m_{2}]}{x(1-x)(M_{0}^{2}+\mathbf{P}_{\perp}^{2})}$ | $\dfrac{2(m_{1}\Delta_{2}+m_{2}\Delta_{1}+\mathcal{A}\mathbf{P}^{2}_{\perp})}{(M^{2}_{0}+\mathbf{P}^{2}_{\perp})}$
$\mathcal{O}_{\rm P}$ | $\Gamma_{\rm P}$ | $\dfrac{\mathbf{k}_{\perp}^{2}}{x(1-x)\mu^{0}_{M}}$ | $\dfrac{\mathcal{A}^{2}}{x(1-x)\mu^{0}_{M}}$ | $\dfrac{\tilde{M}_{0}^{2}}{\mu^{0}_{M}}$
$\mathcal{O}_{\rm T}^{\mu\nu}$ | $\Gamma^{+-}_{\rm T},\Gamma^{\perp-}_{\rm T}$ | $\dfrac{(1-2x^{\prime})\mathbf{k}_{\perp}^{2}}{2x^{\prime}(1-x^{\prime})(1-\rho^{\prime 0}_{+})\mu^{\prime 0}_{M}}$ | $\dfrac{\mathcal{A^{\prime}}\mathcal{B^{\prime}}}{2x^{\prime}(1-x^{\prime})(1-\rho^{\prime 0}_{+})\mu^{\prime 0}_{M}}$ | $\dfrac{-12M^{\prime}_{0}k^{\prime}_{z}}{\mu^{\prime 0}_{M}(1-\rho^{\prime 0}_{+})}$

As for an example, let us choose $\mu\nu=+-$. We first integrate Eq. (III) on
both sides using the dummy variable $x^{\prime}$ (and
$\zeta^{\prime}=2x^{\prime}-1$) with respect to $z^{-}$ as

$\displaystyle\int_{-\infty}^{\infty}\frac{{\rm d}z^{-}}{2\pi}{\rm
e}^{-i\zeta^{\prime}P\cdot z}\bra{0}\bar{q}(z)\Gamma^{+-}_{\rm T}q(-z)\ket{P}$
$\displaystyle=CP^{+}\int_{0}^{1}{\rm d}x\ \int_{-\infty}^{\infty}\frac{{\rm
d}z^{-}}{2\pi}z^{-}{\rm e}^{-i(x^{\prime}-x)P^{+}z^{-}}\psi_{3;{\rm
P}}(x),\quad\quad$ (21)

where $C=-\frac{i}{3}f_{\rm T}\left(1-\rho_{+}\right)\mu_{M}$. Then, we obtain
the RHS of Eq. (III.1) via

$\displaystyle\int_{-\infty}^{\infty}\frac{{\rm d}z^{-}}{2\pi}z^{-}{\rm
e}^{-i(x^{\prime}-x)P^{+}z^{-}}\psi_{3;{\rm P}}(x)$
$\displaystyle=\frac{i}{P^{+}}\frac{\partial}{\partial
x^{\prime}}\int_{-\infty}^{\infty}\frac{{\rm d}z^{-}}{2\pi}{\rm
e}^{-i(x^{\prime}-x)P^{+}z^{-}}\psi_{3;{\rm P}}(x)$
$\displaystyle=\frac{i}{P^{+}}\frac{\partial}{\partial
x^{\prime}}\left[\delta((x^{\prime}-x)P^{+})\psi_{3;{\rm P}}(x)\right]$ (22)

as follows

${\rm RHS\;of\;Eq.~{}\eqref{eq:ten1}}=\frac{1}{3P^{+}}f_{\rm
T}\left(1-\rho_{+}\right)\mu_{M}\frac{\partial}{\partial
x^{\prime}}\psi_{3;{\rm P}}(x^{\prime}).$ (23)

On the other hand, the left-hand side (LHS) of Eq. (III) can be rewritten as

$\displaystyle{\rm LHS\;of\;Eq.~{}\eqref{eq:ten1}}$
$\displaystyle=\sqrt{N_{c}}\int^{1}_{0}{\rm d}x\int\frac{{\rm
d}^{2}\mathbf{k}_{\perp}}{16\pi^{3}}\int_{-\infty}^{\infty}\frac{{\rm
d}z^{-}}{2\pi}{\rm e}^{-i\zeta^{\prime}P\cdot z}{\rm e}^{-i(p_{2}-p_{1})\cdot
z}$
$\displaystyle\times\sum_{\lambda_{1},\lambda_{2}}\Psi_{\lambda_{1}\lambda_{2}}^{00}(x,{\bf
k}_{\perp})\left[\frac{\bar{v}_{\lambda_{2}}(p_{2})}{\sqrt{x_{2}}}\Gamma^{+-}_{\rm
T}\frac{u_{\lambda_{1}}(p_{1})}{\sqrt{x_{1}}}\right],$ (24)

where ${\rm e}^{-i\zeta^{\prime}P\cdot z}{\rm e}^{-i(p_{2}-p_{1})\cdot z}={\rm
e}^{-i(x^{\prime}-x)P^{+}z^{-}}$ and the $z^{-}$ integration gives the
$\delta[(x^{\prime}-x)P^{+}]$ and then it is trivially integrated by ${\rm
d}x$. Integrating Eqs. (23) and (III.1) over $x^{\prime}$, we obtain
$\psi_{3;P}(x)$ as

$\displaystyle\psi_{3;P}(x)$ $\displaystyle=$
$\displaystyle\frac{3\sqrt{N_{c}}}{f_{\rm T}}\int_{0}^{x}{\rm
d}x^{\prime}\int\frac{{\rm
d}^{2}\mathbf{k}_{\perp}}{16\pi^{3}}\Phi(x^{\prime},\mathbf{k}_{\perp})$
$\displaystyle\times\frac{1}{{\cal P}_{\rm
T}}\sum_{\lambda_{1},\lambda_{2}}\mathcal{R^{\prime}}_{\lambda_{1}\lambda_{2}}^{00}\left[\frac{\bar{v}_{\lambda_{2}}(p^{\prime}_{2})}{\sqrt{x_{2}^{\prime}}}\Gamma^{+-}_{\rm
T}\frac{u_{\lambda_{1}}(p^{\prime}_{1})}{\sqrt{x_{1}^{\prime}}}\right],\quad$

where ${\cal P}_{\rm T}=(1-\rho_{+})\mu_{M}$ and the prime($\prime$) in
$({\cal R},p_{i})$ implies that they are functions of $x^{\prime}$. By
integrating both sides with respect to ${\rm d}x$ and using the normalization
of the DA, $\int_{0}^{1}{\rm d}x\ \psi_{3;P}(x)=1$, we obtain the decay
constant $f_{{\rm T}}$ from the pseudotensor channel as

$\displaystyle f_{\rm T}$ $\displaystyle=$ $\displaystyle
3\sqrt{N_{c}}\int_{0}^{1}{\rm d}x\int_{0}^{x}{\rm
d}x^{\prime}\int\frac{d^{2}\mathbf{k}_{\perp}}{16\pi^{3}}\Phi(x^{\prime},\mathbf{k}_{\perp})$
(26) $\displaystyle\times\frac{1}{{\cal P}_{\rm
T}}\sum_{\lambda_{1},\lambda_{2}}\mathcal{R^{\prime}}_{\lambda_{1}\lambda_{2}}^{00}\left[\frac{\bar{v}_{\lambda_{2}}(p^{\prime}_{2})}{\sqrt{x_{2}^{\prime}}}\Gamma^{+-}_{\rm
T}\frac{u_{\lambda_{1}}(p^{\prime}_{1})}{\sqrt{x_{1}^{\prime}}}\right].\quad$

We should note that the term ${\cal P}_{\rm T}$ including the physical mass
$M$ is also incorporated into the integral so that $M$ is replaced with
$M_{0}$. The same results for $\psi_{3;P}(x)$ in Eq. (III.1) and $f_{\rm T}$
in Eq. (26) can be obtained with $\Gamma^{\perp-}_{\rm T}$. We also note that
the main update on the calculation of pseudotensor current compared to Ref.
Choi:2017uos  is the inclusion of the term $(1-\rho_{+})$, which leads to the
process independence of the decay constant, i.e. $f_{\rm A}=f_{\rm P}=f_{\rm
T}$, as we discuss below.

Here, we explicitly demonstrate that all three decay constants, $(f_{\rm
A},f_{\rm P},f_{\rm T})$ as defined by Eqs. (III.1) and (26), yield identical
numerical results. Using the Dirac helicity spinors Lepage:1980fj ; Jaus90
and the SO wave function defined in Eq. (4), it is straightforward to compute
$(f_{\rm A},f_{\rm P},f_{\rm T})$, especially, in terms of different helicity
contributions for different usage of current operators. The final results of
$(f_{\rm A},f_{\rm P},f_{\rm T})$ in the most general ${\bf P}_{\perp}\neq 0$
frame are summarized as follows

$f_{\rm A(P)}=\sqrt{6}\int_{0}^{1}{\rm d}x\int\frac{{\rm
d}^{2}\mathbf{k}_{\bot}}{16\pi^{3}}\
\frac{{\Phi}(x,\mathbf{k}_{\bot})}{\sqrt{{\mathcal{A}}^{2}+\mathbf{k}_{\bot}^{2}}}~{}\mathcal{O}_{\rm
A(P)}(x,{\bf k}_{\perp}),$ (27)

and

$f_{\rm T}=\sqrt{6}\int_{0}^{1}{\rm d}x\int_{0}^{x}{\rm
d}x^{\prime}\int\frac{{\rm
d}^{2}\mathbf{k}_{\perp}}{16\pi^{3}}\frac{\Phi({x^{\prime}},\mathbf{k}_{\perp})}{\sqrt{{\mathcal{A^{\prime}}}^{2}+\mathbf{k}_{\perp}^{2}}}{\mathcal{O}}_{\rm
T}(x^{\prime},{\bf k}_{\perp}),$ (28)

where ${\cal A^{\prime}}={\cal A}(x\to x^{\prime})$. The operators
$\mathcal{O}$ given by Eqs. (27) and (28) are obtained from the sum of each
helicity contribution $H_{\lambda_{1}\lambda_{2}}$, i.e.,

$\displaystyle\mathcal{O}=\sum_{\lambda_{1},\lambda_{2}}H_{\lambda_{1}\lambda_{2}}.$
(29)

The results of each helicity contributions $H_{\lambda_{1}\lambda_{2}}$ and
their sum ${\cal O}$ defined by Eq. (29) for different current operators
$\Gamma=(\Gamma_{\rm A},\Gamma_{\rm P},\Gamma_{\rm T})$ together with
different components of the currents for $\Gamma=(\Gamma_{\rm A},\Gamma_{\rm
T})$ are summarized in Table 1. We should note that all the physical masses
$M$ are replaced with the invariant mass $M_{0}$ in the final results
presented in Table 1. We confirmed that the three decay constants given by
Eqs. (27) and (28) are the same as each other, i.e. the pseudoscalar meson
decay constant in our LFQM can be obtained in the process-independent manner
(i.e. $f_{\rm A}=f_{\rm P}=f_{\rm T}$) regardless of the current operators
$(\Gamma=\Gamma_{\rm A},\Gamma_{\rm P},\Gamma_{\rm T})$ used for the local and
nonlocal matrix elements.

In the Appendix, we also discuss the “Type II” link C13  between the covariant
BS model and our LFQM, which is the alternative method to obtain the self-
consistent LFQM results for the decay constants given by Eqs. (27) and (28).

###  III.2 Lorentz and Rotation Invariance

In this work, we also compute the decay constant with the nonvanishing
$\mathbf{P}_{\perp}$ frame. As one can see from Table 1, the operators ${\cal
O}_{\rm P}$, ${\cal O}^{\mu\nu}_{\rm T}$, and ${\cal O}^{+,\perp}_{\rm A}$ are
completely independent of $\mathbf{P}_{\perp}$. Although the operator ${\cal
O}_{\rm A}$ obtained from the minus component of the current $\Gamma^{-}_{\rm
A}$ depends on ${\bf P}_{\perp}$, which is originated from $P^{-}$ associated
with the Lorentz factor $P^{\mu}$ on the RHS of Eq. (16), we confirm that the
decay constant itself is ${\bf P}_{\perp}$-independent as long as the
replacement $M\to M_{0}$ is made in $P^{-}$.

In this subsection, we shall explicitly prove not only the ${\bf
P}_{\perp}$-independence but also the rotational invariance of the decay
constant $f_{\rm A(P)}$ given by Eq. (27). This can be shown explicitly by
converting Eq. (27) into the integral form of the ordinary three vector
${\vec{k}}=(k_{z},{\bf k}_{\perp})$ using Eqs. (9) and (10) together with the
Jacobi factor given by Eq. (14), which results in

$f_{\rm A(P)}=\sqrt{6}\int\frac{{\rm d}^{3}{\vec{k}}}{16\pi^{3}}\
\sqrt{\frac{M_{0}}{E_{1}E_{2}}}\frac{{\Phi}({\vec{k}})}{{\tilde{M}_{0}}}~{}\mathcal{O}_{\rm
A(P)}({\vec{k}}),$ (30)

where ${\Phi}({\vec{k}})$ now becomes the wave function mixed with
$\phi_{1S}({\vec{k}})$ and $\phi_{2S}({\vec{k}})$ given by Eq. (II). For the
pseudoscalar current case, the rotational invariance of the operator ${\cal
O}_{\rm P}={\tilde{M}}^{2}_{0}/\mu^{0}_{M}$ is evident. For the axial-vector
current case, the operators ${\cal O}^{(+,\perp)}_{\rm A}=2{\cal A}$ can be
converted into

${\cal O}^{(+,\perp)}_{\rm
A}({\vec{k}})=\frac{2}{M_{0}}\left[m_{1}E_{2}+m_{2}E_{1}+(m_{1}-m_{2})k_{z}\right],$
(31)

where the last $k_{z}$-term vanishes for the $k_{z}$-integration in Eq. (30)
and the rest of terms are rotationally invariant. Finally, the operator ${\cal
O}^{-}_{\rm A}$ satisfies

$\displaystyle\mathcal{O}^{+}_{\rm A}({\vec{k}})-\mathcal{O}^{-}_{\rm
A}({\vec{k})}=\frac{4(m_{2}-m_{1})M_{0}}{(\mathbf{P}_{\perp}^{2}+M_{0}^{2})}k_{z},$
(32)

which is an odd function of $k_{z}$. Equation (32) indicates that the decay
constant obtained from the minus current is not only independent of ${\bf
P}_{\perp}$ but also completely equivalent to the one obtained from the plus
(and perpendicular) component of the current. It is worth noting that the
utilization of a factorized form of the LFWFs, such as
$\Psi(x)=\Psi_{1}(x)\Psi_{2}({\bf k}_{\perp})$ Li:2021cwv , would result in
breaking this rotational invariance.

In the numerical section, we will conduct a quantitative analysis to examine
the ${\bf P}_{\perp}$-independence of ${\cal O}^{-}_{\rm A}$ through the ${\bf
P}_{\perp}$-dependence of the helicity contributions.

##  IV Distribution amplitudes

In this section, we discuss the two-particle DAs up to twist 4 obtained from
the three different pseudoscalar meson decay modes. We summarize in Table 2
the twist classification based on the choice of the currents
$(\gamma^{\mu}\gamma_{5},i\gamma_{5},\sigma^{\mu\nu}\gamma_{5})$ and all
possible components of the currents.

Table 2: Twist classification based on the choice of the current and its component. Current | Comp | Twist | DAs
---|---|---|---
$\gamma^{\mu}\gamma_{5}$ | $+,\perp$ | 2 | $\phi_{2;{\rm P}}$
| $-$ | 4 | $\phi_{4;{\rm P}}$
$i\gamma_{5}$ | … | 3 | $\phi_{3;{\rm P}}$
$\sigma^{\mu\nu}\gamma_{5}$ | $+-,\perp-$ | 3 | $\psi_{3;{\rm P}}$

The DAs up to twist-4 accuracy for the pseudoscalar meson with axial-vector
current $\Gamma_{\rm A}$ are defined in terms of the following matrix element
of gauge invariant nonlocal operators as Ball:2006wn

$\displaystyle A_{\rm A}^{\mu}$ $\displaystyle=$
$\displaystyle\bra{0}\bar{q}(z)\gamma^{\mu}\gamma_{5}q(-z)\ket{{\rm P}(P)},$
(33) $\displaystyle=$ $\displaystyle if_{\rm A}\int^{1}_{0}{\rm d}x\ {\rm
e}^{i\zeta P\cdot z}\biggl{[}P^{\mu}\left(\phi_{2;{\rm
P}}(x)+z^{2}(\cdots)\right)$ $\displaystyle\hskip
14.22636pt+\frac{M^{2}}{2}\frac{z^{\mu}}{P\cdot z}\biggl{(}\phi_{4;{\rm
P}}(x)-\phi_{2;{\rm P}}(x)\biggr{)}\biggr{]}.$

In order to make a connection between the DAs and the LFWFs of the meson, we
utilize the equal LF time condition on the light-like vector $z^{\mu}$ (i.e.
$z^{2}=z^{-}z^{+}-\mathbf{z}_{\perp}^{2}=0$) with
$z^{+}=\mathbf{z}_{\perp}=0$. We then obtain

$\displaystyle A^{\mu}_{\rm A}\biggr{|}_{z^{+}=\mathbf{z}_{\perp}=0}$
$\displaystyle=$ $\displaystyle if_{\rm A}\int^{1}_{0}{\rm d}x\ {\rm
e}^{i\zeta P\cdot z}\biggl{[}P^{\mu}\phi_{2;{\rm P}}(x)$ (34)
$\displaystyle+\frac{M^{2}z^{\mu}}{P^{+}z^{-}}\biggl{(}\phi_{4;{\rm
P}}(x)-\phi_{2;{\rm P}}(x)\biggr{)}\biggr{]}.\quad\quad$

To isolate the twist-2 DA, $\phi_{2;{\rm P}}(x)$, one may take either the plus
or transverse component of the current and obtain

$A^{(+,\perp)}_{\rm A}=if_{\rm A}P^{(+,\perp)}\int^{1}_{0}{\rm d}x\ {\rm
e}^{i\zeta P\cdot z}\phi_{2;{\rm P}}(x).\quad\quad$ (35)

This explains why the two decay constants, $f^{(+)}_{\rm A}$ and
$f^{(\perp)}_{\rm A}$, have the same operator ${\cal O}^{(+)}_{\rm P}={\cal
O}^{(\perp)}_{\rm P}$. On the other hand, the twist-4 DA $\phi_{4;{\rm P}}(x)$
can be obtained from the minus component of the current in the ${\bf
P}_{\perp}=0$ frame as

$A^{-}_{\rm A}=if_{\rm A}P^{-}\int^{1}_{0}{\rm d}x\ {\rm e}^{i\zeta P\cdot
z}\phi_{4;{\rm P}}(x).\quad$ (36)

Here it is shown that the higher twist DAs are associated to the bad current
while the leading twist DAs correspond to the good current.

Table 3: The constituent quark masses $m$ [GeV] and variational parameters $\beta$ [GeV] for $\theta=12^{\circ}$ adapted from our previous work Arifi22 .  $m_{q}$ | $m_{s}$ | $m_{c}$ | $m_{b}$ | $\beta_{q\bar{c}}$ | $\beta_{s\bar{c}}$ | $\beta_{q\bar{b}}$ | $\beta_{s\bar{b}}$ | $\beta_{c\bar{c}}$ | $\beta_{c\bar{b}}$ | $\beta_{b\bar{b}}$
---|---|---|---|---|---|---|---|---|---|---
0.22 | 0.45 | 1.68 | 5.10 | 0.424 | 0.455 | 0.495 | 0.538 | 0.592 | 0.767 | 1.167

For the twist-3 case, there are two different DAs that are related to
pseudoscalar ($\Gamma_{\rm P}$) and pseudotensor ($\Gamma_{\rm T}$) currents.
For pseudoscalar current, the twist-3 DA $\phi_{3;{\rm P}}(x)$ is uniquely
determined by Ball:2006wn

$\displaystyle A_{\rm P}\biggr{|}_{z^{+}=\mathbf{z}_{\perp}=0}$
$\displaystyle=$ $\displaystyle\bra{0}\bar{q}(z)i\gamma_{5}q(-z)\ket{{\rm
P}(P)},$ (37) $\displaystyle=$ $\displaystyle f_{\rm P}\mu_{M}\int_{0}^{1}{\rm
d}x\ {\rm e}^{i\zeta P\cdot z}\phi_{3;{\rm P}}(x),$

without choosing a particular component of current. For pseudotensor current,
the DA is computed as Ball:2006wn

$\displaystyle A^{\mu\nu}_{\rm T}\biggr{|}_{z^{+}=\mathbf{z}_{\perp}=0}$
$\displaystyle=$
$\displaystyle\bra{0}\bar{q}(z)\sigma^{\mu\nu}\gamma_{5}q(-z)\ket{{\rm
P}(P)},$ (38) $\displaystyle=$ $\displaystyle-\frac{i}{3}f_{\rm
T}\left(1-\rho_{+}\right)\mu_{M}(P^{\mu}z^{\nu}-P^{\nu}z^{\mu})$
$\displaystyle\times\int_{0}^{1}{\rm d}x\ {\rm e}^{i\zeta P\cdot
z}\psi_{3;{\rm P}}(x).\quad\quad$

In this case, the nonvanishing components are $\mu\nu=+-$ and $\perp-$ and
they give the same $\psi_{3;{\rm P}}(x)$ as we have shown for the computation
of the decay constant $f_{\rm T}.$

In our notation, all DAs $\phi_{n;{\rm P}}(x)(n=2,3,4)$ and $\psi_{3;{\rm
P}}(x)$ are normalized to unity as

$\displaystyle\int_{0}^{1}{\rm d}x\ \left\\{\phi_{n;{\rm P}}(x),\psi_{3;{\rm
P}}(x)\right\\}=1.$ (39)

From Eqs. (27) and (28) together with Eq. (39), we obtain $\phi_{n;{\rm
P}}(x)(n=2,3,4)$ from the axial-vector ($n=2,4$) and pseudoscalar $(n=3)$
channels as

$\displaystyle\phi_{n;\rm P}(x)=\frac{\sqrt{6}}{f_{\rm P}}\int\frac{{\rm
d}^{2}\mathbf{k}_{\bot}}{16\pi^{3}}\
\frac{\Phi(x,\mathbf{k}_{\perp})}{\sqrt{\mathcal{A}^{2}+\mathbf{k}_{\perp}^{2}}}\
\mathcal{O}_{\rm A(P)}.$ (40)

Here, we have $\mathcal{O}^{+}_{\rm A}=\mathcal{O}^{\perp}_{\rm A}$
corresponding to $n=2$, $\mathcal{O}^{-}_{\rm A}$ corresponding to $n=4$, and
$\mathcal{O}_{\rm P}$ corresponding to $n=3$. For $\psi_{3;{\rm P}}(x)$ from
the pseudotensor channel, we obtain

$\displaystyle\psi_{3;\rm P}(x)=\frac{\sqrt{6}}{f_{\rm P}}\int^{x}_{0}{\rm
d}x^{\prime}\int\frac{{\rm d}^{2}\mathbf{k}_{\bot}}{16\pi^{3}}\
\frac{\Phi(x^{\prime},\mathbf{k}_{\perp})}{\sqrt{\mathcal{A}^{\prime
2}+\mathbf{k}_{\perp}^{2}}}\ \mathcal{O}_{\rm T},$ (41)

where $\mathcal{O}_{\rm T}=\mathcal{O}^{+-}_{\rm T}=\mathcal{O}^{\perp-}_{\rm
T}$.

##  V Numerical results and discussion

The model parameters $m$ and $\beta$ used in the present work are summarized
in Table 3, which were determined from the spectroscopic study in our previous
work Arifi22 .

###  V.1 Light-front wave function

We first discuss the LFWFs $\Psi^{00}_{\lambda_{1}\lambda_{2}}(x,{\bf
k}_{\perp})$ defined in Eq. (3) for the ground ($1S$) state and the radially
excited ($2S$) state heavy pseudoscalar mesons. Figure 1(a) shows the two-
dimensional (2D) plots of $D(1S)$ and $D(2S)$ mesons as a function of
$(x,k_{\perp})$, respectively, as an example of the unequal-mass case while
Figure 1(b) shows the 2D plots of $\eta_{c}(1S)$ and $\eta_{c}(2S)$ heavy
quarkonia, respectively, as an example of the equal-mass case. In Fig. 1, the
LFWFs $\Psi^{00}_{\lambda_{1}\lambda_{2}}$ are presented in term of the
helicity configuration $\lambda_{1}\lambda_{2}$ where we denote $\lambda=+1/2$
and $-1/2$ as $\uparrow$ and $\downarrow$, respectively. Note that the
longitudinal momentum fraction $x$ is carried by the lighter quark with mass
$m_{1}$. The LFWFs can be compared with those obtained in Ref. Li:2021cwv .

Figure 1:  Two-dimensional plot of LFWF of (a) $D$ and (b) $\eta_{c}$ mesons
for each helicity configuration. Note that the longitudinal momentum $x$ is
carried by the light quark.

Figure 2:  Helicity contributions to decay constant for 1S and 2S heavy
pseudoscalar mesons. The pattern histogram represents subtracted contribution.
Since the contribution from $H_{\uparrow\downarrow}=H_{\downarrow\uparrow}$,
we sum them up for simplicity. It also applies to
$H_{\uparrow\uparrow}=H_{\downarrow\downarrow}$. Here we set
$\mathbf{P}_{\perp}=0$ for the case of $\gamma^{-}\gamma_{5}$. For the
$\sigma^{\mu\nu}\gamma_{5}$ case, the helicity contribution depends on the
choice of assigning $x$ to light or heavy quark, see Sec. V.3, for the detail.
Figure 3:  The $P_{\perp}$ dependence of the helicity contribution of
$\gamma^{-}\gamma_{5}$ for 1S and 2S pseudoscalar heavy mesons with
$|\mathbf{P}_{\perp}|=P_{\perp}$. The solid line represents the contributions
from $H_{\uparrow\downarrow}+H_{\downarrow\uparrow}$ while the dashed line
represents $H_{\uparrow\uparrow}+H_{\downarrow\downarrow}$.

There are several salient features related to the LFWFs in Fig. 1. (i) The
center of the LFWF ($k_{\perp}\to 0$ and $k_{z}\to 0$), which is associated
with its extremum point, is located at

$\displaystyle x=\frac{m_{1}}{m_{1}+m_{2}},$ (42)

which is obtained by solving Eq. (13). For the equal-mass case,
$\Psi(k_{\perp}\to 0,k_{z}\to 0)$ is located at $x=1/2$ as can be seen in Fig.
1. But, for the unequal-mass case (i.e. $m_{1}=m_{u(d)},m_{2}=m_{c}$) in Fig.
1, the center moves to the value of $x<1/2$ and therefore, the LFWF is
somewhat distorted on $(x,k_{\perp})$ plane. (ii) The LFWF correctly
represents the pseudoscalar meson as
$\Psi_{\uparrow\downarrow}^{00}(x,\mathbf{k}_{\perp})=-\Psi_{\downarrow\uparrow}^{00}(x,\mathbf{k}_{\perp}).$
In addition to the ordinary helicity
$(\uparrow\downarrow,\downarrow\uparrow)$, there is also a nonvanishing
contribution from the higher helicities
$(\uparrow\uparrow,\downarrow\downarrow)$ that couple to the quark orbital
angular momentum as the sign is different in the positive and negative domain
of $k_{\perp}$. This configuration is possible in the relativistic dynamics.
However, such contribution is suppressed as the quark mass increases and
vanishes in the heavy-quark limit ($m\to\infty$). Therefore, at the heavy-
quark or nonrelativistic limit, the LFWF will take only the contribution from
the ordinary helicity without involving the orbital angular momentum. (iii) It
is also shown that the $2S$ state has a nodal structure represented as a white
circle/oval where the center of LFWF has a dip represented as a blue region
for the case of ordinary helicity. One may notice that the shape of the LFWFs
is largely reflected in the DAs.

###  V.2 Decay constant

Figure 4:  The 3D plots of wave functions for $D(1S)$ mesons defined by $f_{A}=\int{\rm d}^{3}\vec{k}\ \psi^{\mu}(\vec{k})$: (a) $\psi^{\mu}(\vec{k})$ with various current component, (b) $\psi^{\mu}(\vec{k})$ without the $k_{z}$ term, and (c) $\tilde{\psi}(\vec{k})=\psi^{+,\perp}(\vec{k})-\psi^{-}(\vec{k})$. Displayed in (d), (e), and (f) is for those of $D(2S)$ mesons. Table 4: Decay constants of heavy pseudoscalar mesons predicted in the LFQM Arifi22 . The results are given in MeV.  State | $f_{theo}$ | $f_{exp}$ | State | $f_{theo}$ | $f_{exp}$
---|---|---|---|---|---
$D(1S)$ | 208 | 206.7(8.9) | $D(2S)$ | 110 | …
$D_{s}(1S)$ | 246 | 257.5(6.1) | $D_{s}(2S)$ | 133 | …
$\eta_{c}(1S)$ | 348 | 335(75) | $\eta_{c}(2S)$ | 214 | …
$B(1S)$ | 190 | 188(25) | $B(2S)$ | 126 | …
$B_{s}(1S)$ | 228 | … | $B_{s}(2S)$ | 150 | …
$B_{c}(1S)$ | 394 | … | $B_{c}(2S)$ | 268 | …
$\eta_{b}(1S)$ | 628 | … | $\eta_{b}(2S)$ | 443 | …

First of all, the numerical values of decay constants for $1S$ and $2S$ state
heavy pseudoscalar mesons obtained from the axial-vector current are presented
in our previous work Arifi22 . In this work, we confirm that the decay
constants obtained from the pseudoscalar and pseudotensor channels also
produce the same results as those obtained from the axial-vector channel
regardless of the currents as well as all possible current components. Namely,
we obtain the process-independent pseudoscalar meson decay constants in the
LFQM. For the sake of completeness, we display again the results of $1S$ and
$2S$ state heavy pseudoscalar mesons for the case of mixing angle
$\theta=12^{\circ}$ in Table 4.

In addition, we examine in Fig. 2 the contributions of helicity to the decay
constants for $1S$ and $2S$ state heavy pseudoscalar mesons, as they are
contingent upon the current component, as indicated in Table 1. Notably, as
observed in Fig. 2, the helicity contributions exhibit variation across
different currents and current components. Despite these variations, however,
the resulting decay constant remains unaltered.

For the $\gamma^{+}\gamma_{5}$, the contribution denoted by
$H_{\uparrow\downarrow}+H_{\downarrow\uparrow}$ is entirely from the ordinary
helicity wave function
$\Psi_{\uparrow\downarrow-\downarrow\uparrow}^{00}(x,\mathbf{k}_{\perp})$
without involving the orbital angular momentum. This is one of the reasons why
we call plus current $(\mu=+)$ as the good current where the dynamics becomes
much simpler and it is also related to the leading twist DAs as explained in
Sec. IV. For the $\gamma^{\perp}\gamma_{5}$, the contribution is still
entirely from the ordinary helicity. However, when we use
$\gamma^{-}\gamma_{5}$ or the bad current, the higher helicity contributions
denoted by $H_{\uparrow\uparrow}+H_{\downarrow\downarrow}$ arise and the
dynamics becomes more complicated. It is clearly shown that the higher
helicity contribution is suppressed when the constituent quark mass becomes
heavier. A similar behavior is also observed for the $i\gamma_{5}$ case. When
we use $\sigma^{+-}\gamma_{5}$ (or $\sigma^{\perp-}\gamma_{5}$), the ordinary
helicity contribution appears more than expected as shown in the bottom left
panel of Fig. 2 for some cases. However, the higher helicity contribution
compensates for it and keeps the decay constant remains the same. It is also
worth mentioning that the behaviors of helicity contribution for the ground
state and the radially excited state are similar. The difference is that the
higher helicity contribution is more pronounced for the radially excited
state.

In Fig. 3, we show the ${\bf P}_{\perp}$-independence of the decay constants
for the $1S$ and $2S$ state heavy pseudoscalar mesons. While each helicity
contribution shows the $\mathbf{P}_{\perp}$ dependence when one uses the minus
component of the axial-vector current, the sum of all helicity contributions
is completely independent of ${\bf P}_{\perp}$ as it should be. It is also
evident from Fig. 3 that the higher (ordinary) helicity contributions dominate
at low (high) $\mathbf{P}_{\perp}$ region consistent with the previous
observation for the equal quark mass case Arifi:2022qnd . We can also see that
higher helicity is more enhanced for the $2S$ state, similar to that in Fig.
2.

Figure 5:  Two-particle DAs with various twists for 1S and 2S heavy
pseudoscalar mesons with various quark flavor contents where the longitudinal
momentum $x$ is carried by the light quark.

Finally, we also examine the rotational invariance of the decay constant by
investigating the wave function $\psi^{\mu}(\vec{k})$ defined in Eq. (30) with
$f_{A}=\int{\rm d}^{3}\vec{k}\ \psi^{\mu}(\vec{k})$ for the axial-vector
$(\Gamma^{\mu}_{\rm A}=\gamma^{\mu}\gamma_{5})$ current. For the equal-mass
case, the $\psi^{\mu}(\vec{k})$ has a spherical shape since the obtained
operator $\mathcal{O}_{A}^{\mu}=2m$ regardless of the current component $\mu$.
For the unequal-mass case such as $D$ meson, the wave functions
$\psi_{D}^{+,\perp}(\vec{k})$ and $\psi_{D}^{-}(\vec{k})$ are slightly
deformed and shifted to the negative and positive $k_{z}$ domains222The
shifting to either positive or negative domain depends on the choice of
$m_{1}$ and $m_{2}$ since
$\mathcal{O}_{A}^{-}(\vec{k})\propto(m_{1}-m_{2})k_{z}$ in Eq. (31). ,
respectively, as depicted in the upper panels of Fig. 4. It is also shown that
the wave functions $\psi^{\mu}(\vec{k})$ are more separated in the $k_{z}$
direction for the $D(1S)$ state compared to those of the $D(2S)$ state. The
shifting in $k_{z}$ direction can be understood by the appearance of the odd
function of $k_{z}$ for ${\cal O}_{A}^{+,\perp}(\vec{k})$ in Eq. (31) although
such an odd $k_{z}$ term does not actually contribute to the decay constant.
The wave functions $\psi^{+,\perp}(\vec{k})$ and $\psi^{-}(\vec{k})$ become a
sphere centered at the origin and coincide with each other if the $k_{z}$ term
is removed as shown in the middle panels of Fig. 4. Moreover, the difference
defined as
$\tilde{\psi}_{D}(\vec{k})=\psi^{+}_{D}(\vec{k})-\psi^{-}_{D}(\vec{k})$ is
displayed in lower panels of Fig. 4 showing that the integration over $k_{z}$
will give a vanishing result. Therefore, it is evident that the decay constant
with $\mu=+,\perp,-$ is the same. As for the pseudoscalar current, the wave
functions are spherical as also implied from its operator ${\cal
O}_{P}(\vec{k})$.

###  V.3 Distribution amplitude

Figure 6:  Two-dimensional plot of the the DAs of $D(2S)$ for various twist.
Here we define ${\bf k}_{\perp}^{2}=y/(1-y)$ to make a rectangular domain. The
longitudinal momentum $x$ is carried by the light quark. Figure 7:  Helicity
contribution to DAs with various twist. Although the total DAs of twist 2,
twist 3, and twist 4 are the same, they consist of different helicity
contributions. We note that the total DA of twist 3 with
$\sigma^{+-}\gamma_{5}$ is different from the others.  Figure 8:  The
difference in the helicity contribution depending on the choice of assigning
the LF longitudinal momentum fraction $x$ to the light or heavy quark. The
dashed and dotted lines represent the ordinary and higher helicity,
respectively.  Figure 9:  Comparison of $\psi_{3;P}(x)$ with and without the
inclusion of $(1-\rho_{+})$ factor.

Figure 5 presents the DAs of different twists for the $1S$ and $2S$ state
heavy pseudoscalar mesons. Note here that the longitudinal momentum $x$ is
carried by the lighter quark. As a result, the DAs for $D_{(s)}$ and $B_{(s)}$
are more concentrated in the lower $x$ region. For the equal mass case such as
$\eta_{c}$ and $\eta_{b}$, the $\phi_{2;{\rm P}}(x),$ $\phi_{3;{\rm P}}(x)$
and $\phi_{4;{\rm P}}(x)$ have the same shape since the corresponding
operators are the same as shown in Table 7. The key reason for this is the
utilization of the self-consistent condition for the replacement of $M$ with
$M_{0}$ when obtaining the aforementioned results. But, the $\psi_{3;{\rm
P}}(x)$ has a different shape with a narrower and higher peak. In addition,
the distance between the peaks becomes closer for $\psi_{3;{\rm P}}(x)$ of the
$2S$ state. For the unequal-mass case such as $D$ or $B$ meson, the peak is
shifted to the lower $x$ region for the higher twist DAs where the
$\phi_{4;{\rm P}}(x)$ has the highest peak. For the $2S$ state, the peak near
$x=0$ is more enhanced for the higher twist while the peak near $x=0.5$ is
suppressed and shifted to the lower $x$ region. Furthermore, the dip between
the peaks is enhanced for the higher twist.

In order to gain a deeper understanding of the structure of DAs, we construct
a 2D plot illustrating the DAs, $\phi(x)\equiv\\{\phi_{n;{\rm P}},\psi_{3;{\rm
P}}\\}$, using the following definition:

$\displaystyle\phi(x)=\int^{\infty}_{0}{\rm d}^{2}{\bf k}_{\perp}\psi(x,{\bf
k}_{\perp})=\int_{0}^{1}{\rm d}y\ \phi(x,y),$ (43)

where the wave function $\phi(x,y)=\pi\psi(x,y)/(1-y)^{2}$ is obtained by
using the variable transformation ${\bf k}_{\perp}^{2}=y/(1-y)$ so that $y$
ranges from 0 to 1. For the sake of demonstration, we show the 2D plot of
$\phi(x,y)$ only for $D(2S)$ as shown in Fig. 6. One can clearly see that the
$\phi(x,y)$ bears resemblance with the LFWF shown in Fig. 1(a). The two-peak
structure in the DAs of the $2S$ state clearly originated from the nodal
structure shown as the white bands. It appears that the $\phi_{2;{\rm
P}}(x,y)$, $\phi_{3;{\rm P}}(x,y)$, and $\phi_{4;{\rm P}}(x,y)$ have a similar
shape. But, the DAs with the higher twist are more concentrated in the lower
$x$ region. On the other hand, $\psi_{3;{\rm P}}(x,y)$ has a smaller nodal
structure so that it shows slightly different behavior in Fig. 5.

Figure 7 shows the helicity contributions to DAs of $\eta_{c}(1S)$ for various
twist. The dashed and dotted lines represent the ordinary
$H_{\uparrow\downarrow}+H_{\downarrow\uparrow}$ and higher
$H_{\uparrow\uparrow}+H_{\downarrow\downarrow}$ helicity contributions,
respectively. The solid lines represent the full results. As mentioned
earlier, the DAs of twist 2, twist 3, and twist 4 for equal-mass cases are the
same. However, these DAs exhibit distinct helicity contributions.
Specifically, the twist-2 DAs are exclusively composed of the ordinary
helicity component, while the higher twist DAs incorporate a finite
contribution from the higher helicity component. This observation is in
accordance with the result presented for helicity contribution to the decay
constant in Fig. 2.

Although the helicity contribution to DAs generally remains unchanged
regardless of whether the light quark is assigned to $x$ or $(1-x)$, it is
crucial to acknowledge that the specific choice between $x$ and $1-x$ does
impact the helicity contribution to the DA $\psi_{3;{\rm P}}(x)$ obtained from
the nonlocal matrix element. In particular, when assigning the light quark to
either $x$ or $(1-x)$, the helicity contribution to $\psi_{3;{\rm P}}(x)$
exhibits markedly distinct behaviors. Figure 8 illustrates, as an example, the
discrepancy in helicity contributions to $\psi_{3;{\rm P}}(x)$ for $D(1S)$
when the light quark is assigned to carry either $x$ or $(1-x)$. The upper and
lower panels in Fig. 8 represent the results obtained when the light quark is
assigned to $x$ and $1-x$, respectively, and the same line codes are used as
in Fig. 7. One can clearly see from Fig. 8 that the ordinary helicity exhibits
a negative contribution when the light quark carries the value of $x$ (upper
panel), while it yields a positive contribution when the heavy quark carries
the value of $x$ (lower panel). This distinct behavior observed can be
attributed to the integration over $x^{\prime}$, where the behavior depends on
the specific choice of $x$. It indicates that each individual ordinary and
higher helicity contribution to $\psi_{3;{\rm P}}(x)$ depends on whether the
light quark or the heavy quark carries the specific light-front longitudinal
momentum fraction $x$. However, it is important to note that the total DA
remains unchanged whether we assign $x$ to the light or heavy quark.
Furthermore, it is worth noting that there is a substantial cancellation
between the ordinary and higher helicity contributions when $x$ is associated
with the light quark, whereas the cancellation is comparatively smaller when
$x$ is assigned to the heavy quark.

Table 5:  The $\xi$-moment up to $n=6$ for the $1S$ state heavy pseudoscalar mesons. Here we define the $x$ carried by the lighter quark. Therefore, the odd-power of $\expectationvalue{\xi}$ has an opposite sign to our previous work Arifi22 .  $(1S)$ | tw | $D$ | $D_{s}$ | $\eta_{c}$ | $B$ | $B_{s}$ | $B_{c}$ | $\eta_{b}$
---|---|---|---|---|---|---|---|---
$\expectationvalue{\xi^{1}}$ | 2 | $-0.337$ | $-0.294$ | … | $-0.644$ | $-0.614$ | $-0.390$ | …
| 3p | $-0.445$ | $-0.365$ | … | $-0.713$ | $-0.670$ | $-0.419$ | …
| 4 | $-0.553$ | $-0.436$ | … | $-0.781$ | $-0.726$ | $-0.447$ | …
| 3t | $-0.445$ | $-0.365$ | … | $-0.713$ | $-0.670$ | $-0.417$ | …
$\expectationvalue{\xi^{2}}$ | 2 | 0.226 | 0.197 | 0.088 | 0.453 | 0.417 | 0.201 | 0.049
| 3p | 0.312 | 0.242 | 0.088 | 0.545 | 0.487 | 0.223 | 0.049
| 4 | 0.397 | 0.288 | 0.088 | 0.636 | 0.558 | 0.246 | 0.049
| 3t | 0.273 | 0.202 | 0.053 | 0.533 | 0.475 | 0.205 | 0.029
$\expectationvalue{\xi^{3}}$ | 2 | $-0.145$ | $-0.114$ | … | $-0.337$ | $-0.302$ | $-0.113$ | …
| 3p | $-0.222$ | $-0.154$ | … | $-0.435$ | $-0.373$ | $-0.129$ | …
| 4 | $-0.299$ | $-0.193$ | … | $-0.534$ | $-0.445$ | $-0.146$ | …
| 3t | $-0.177$ | $-0.116$ | … | $-0.413$ | $-0.350$ | $-0.108$ | …
$\expectationvalue{\xi^{4}}$ | 2 | 0.108 | 0.083 | 0.018 | 0.261 | 0.228 | 0.068 | 0.006
| 3p | 0.173 | 0.112 | 0.018 | 0.360 | 0.296 | 0.080 | 0.006
| 4 | 0.238 | 0.142 | 0.018 | 0.458 | 0.364 | 0.092 | 0.006
| 3t | 0.124 | 0.072 | 0.008 | 0.329 | 0.265 | 0.061 | 0.003
$\expectationvalue{\xi^{5}}$ | 2 | $-0.082$ | $-0.058$ | … | $-0.209$ | $-0.178$ | $-0.043$ | …
| 3p | $-0.138$ | $-0.082$ | … | $-0.304$ | $-0.241$ | $-0.052$ | …
| 4 | $-0.195$ | $-0.106$ | … | $-0.399$ | $-0.304$ | $-0.060$ | …
| 3t | $-0.090$ | $-0.048$ | … | $-0.267$ | $-0.206$ | $-0.035$ | …
$\expectationvalue{\xi^{6}}$ | 2 | 0.065 | 0.044 | 0.005 | 0.172 | 0.143 | 0.029 | 0.001
| 3p | 0.115 | 0.063 | 0.005 | 0.262 | 0.200 | 0.035 | 0.001
| 4 | 0.165 | 0.083 | 0.005 | 0.352 | 0.258 | 0.041 | 0.001
| 3t | 0.068 | 0.033 | 0.002 | 0.220 | 0.162 | 0.021 | 0.0004
Table 6:  The $\xi$-moment up to $n=6$ for the $2S$ state heavy pseudoscalar mesons. Here we define the $x$ carried by the lighter quark. Therefore, the odd-power of $\expectationvalue{\xi}$ has an opposite sign to our previous work Arifi22 .  $(2S)$ | tw | $D$ | $D_{s}$ | $\eta_{c}$ | $B$ | $B_{s}$ | $B_{c}$ | $\eta_{b}$
---|---|---|---|---|---|---|---|---
$\expectationvalue{\xi^{1}}$ | 2 | 0.042 | $-0.015$ | … | $-0.426$ | $-0.411$ | $-0.275$ | …
| 3p | $-0.161$ | $-0.155$ | … | $-0.549$ | $-0.518$ | $-0.333$ | …
| 4 | $-0.363$ | $-0.296$ | … | $-0.672$ | $-0.625$ | $-0.390$ | …
| 3t | $-0.161$ | $-0.155$ | … | $-0.549$ | $-0.518$ | $-0.333$ | …
$\expectationvalue{\xi^{2}}$ | 2 | 0.052 | 0.132 | 0.179 | 0.198 | 0.202 | 0.160 | 0.099
| 3p | 0.164 | 0.199 | 0.179 | 0.344 | 0.322 | 0.201 | 0.099
| 4 | 0.275 | 0.266 | 0.179 | 0.491 | 0.442 | 0.241 | 0.099
| 3t | 0.081 | 0.115 | 0.107 | 0.320 | 0.295 | 0.163 | 0.059
$\expectationvalue{\xi^{3}}$ | 2 | 0.004 | $-0.055$ | … | $-0.094$ | $-0.112$ | $-0.101$ | …
| 3p | $-0.099$ | $-0.121$ | … | $-0.237$ | $-0.225$ | $-0.131$ | …
| 4 | $-0.202$ | $-0.187$ | … | $-0.380$ | $-0.337$ | $-0.161$ | …
| 3t | $-0.029$ | $-0.061$ | … | $-0.194$ | $-0.180$ | $-0.092$ | …
$\expectationvalue{\xi^{4}}$ | 2 | 0.012 | 0.065 | 0.048 | 0.043 | 0.070 | 0.071 | 0.016
| 3p | 0.091 | 0.113 | 0.048 | 0.175 | 0.171 | 0.093 | 0.016
| 4 | 0.171 | 0.161 | 0.048 | 0.307 | 0.272 | 0.115 | 0.016
| 3t | 0.018 | 0.048 | 0.021 | 0.120 | 0.116 | 0.057 | 0.007
$\expectationvalue{\xi^{5}}$ | 2 | $-0.006$ | $-0.047$ | … | $-0.017$ | $-0.048$ | $-0.051$ | …
| 3p | $-0.076$ | $-0.089$ | … | $-0.137$ | $-0.138$ | $-0.067$ | …
| 4 | $-0.146$ | $-0.131$ | … | $-0.257$ | $-0.228$ | $-0.084$ | …
| 3t | $-0.011$ | $-0.035$ | … | $-0.075$ | $-0.078$ | $-0.037$ | …
$\expectationvalue{\xi^{6}}$ | 2 | 0.010 | 0.044 | 0.016 | 0.004 | 0.036 | 0.038 | 0.003
| 3p | 0.070 | 0.079 | 0.016 | 0.112 | 0.116 | 0.050 | 0.003
| 4 | 0.129 | 0.113 | 0.016 | 0.220 | 0.196 | 0.063 | 0.003
| 3t | 0.009 | 0.028 | 0.005 | 0.047 | 0.056 | 0.025 | 0.001

As we previously explained in the definition of $\psi_{3;{\rm P}}(x)$ provided
by Eq.(III), there exist two variations in the QCD sum-rules for defining
$\psi_{3;{\rm P}}(x)$, namely, with the inclusion of $(1-\rho_{+})$Ball90  or
without it Ball:2006wn . In the previous work Choi:2017uos , the $\psi_{3;{\rm
P}}(x)$ DA was computed without the term $(1-\rho_{+})$. However, we found
that the inclusion of $(1-\rho_{+})$ is pivotal in obtaining the identical
decay constants from the nonlocal pseudotensor channel as those derived from
the local axial-vector and pseudoscalar channels.

Figure 9 depicts a comparison of $\psi_{3;P}(x)$ obtained with the term
$(1-\rho_{+})$ (solid lines) and without it (dashed lines), for the cases of
heavy $\eta_{c}(1S)$ (upper panel) and the light $\pi(1S)$ (lower panel)
mesons. We note that the same model parameters are used as in Choi:2017uos
for the plots of $\pi$. The analysis reveals that the inclusion of the term
$(1-\rho_{+})$ in $\psi_{3;{\rm P}}(x)$ leads to the narrower and higher shape
compared to the case where $(1-\rho_{+})$ is absent. The quantitative impact
of $(1-\rho_{+})$ on $\psi_{3;P}(x)$ is found to be more significant in the
heavy quark sector compared to the light quark sector.

Finally, we also compute the $\xi$-moment up to $n=6$ defined by

$\displaystyle\expectationvalue{\xi^{n}}=\int_{0}^{1}{\rm d}x\ \xi^{n}\
\phi(x),$ (44)

where $\xi=x-(1-x)=2x-1$. The results are shown in Tables 5 and 6 for the $1S$
and $2S$ state heavy pseudoscalar mesons, respectively. Here, we note that $x$
is carried by the light quark, so the sign is opposite for the odd
$\xi$-moment as compared to our previous work Arifi22 . For the unequal-mass
case, we observe that the absolute value of odd $\xi$-moment is getting larger
for the higher twist, indicating the DAs have more deviated from the center.
The minus sign shows that the DAs are shifted to the lower $x$ region. We also
note that the absolute value of the even $\xi$ moments increases for the
higher twist. The absolute values of $\xi$ moments for the $\psi_{3;{\rm
P}}(x)$ are generally smaller than those of $\phi_{3;{\rm P}}(x)$.

##  VI Summary

We have investigated the decay constants and the DAs up to the twist-4 for the
$1S$ and $2S$ state heavy pseudoscalar mesons in the LFQM, computing the local
and nonlocal matrix elements $\langle 0|{\bar{q}}{\Gamma}q|P\rangle$ with
three different current operators
${\Gamma}=(\gamma^{\mu}\gamma_{5},i\gamma_{5},\sigma^{\mu\nu}\gamma_{5})$.

In our LFQM, we performed a comprehensive analysis utilizing a general
reference frame where ${\bf P}_{\perp}\neq 0$ and explored all possible
components of the currents. Our explicit results demonstrate the equality of
the three pseudoscalar meson decay constants derived from the three distinct
current operators $\Gamma$. This remarkable consistency in decay constants is
achieved when we enforce the self-consistency condition, i.e. the replacement
of the physical mass $M$ with the invariant mass $M_{0}$, within the LFQM.
This condition stems from the Bakamjian-Thomas (BT) construction, in which the
meson state is based on a noninteracting quark-antiquark representation. It is
important to note that the inclusion of the $(1-\rho_{+})$ factor in the
definition of the nonlocal matrix elements $\langle
0|{\bar{q}}(z)\sigma^{\mu\nu}\gamma_{5}q(-z)|P\rangle$ is crucial in order to
obtain the same decay constant as those derived from the axial-vector and
pseudoscalar currents. In addition to secure the process-independent
pseudoscalar meson decay constant, regardless of the choice of current
operators $\Gamma$, we also demonstrated its explicit Lorentz and rotation
invariance.

We also examined the helicity contributions to the decay constants, offering
additional insights into the structural aspects of the decay constant. While
the decay constant remains unchanged regardless of the choice of currents, the
helicity contributions to the decay constant exhibit variations depending on
the specific current and its components, as illustrated in Table 1. As
illustrated in Fig. 2, while the good (plus) current only receives the
ordinary helicity contributions ($\uparrow\downarrow,\downarrow\uparrow$), the
other components including the bad (minus) current receive the higher helicity
contributions $(\uparrow\uparrow,\downarrow\downarrow)$. We further explored
the impact of $\mathbf{P}_{\perp}$ dependence on the helicity contributions
when considering the axial-vector current with the minus current component.
Notably, it becomes evident that the higher helicity contributions play a more
prominent role in the low $\mathbf{P}_{\perp}$ region, particularly for the
$2S$ state. These observations are depicted in Fig. 3.

According to the classification provided in Table 2, employing various current
operators and different components of the currents leads to distinct twists in
the DAs. In particular, we explored the twist-4 DA derived from the minus
component of the axial-vector current.

The various twist DAs for the $1S$ and $2S$ heavy pseudoscalar mesons are
exhibited in Fig. 5. It is evident that the higher twist DAs for the unequal-
mass case is more concentrated in the lower $x$ region. The $\xi$ moments are
also computed for the various twist DAs. One of the notable results is that
the odd $\xi$-moment is getting larger for the higher twist.

We expect that our result is useful for the calculation of hard exclusive
processes based on the QCD factorization. Especially, the higher twist DAs may
be important in the low-$Q^{2}$ region Ryu:2018egt ; Raha:2010kz . It would be
of great importance to extend our analysis to the vector mesons with the
longitudinal and transverse polarizations for further test of our methodology
Arifi:2022qnd . Moreover, the investigation of the decay constant for the
excited scalar, axial-vector, and tensor mesons would be also interesting to
confirm whether our LFQM based on the BT construction can universally be
applicable regardless of the meson quantum numbers Chang:2018zjq ; Cheng04 .
Finally, the extension of our approach to encompass three-point functions,
such as elastic or transition form factors, deserves a thorough investigation
to explore the LF zero-mode effects.

## Acknowledgement

The work of A.J.A. is supported by the RIKEN special postdoctoral researcher
program and the Young Scientist Training (YST) Program at the Asia Pacific
Center for Theoretical Physics (APCTP) through the Science and Technology
Promotion Fund and Lottery Fund of the Korean Government and also by the
Korean Local Governments – Gyeongsangbuk-do Province and Pohang City. The work
of H.-M.C. was supported by the National Research Foundation of Korea (NRF)
under Grant No. NRF- 2023R1A2C1004098. The work of C.-R.J. was supported in
part by the U.S. Department of Energy (Grant No. DE-FG02-03ER41260). The
National Energy Research Scientific Computing Center (NERSC) supported by the
Office of Science of the U.S. Department of Energy under Contract No. DE-
AC02-05CH11231 is also acknowledged.

Table 7:  BS and LFQM operators of decay constant with various currents and their corresponding DA for both unequal and equal masses of the constitutes. Here we have defined $\mathcal{A}=(1-x)m_{1}+xm_{2}$, $\Delta_{1}=(m_{1}^{2}+\mathbf{k}_{\perp}^{2})/x$, $\Delta_{2}=(m_{2}^{2}+\mathbf{k}_{\perp}^{2})/(1-x)$, $\mu^{0}_{M}=M_{0}^{2}/(m_{1}+m_{2})$, and $\rho^{0}_{+}=(m_{1}+m_{2})^{2}/M_{0}^{2}$. It is evident that the BS operator will be the same with the LFQM operator if we apply the $M\to M_{0}$ replacement.  $\mathcal{O}$ | Current | $\mathcal{O}_{\rm BS}^{\rm on}(m_{1}\neq m_{2})$ |  $\mathcal{O}_{\rm LFQM}$($m_{1}\neq m_{2})$ | $\mathcal{O}_{\rm BS}^{\rm on}(m)$ | $\mathcal{O}_{\rm LFQM}(m)$
---|---|---|---|---|---
$\mathcal{O}_{\rm A}^{\mu}$ | $\Gamma_{\rm A}^{+},\Gamma_{\rm A}^{\perp}$ | $2\mathcal{A}$ | $2\mathcal{A}$ | $2m$ | $2m$
$\Gamma_{\rm A}^{-}$ | $\dfrac{2(m_{1}\Delta_{2}+m_{2}\Delta_{1}+\mathcal{A}\mathbf{P}_{\perp}^{2})}{(M^{2}+\mathbf{P}_{\perp}^{2})}$ | $\dfrac{2(m_{1}\Delta_{2}+m_{2}\Delta_{1}+\mathcal{A}\mathbf{P}_{\perp}^{2})}{(M_{0}^{2}+\mathbf{P}_{\perp}^{2})}$ | $2m\dfrac{(M_{0}^{2}+\mathbf{P}_{\perp}^{2})}{(M^{2}+\mathbf{P}_{\perp}^{2})}$ | $2m$
$\mathcal{O}_{\rm P}$ | $\Gamma_{\rm P}$ | $\dfrac{{\tilde{M}_{0}}^{2}}{\mu_{M}}$ | $\dfrac{{\tilde{M}_{0}}^{2}}{\mu^{0}_{M}}$ | $2m\dfrac{M_{0}^{2}}{M^{2}}$ | $2m$
$\mathcal{O}_{\rm T}^{\mu\nu}$ | $\Gamma_{\rm T}^{+-}$, $\Gamma_{\rm T}^{\perp-}$ | $\dfrac{-12M^{\prime}_{0}k^{\prime}_{z}}{\mu_{M}(1-\rho_{+})}$ | $\dfrac{-12M^{\prime}_{0}k^{\prime}_{z}}{\mu^{\prime 0}_{M}(1-\rho^{\prime 0}_{+})}$ | $\dfrac{12m(1-2x^{\prime}){M^{\prime}}_{0}^{2}}{M^{2}-4m^{2}}$ | $\dfrac{12m(1-2x^{\prime}){M}^{\prime 2}_{0}}{{M}^{\prime 2}_{0}-4m^{2}}$

##  Appendix A Link between the Covariant BS Model and the LFQM

As we explained in the Introduction, our self-consistent LFQM results, e.g.
Eqs. (27) and (28) in this work, can also be obtained from the “Type II” link
between the manifestly covariant BS model and the LFQM, which was first
introduced in C13 . This is another approach to arrive at the self-consistent
LFQM. Since the detailed analysis for the link between the manifestly
covariant BS model and the LFQM has already been made in the previous works
C15 ; Choi:2017uos , we shall briefly discuss the essential feature of the
“Type II” link starting from the covariant BS model in this Appendix.

The matrix element $A_{\rm A(P)}\equiv\bra{0}\bar{q}\Gamma_{\rm A(P)}q\ket{P}$
for the local operators $\Gamma_{\rm A(P)}$ in the covariant BS model is given
in the one-loop approximation as

$A_{\rm A(P)}=iN_{c}\int\frac{{\rm d}^{4}k}{(2\pi)^{4}}\frac{H_{0}S_{\rm
A(P)}}{(p^{2}_{1}-m^{2}_{1}+i\epsilon)(p^{2}_{2}-m^{2}_{2}+i\epsilon)},$ (45)

where $S_{\rm A(P)}={\rm Tr}[\Gamma_{\rm
A(P)}(\not{p}_{1}+m_{1})\gamma_{5}(-\not{p}_{2}+m_{2})]$ is the trace term
with $p_{1}=P-k$ and $p_{2}=k$. To regularize the loop, we use the usual
multipole ansatz $H_{0}=\frac{g}{D_{\Lambda}^{2}}$ with
$D_{\Lambda}=p_{1}^{2}-m_{\Lambda}^{2}+i\epsilon$, where $m_{\Lambda}$ plays
the role of the momentum cut-off.

In this exactly solvable manifestly covariant BS model, the decay constants
for the axial-vector and pseudoscalar currents can be obtained from the
manifestly covariant calculation using the Feynman parametrization and the
final results are given by

$\displaystyle f_{\rm A}$ $\displaystyle=$
$\displaystyle\frac{gN_{c}}{4\pi^{2}}\int_{0}^{1}{\rm d}x\int_{0}^{1-x}{\rm
d}y\ \frac{(1-x-y)B_{1}}{C^{2}},$ (46) $\displaystyle f_{\rm P}$
$\displaystyle=$ $\displaystyle\frac{gN_{c}}{4\pi^{2}\mu_{M}}\int_{0}^{1}{\rm
d}x\int_{0}^{1-x}{\rm d}y\ \frac{(1-x-y)(B_{2}-2C)}{C^{2}},$

where $C=y(1-y)M^{2}-xm_{1}^{2}-ym_{2}^{2}-(1-x-y)m_{\Lambda}^{2}$,
$B_{1}=m_{2}+(1-y)(m_{1}-m_{2})$, and $B_{2}=y(1-y)M^{2}+m_{1}m_{2}$. We
should note, at this point, that the two pseudoscalar meson decay constants
$f_{\rm A}$ and $f_{\rm P}$ obtained in the BS model are not the same each
other, e.g. $f_{\rm A}=208$ MeV vs. $f_{\rm P}=225$ MeV for $D(1S)$ meson with
the value of $m_{\Lambda}=1.673$ GeV. This contrasts to our LFQM in which we
obtain the process-independent decay constant.

In parallel with the manifestly covariant calculation, we perform the LF
calculation of Eq. (45) by doing the LF energy integration $p^{-}_{2}$ picking
up the on-mass shell pole $p^{2}_{2}=p^{2}_{2\rm on}=m^{2}_{2}$ and obtain

$f_{\rm A(P)}=N_{c}\int^{1}_{0}{\rm d}x\int\frac{{\rm
d}^{2}\mathbf{k}_{\bot}}{8\pi^{3}}\
\frac{\chi(x,\mathbf{k}_{\perp})}{1-x}{\cal O}^{\rm A(P)}_{\rm BS},$ (48)

where
$\chi(x,\mathbf{k}_{\perp})=1/([x(M^{2}-M_{0}^{2})][x(M^{2}-M_{\Lambda}^{2})]^{2})$
is the vertex function with $M^{2}_{\Lambda}=M^{2}_{0}(m_{1}\to m_{\Lambda})$
and ${\cal O}^{\rm A(P)}_{\rm BS}=iS_{\rm A(P)}/2\cal P_{\rm A(P)}$.

In contrast to the LFQM constrained by the on-mass shellness of the
constituents, the LF calculation of the BS model allows the off-mass shell
quark propagators. For the axial-vector current $\Gamma^{\mu}_{\rm
A}=\gamma^{\mu}\gamma_{5}$ with the current component $\mu=(+,\perp)$, we find
that only the on-mass shell quark propagators contribute and the full result
of the operator is obtained as $[{\cal O}^{\rm A}_{\rm BS}]_{\rm full}=[{\cal
O}^{\rm A}_{\rm BS}]^{+}_{\rm on}=[{\cal O}^{\rm A}_{\rm BS}]^{\perp}_{\rm
on}=2{\cal A}$. On the other hand, the minus component of the axial-vector
current receives not only the instantaneous but also the zero-mode
contributions in addition to the on-mass shell contribution, i.e. $[{\cal
O}^{\rm A}_{\rm BS}]_{\rm full}=[{\cal O}^{\rm A}_{\rm BS}]^{-}_{\rm
on}+[{\cal O}^{\rm A}_{\rm BS}]^{-}_{\rm inst}+[{\cal O}^{\rm A}_{\rm
BS}]^{-}_{\rm z.m.}=2{\cal A}$, where

$\displaystyle\left[{\cal O}^{\rm A}_{\rm BS}\right]^{-}_{\rm on}$
$\displaystyle=$ $\displaystyle\frac{2(m_{1}\Delta_{1}+m_{2}\Delta_{2}+{\cal
A}{\bf P}^{2}_{\perp})}{M^{2}+{\bf P}^{2}_{\perp}},$ $\displaystyle\left[{\cal
O}^{\rm A}_{\rm BS}\right]^{-}_{\rm inst}$ $\displaystyle=$
$\displaystyle\frac{2m_{2}(M^{2}-M^{2}_{0})}{M^{2}+{\bf P}^{2}_{\perp}},$
$\displaystyle\left[{\cal O}^{\rm A}_{\rm BS}\right]^{-}_{\rm z.m.}$
$\displaystyle=$ $\displaystyle\frac{2(m_{1}-m_{2})Z_{2}}{M^{2}+{\bf
P}^{2}_{\perp}},$ (49)

with $\Delta_{j}=(m^{2}_{j}+{\bf k}^{2}_{\perp})/x_{j}(j=1,2)$ and
$Z_{2}=x(M^{2}-M^{2}_{0})+m^{2}_{1}-m^{2}_{2}+(1-2x)M^{2}$.

For the pseudoscalar current $\Gamma_{\rm P}=i\gamma_{5}$, the full operator
is obtained from the sum of the three nonvanishing contributions, i.e. $[{\cal
O}^{\rm P}_{\rm BS}]_{\rm
full}=2[m_{1}(m_{2}-m_{1})+xM^{2}_{0}]/\mu_{M}=[{\cal O}^{\rm P}_{\rm
BS}]_{\rm on}+[{\cal O}^{\rm P}_{\rm BS}]_{\rm inst}+[{\cal O}^{\rm P}_{\rm
BS}]_{\rm z.m.}$, where

$\displaystyle\left[{\cal O}^{\rm P}_{\rm BS}\right]_{\rm on}$
$\displaystyle=$ $\displaystyle\frac{{\tilde{M}}^{2}_{0}}{\mu_{M}},$
$\displaystyle\left[{\cal O}^{\rm P}_{\rm BS}\right]_{\rm inst}$
$\displaystyle=$ $\displaystyle\frac{(1-x)(M^{2}-M^{2}_{0})}{\mu_{M}},$
$\displaystyle\left[{\cal O}^{\rm P}_{\rm BS}\right]_{\rm z.m.}$
$\displaystyle=$ $\displaystyle-\frac{Z_{2}}{\mu_{M}},$ (50)

With those full operators $[{\cal O}^{\rm A(P)}_{\rm BS}]_{\rm full}$, the LF
results for $f_{\rm A(P)}$ given by Eq. (48) are the same as the corresponding
covariant ones given by Eq. (46).

The basic idea of Ref. C13  in obtaining the LFQM result given by Eq. (27)
from the BS model amplitude given by Eq. (48) is to replace not only vertex
function $\chi(x,{\bf k}_{\perp})$ in Eq. (48) with the Gaussian wave function
$\Phi(x,{\bf k}_{\perp})$ but also all the physical mass $M$ appeared in the
BS model with the invariant mass $M_{0}$ via the “Type II” link between the BS
model and the LFQM as we coined in C13 :

$\displaystyle\sqrt{2N_{c}}\frac{\chi(x,\mathbf{k}_{\perp})}{(1-x)}$
$\displaystyle\to$
$\displaystyle\frac{\Phi(x,\mathbf{k}_{\perp})}{\sqrt{\mathcal{A}^{2}+\mathbf{k}_{\perp}^{2}}},$
$\displaystyle M$ $\displaystyle\to$ $\displaystyle M_{0}.$ (51)

The first immediate action for the replacement of $M\to M_{0}$ in the LFQM is
to remove the instantaneous contribution ($\propto M^{2}-M^{2}_{0}$), which
may appear in the covariant BS model but absent in the LFQM consistent with
the BT construction. The only spurious effect, which may appear in the LFQM,
is the LF zero mode. Furthermore, a crucial feature when utilizing the “Type
II” link between the BS model and the LFQM is solely to employ the on-mass
shell BS operator. In other words, the LFQM operator, denoted as ${\cal
O}\equiv{\cal O}_{\rm LFQM}$ and defined by Eq. (27), can be directly derived
by substituting $\left[{\cal O}_{\rm BS}\right]_{\rm on}(M\to M_{0})$.

For the axial-vector current case, the full BS operator $\left[{\cal O}^{\rm
A}_{\rm BS}\right]_{\rm full}=2{\cal A}$ obtained from $\left[{\cal O}^{\rm
A}_{\rm BS}\right]^{+}_{\rm on}=\left[{\cal O}^{\rm A}_{\rm
BS}\right]^{\perp}_{\rm on}=\left[{\cal O}^{\rm A}_{\rm BS}\right]^{-}_{\rm
full}$ is shown exactly the same as the ${\cal O}^{\rm A}_{\rm LFQM}$. The
plus and perpendicular components of the current are free from the
instantaneous and zero-mode contributions, and they are indeed the “good”
components of the current. While the full operator derived from the minus
component of the current shares the exact same form as those derived from the
“good” currents, this feature can be considered extremely rare. It should be
noted that the full operator obtained by including the zero mode does not
generally align with the full operator derived solely from the on-mass shell
contribution, as evident in the case of the pseudoscalar current in which the
full operator $[{\cal O}^{\rm P}_{\rm BS}]_{\rm full}$ is different from the
on-mass shell operator $[{\cal O}^{\rm P}_{\rm BS}]_{\rm on}$ and one finds
that $f_{\rm P}$ in Eq. (27) does not match with $f_{\rm A}$ if the full
operator $[{\cal O}^{\rm P}_{\rm BS}]_{\rm full}$ is used for the replacement
of $M\to M_{0}$ instead of $[{\cal O}^{\rm P}_{\rm BS}]_{\rm on}$. A similar
observation has also been made in the previous analysis for the vector meson
decay constant C13 . This indicates that the zero mode found in the BS model
is no longer applicable to the LFQM. Instead, the replacement of $M\to M_{0}$
for the on-mass shell operator $\left[{\cal O}_{\rm BS}\right]_{\rm on}$, e.g.
$\left[{\cal O}^{\rm A}_{\rm BS}\right]^{-}_{\rm on}(M\to M_{0})$ and
$\left[{\cal O}^{\rm P}_{\rm BS}\right]_{\rm on}(M\to M_{0})$, can be regarded
as an effective zero-mode inclusion in the LFQM.

For the pseudotensor current, the local matrix element
$\bra{0}\bar{q}(0)\Gamma_{\rm T}q(0)\ket{P}$ defined by Eq. (45) in the
covariant BS model is zero in the manifestly covariant calculation. Performing
the LF calculation for this local matrix element, we also confirm that the
full operator $[{\cal O}^{\rm T}_{\rm BS}]_{\rm full}$ is zero if and only if
we include both nonvanishing instantaneous and zero-mode contributions, i.e.,
$[{\cal O}^{\rm T}_{\rm BS}]_{\rm full}=[{\cal O}^{\rm T}_{\rm BS}]_{\rm
on}+[{\cal O}^{\rm T}_{\rm BS}]_{\rm inst}+[{\cal O}^{\rm T}_{\rm BS}]_{\rm
z.m.}=0$. Because of this, the decay constant for the pseudotensor current
needs to be defined only through the nonlocal matrix element
$\bra{0}\bar{q}(z)\Gamma_{\rm T}q(-z)\ket{P}$ defined by Eq. (III).

Defining $z^{\mu}=\tau\eta^{\mu}$ using the lightlike vector $\eta=(1,0,0,-1)$
and multiplying $(P_{\mu}\eta_{\nu}-P_{\nu}\eta_{\mu})$ on both sides of Eq.
(III), one can rewrite Eq. (III) as [see Ref. Choi:2017uos  for more detailed
derivation]

$\displaystyle\bra{0}{\bar{q}}(\tau\eta)i(\not{P}\not{\eta}-P\cdot\eta)\gamma_{5}q(-\tau\eta)\ket{{\rm
P}(P)}$ $\displaystyle=\frac{i}{3}f_{\rm
T}\tilde{\mu}_{M}(P\cdot\eta)^{2}\int^{1}_{0}{\rm d}xe^{i\zeta\tau
P\cdot\eta}\psi_{3;{\rm P}}(x),$ (52)

where $\tilde{\mu}_{M}=\mu_{M}(1-\rho_{+}).$333 In Ref. Choi:2017uos , the
term $1-\rho_{+}$ in Eq. (53) is absent but the inclusion of this term in this
work guarantees the process-independent decay constant in the LFQM. We then
obtain

$\psi_{3;{\rm P}}(x)=-\frac{12}{f_{\rm
T}{\tilde{\mu}_{M}}}\int^{\infty}_{-\infty}\frac{d\tau}{2\pi}\int^{x}_{0}dx^{\prime}e^{-i\zeta^{\prime}\tau(P\cdot\eta)}{\mathcalligra
M}_{\rm T},$ (53)

where ${\mathcalligra M}_{\rm T}=\langle
0|{\bar{q}}(\tau\eta)i(\not{P}\not{\eta}-P\cdot\eta)\gamma_{5}q(-\tau\eta)|M(P)\rangle$
is given by the following momentum integral in the same covariant BS model as
Eq. (45)

${\mathcalligra M}_{\rm
T}=N_{c}\int\frac{d^{4}k}{(2\pi)^{4}}\frac{e^{-i\tau(p_{2}-p_{1})\cdot\eta}H_{0}}{(p^{2}_{1}-m^{2}_{1}+i\epsilon)(p^{2}_{2}-m^{2}_{2}+i\epsilon)}S_{\rm
T},$ (54)

with the trace term $S_{\rm T}={\rm
Tr}[i(\not{P}\not{\eta}-P\cdot\eta)\gamma_{5}(\not{p}_{1}+m_{1})\gamma_{5}(-\not{p}_{2}+m_{2})]$.
It is worth noting that the explicit covariant calculation of Eq. (53) is
challenging due to the nonlocal nature of the matrix element. For the LF
calculation of Eq. (53), we apply the equal LF time condition, $z^{+}=0$, and
choose the LF gauge $A^{+}=0$ so that the path-ordered gauge factor becomes
unity. We should note that the valence contribution, i.e. $[S_{\rm T}]_{\rm
val}(x^{\prime},{\bf k}_{\perp})=[S_{\rm T}]_{\rm on}+[S_{\rm T}]_{\rm inst}$,
has the same form as the one for the local current matrix element. However,
the possible zero-mode contribution $[S_{\rm T}]_{\rm z.m.}$ is different from
the one obtained in the local current case since the trace term, as well as
the vertex function, should be integrated over $x^{\prime}$ before the
integration over $x$. The nonlocal nature of the matrix element introduces a
discrepancy in the power counting of the singular term ($1/x$), which is an
essential procedure to identify any possible zero modes, compared to the local
current matrix element calculation. As a consequence, this discrepancy gives
rise to distinct zero modes in the nonlocal current case. The identification
of the zero mode in this nonlocal matrix element calculation is not yet known.
However, our “Type II” link between the covariant BS model and the LFQM
applying only to the on-mass contribution works also for the nonlocal matrix
element calculation regardless of the existence of the LF zero mode.

Thus, considering only the on-mass shell contribution to the trace term, we
obtain from Eqs. (53) and (54)

$\psi_{3;{\rm P}}(x)=-\frac{3N_{c}}{f_{\rm
T}}\int^{x}_{0}dx^{\prime}\int\frac{d^{2}{\bf
k}_{\perp}}{8\pi^{3}}\frac{\chi(x^{\prime},{\bf
k}_{\perp})}{(1-x^{\prime})}\frac{[S_{\rm T}]_{\rm
on}}{P^{+}{\tilde{\mu}}_{M}},$ (55)

where $[S_{\rm T}]_{\rm on}=4P^{+}M^{\prime}_{0}k^{\prime}_{z}$. Now, from the
normalization of $\psi_{3;{\rm P}}(x)$, i.e. $\int^{1}_{0}{\rm
d}x\;\psi_{3;{\rm P}}(x)=1$, we obtain

$f_{\rm T}=N_{c}\int^{1}_{0}{\rm d}x\int^{x}_{0}{\rm
d}x^{\prime}\int\frac{d^{2}{\bf k}_{\perp}}{8\pi^{3}}\
\frac{\chi(x,\mathbf{k}_{\perp})}{1-x}\left[{\cal O}^{\rm T}_{\rm
BS}\right]_{\rm on},$ (56)

where $\left[{\cal O}^{\rm T}_{\rm BS}\right]_{\rm
on}=-12M^{\prime}_{0}k^{\prime}_{z}/{\tilde{\mu}}_{M}$.

Finally, applying the “Type II” link given by Eq. (A) to Eqs. (48) and (56),
we obta1in our LFQM results for $(f_{\rm A},f_{\rm P},f_{\rm T})$ defined by
Eqs. (III.1) and (26). The on-mass shell BS and LFQM operators ${\cal O}_{\rm
on}$ of the three decay constants with all possible current components are
summarized in Table 7.

## References

  * (1) G. P. Lepage and S. J. Brodsky, Exclusive Processes in Perturbative Quantum Chromodynamics, Phys. Rev. D 22, 2157 (1980).
  * (2) A. V. Efremov and A. V. Radyushkin, Factorization and asymptotic behaviour of pion form factor in QCD, Phys. Lett. B 94, 245 (1980).
  * (3) V. L. Chernyak and A. R. Zhitnitsky, Asymptotic Behavior of Exclusive Processes in QCD, Phys. Rept. 112 173 (1984).
  * (4) S. J. Brodsky, H.-C. Pauli, S. S. Pinsky, Quantum chromodynamics and other field theories on the light cone, Phys. Rept. 301 299 (1998).
  * (5) P. Ball, V. M. Braun and A. Lenz, Higher-twist distribution amplitudes of the $K$ meson in QCD, JHEP 05, 004 (2006).
  * (6) A. Accardi, J. L. Albacete, M. Anselmino, N. Armesto, E. C. Aschenauer, A. Bacchetta, D. Boer, W. K. Brooks, T. Burton and N. B. Chang, et al. Electron Ion Collider: The Next QCD Frontier: Understanding the glue that binds us all, Eur. Phys. J. A 52, 268 (2016).
  * (7) R. Abdul Khalek, et al. Science Requirements and Detector Concepts for the Electron-Ion Collider: EIC Yellow Report, Nucl. Phys. A 1026, 122447 (2022).
  * (8) V. M. Braun, Higher Twists, EPJ Web Conf. 274, 01012 (2022).
  * (9) P. Kroll, K. Passek-Kumerički, Twist-3 contributions to wide-angle photoproduction of pions, Phys. Rev. D 97, 074023 (2018).
  * (10) M. Beneke, G. Buchalla, M. Neubert, and C. T. Sachrajda, QCD factorization in $B\to\pi K,\pi\pi$ decays and extraction of Wolfenstein parameters, Nucl. Phys. B 606, 245 (2001).
  * (11) P. Ball, Theoretical update of pseudoscalar meson distribution amplitudes of higher twist: the nonsinglet case, J. High Energy Phys. 01, 010 (1999).
  * (12) T. Huang, X. H. Wu, and M. Z. Zhou, Twist-3 distribution amplitude of the pion in QCD sum rules, Phys. Rev. D 70, 014013 (2004).
  * (13) S. S. Agaev, Impact of the higher twist effects on the $\gamma\gamma^{*}\to\pi^{0}$ transition form factor, Phys. Rev. D 72, 114010 (2005).; Erratum Phys. Rev. D 73, 059902 (2006).
  * (14) A. P. Bakulev, S. V. Mikhailov, and N. G. Stefanis, Tagging the pion quark structure in QCD, Phys. Rev. D 73, 056002 (2006); QCD-based pion distribution amplitudes confronting experimental data, Phys. Lett. B 508, 279 (2001).
  * (15) S. V. Mikhailov, A. V. Pimikov, and N. G. Stefanis, Endpoint behavior of the pion distribution amplitude in QCD sum rules with nonlocal condensates, Phys. Rev. D 82, 054020 (2010).
  * (16) V. Y. Petrov, M. V. Polyakov, R. Ruskov, C. Weiss, and K. Goeke, Pion and photon light-cone wave-functions from the instanton vacuum, Phys. Rev. D 59, 114018 (1999).
  * (17) S. I. Nam and H.-Ch. Kim, Twist-3 pion and kaon distribution amplitudes from the instanton vacuum with flavor SU(3) symmetry breaking, Phys. Rev. D 74, 096007 (2006).
  * (18) E. R. Arriola and W. Broniowski, Pion light-cone wave-function and pion distribution amplitude in the Nambu–Jona-Lasinio model, Phys. Rev. D 66, 094016 (2002).
  * (19) M. Praszałowicz and A. Rostworowski, Pion light cone wave-function in the nonlocal NJL model, Phys. Rev. D 64, 074003 (2001)..
  * (20) L. Chang, C. D. Roberts, and S. M. Schmidt, Light front distribution of the chiral condensate, Phys. Lett. B 727, 255 (2013).
  * (21) C. Shi, C. Chen, L. Chang, C. D. Roberts, S. M. Schmidt, and H.-S. Zong, Kaon and pion parton distribution amplitudes to twist three, Phys. Rev. D 92, 014035 (2015).
  * (22) C.-W. Hwang, Analyses of decay constants and light-cone distribution amplitudes for $s$-wave heavy meson, Phys. Rev. D 81, 114024 (2010).
  * (23) H.-M. Choi and C.-R. Ji, Distribution amplitudes and decay constants for $(\pi,K,\rho,K^{*})$ mesons in the light-front quark model, Phys. Rev. D 75, 034019 (2007).
  * (24) H. M. Choi and C. R. Ji, Self-consistent covariant description of vector meson decay constants and chirality-even quark-antiquark distribution amplitudes up to twist-3 in the light-front quark model, Phys. Rev. D 89, 033011 (2014).
  * (25) H. M. Choi and C. R. Ji, Consistency of the light-front quark model with chiral symmetry in the pseudoscalar meson analysis, Phys. Rev. D 91, 014018 (2015).
  * (26) H. M. Choi and C. R. Ji, Two-particle twist-3 distribution amplitudes of the pion and kaon in the light-front quark model, Phys. Rev. D 95 056002 (2017).
  * (27) H.-M. Choi and C.-R. Ji,  Mixing angles and electromagnetic properties of ground state pseudoscalar and vector meson nonets in the light-cone quark model,  Phys. Rev. D 59, 074015 (1999).
  * (28) H.-M. Choi and C.-R. Ji,  Light-front quark model analysis of exclusive $0^{-}\to 0^{-}$ semileptonic heavy meson decays,  Phys. Lett. B 460, 461 (1999).
  * (29) H.-M. Choi, C.-R. Ji, Z. Li, and H.-Y. Ryu,  Variational analysis of mass spectra and decay constants for ground state pseudoscalar and vector mesons in the light-front quark model,  Phys. Rev. C 92, 055203 (2015).
  * (30) N. Dhiman, H. Dahiya, C.-R. Ji, and H.-M. Choi, Twist-2 pseudoscalar and vector meson distribution amplitudes in light-front quark model with exponential-type confining potential, Phys. Rev. D 100, 014026 (2019).
  * (31) A. J. Arifi, H. M. Choi, C.-R. Ji and Y. Oh,  Mixing effects on 1S and 2S state heavy mesons in the light-front quark model,  Phys. Rev. D 106, 014009 (2022).
  * (32) W. Jaus, Semileptonic decays of $B$ and $D$ mesons in the light-front formalism, Phys. Rev. D 41, 3394 (1990).
  * (33) W. Jaus, Relativistic constituent-quark model of electroweak properties of light mesons, Phys. Rev. D 44, 2851 (1991).
  * (34) H.-Y. Cheng, C.-Y. Cheung, and C.-W. Hwang, Mesonic form factors and the Isgur-Wise function on the light front, Phys. Rev. D 55, 1559 (1997).
  * (35) H.-M. Choi, Decay constants and radiative decays of heavy mesons in light-front quark model, Phys. Rev. D 75, 073016 (2007).
  * (36) H. J. Melosh,  Quarks: Currents and constituents, Phys. Rev. D 9, 1095 (1988).
  * (37) B. Bakamjian and L. H. Thomas,  Relativistic particle dynamics. II, Phys. Rev. 92, 1300 (1953).
  * (38) B. D. Keister and W. N. Polyzou,  Relativistic Hamiltonian dynamics in nuclear and particle physics,  Adv. Nucl. Phys. 20, 225 (1991).
  * (39) B.D. Keister,  Rotational covariance and light-front current matrix elements,  Phys. Rev. D 49, 1500 (1994).
  * (40) J. P. C. B. de Melo, T. Frederico, H.W.L. Naus, P.U. Sauer, Covariance of light-front models: pair current, Nucl. Phys. A 660, 219 (1999).
  * (41) H.-M. Choi and C.-R. Ji,  Nonvanishing zero modes in the light-front current,  Phys. Rev. D 58, 071901(R) (1998).
  * (42) S. J. Brodsky, D. S. Hwang, Exact light-cone wavefunction representation of matrix elements of electroweak currents, Nucl. Phys. B 543, 239 (1999).
  * (43) W. Jaus,  Covariant analysis of the light-front quark model,  Phys. Rev. D 60, 054026 (1999).
  * (44) J.P.B.C. de Melo, T. Frederico, E. Pace, G. Salmè, Pair term in the electromagnetic current within the Front-Form dynamics: spin-0 case, Nucl. Phys. A 707, 399 (2002).
  * (45) B. L. G. Bakker, H.-M. Choi, and C.-R. Ji, Regularizing the divergent structure of light-front currents, Phys. Rev. D 63, 074014 (2001).; The vector meson form factor analysis in light-front dynamics, Phys. Rev. D 65, 116001 (2002).
  * (46) V. A. Karmanov, The wave functions of relativistic bound systems, Zh. Eksp. Teor. Fiz. 71, 399 (1976) [transl.: JETP 44, 210 (1976)].
  * (47) J. Carbonell, B. Desplanques, V. A. Karmanov, and J. F. Mathiot, Explicitly covariant light-front dynamics and relativistic few-body systems, Phys. Rep. 300, 215 (1998).
  * (48) Q. Chang, X. N. Li, X. Q. Li, F. Su and Y. D. Yang, Self-consistency and covariance of light-front quark models: testing via $P$, $V$ and $A$ meson decay constants, and $P\to P$ weak transition form factors, Phys. Rev. D 98, 114018 (2018).
  * (49) L. Chen, Y. W. Ren, L. T. Wang and Q. Chang, Form factors of $P\to T$ transition within the light-front quark models, Eur. Phys. J. C 82, 451 (2022).
  * (50) A. J. Arifi, H. M. Choi, C. R. Ji and Y. Oh, Independence of current components, polarization vectors, and reference frames in the light-front quark model analysis of meson decay constants, Phys. Rev. D 107, 053003 (2023).
  * (51) H.-M. Choi,  Self-consistent light-front quark model analysis of $B\to Dlv_{l}$ transition form factors, Phys. Rev. D 103, 073004 (2021).
  * (52) H.-Y. Cheng, C.-K. Chua, and C.-W. Hwang,  Covariant light-front approach for $s$-wave and $p$-wave mesons: Its application to decay constants and form factors,  Phys. Rev. D 69, 074025 (2004).
  * (53) H. Y. Cheng and C. K. Chua, Light-front approach for Pentaquark strong decays, JHEP 11, 072 (2004).
  * (54) H.-M. Choi and C.-R. Ji, Relations among the light-cone quark models with the invariant meson mass scheme and the model prediction of the $\eta-\eta^{\prime}$ mixing angle, Phys. Rev. D 56, 6010 (1997).
  * (55) M. Li, Y. Li, G. Chen, T. Lappi and J. P. Vary, Light-front wavefunctions of mesons by design, Eur. Phys. J. C 82, 1045 (2022).
  * (56) H. Y. Ryu, H. M. Choi and C. R. Ji, Systematic twist expansion of $(\eta_{c},\eta_{b})\to\gamma^{*}\gamma$ transition form factors in light-front quark model, Phys. Rev. D 98, 034018 (2018).
  * (57) U. Raha and H. Kohyama, Space-and Time-like Electromagnetic Kaon Form Factors, Phys. Rev. D 82, 114012 (2010).