Radiative transitions between $0^{-+}$ and $1^{--}$ heavy quarkonia on the light front

We present calculations of radiative transitions between vector and pseudoscalar quarkonia in the light-front Hamiltonian approach. The valence sector light-front wavefunctions of heavy quarkonia are obtained from the Basis Light-Front Quantization (BLFQ) approach in a holographic basis. We study the transition form factor with both the traditional"good current"$J^+$ and the transverse current $\vec J_\perp$ (in particular, $J^R=J^x+i J^y$). This allows us to investigate the role of rotational symmetry by considering vector mesons with different magnetic projections ($m_j=0,\pm 1$). We use the $m_j=0$ state of the vector meson to obtain the transition form factor, since this procedure employs the dominant spin components of the light-front wavefunctions and is more robust in practical calculations. While the $m_j=\pm 1$ states are also examined, transition form factors depend on subdominant components of the light-front wavefunctions and are less robust. Transitions between states below the open-flavor thresholds are computed, including those for excited states. Comparisons are made with the experimental measurements as well as with Lattice QCD and quark model results. In addition, we apply the transverse current to calculate the decay constant of vector mesons where we obtain consistent results using either $m_j=0$ or $m_j=1$ light-front wavefunctions. This consistency provides evidence for features of rotational symmetry within the model.


I. INTRODUCTION
Radiative transitions offer insights into the internal structure of quark-antiquark bound states through electromagnetic probes. The radiative transition between 0 −+ (pseudoscalar) and 1 −− (vector) mesons via the emission of a photon is characterized as the magnetic dipole (M1) transition. This transition mode is known to be sensitive to relativistic effects [1], especially for those between different spatial multiplets, such as η c (2S ) → J/ψ(1S ) + γ.
Heavy quarkonium is often dubbed as the "hydrogen atom" of quantum chromodynamics (QCD). It provides an ideal testing ground for various investigations to understand QCD. States below the open flavor threshold (DD for charmonium and BB for bottomonium) have very narrow widths since they cannot decay via any Okubo-Zweig-Iizuka allowed strong decay channels [2]. Electromagnetic transition rates are therefore important. Comparing the theoretical and experimental rates for radiative transitions then provides guidance to improve our understanding of the internal structure of heavy quarkonia. Recently, some of us proposed a model based on the light-front holographic QCD and the one-gluon exchange with a running coupling [3][4][5]. In this model, charmonia and bottomonia are solved as relativistic quark-antiquark bound states using Basis Light-Front Quantization (BLFQ), a nonperturbative Hamiltonian approach within light-front dynamics [6]. The model provides a reasonable description of the mass spectrum and other properties that were studied. Observables including decay constants, r.m.s. radii [7], distribution amplitudes and parton distributions as well as diffractive vector meson productions [8] have been directly calculated from the light-front wavefunctions (LFWFs), and are in reasonable agreement with experiments and with other approaches (see also Ref. [9]). Therefore, we are motivated to investigate radiative transitions within this model. On the one hand, we hope to further test the model by comparing with the existing experimental results. On the other hand, results calculated for transitions that have not yet been measured provide predictions for future experiments.
In this work, we derive the formulae for radiative transitions between 0 −+ and 1 −− mesons on the light front, using both the traditional "good current" J + = J 0 + J z and the transverse current J ⊥ = (J x , J y ). Though, in principle, these two choices should be equivalent due to Lorentz covariance, adoption of certain approximations in the model may lead to violation of the Lorentz symmetry that would be evident through inequivalent results. In nonrelativistic quantum mechanics, magnetic moments and transitions can only be extracted from the current density J = (J x , J y , J z ) rather than the charge density J 0 . Therefore one may expect that for the M1 transitions in nonrelativistic systems such as heavy quarkonia, the transverse current density J ⊥ could be better than the charge density J + . Specifically, as we will see later, the transverse current allows us to extract the transition form factor through the m j = 0 state of the vector meson, which is not accessible with J + . The calculation with the m j = 0 state provides a more robust result by employing the dominant spin components of the two mesons (spin-triplet for the vector and spin-singlet for the pseudoscalar) in the transition, while that with m j = ±1 always requires subdominant components of the LFWFs with relativistic origins. So in this work, we obtain the transition form factors with the m j = 0 states of the vector mesons through the transverse current. The results from m j = 1, though less robust, are also presented for comparison. In addition, as a cross check, we revisit the decay constants by utilizing the transverse current to compare with the previous calculation using J + in Ref. [4]. It follows that different m j components of the vector mesons are involved. This provides us with a different yet pertinent perspective to understand the degree of Lorentz symmetry manifestation in the current model. The layout of the paper is as follows. In Sec. II, we introduce the formalism and methods to calculate the M1 transition form factor on the light front. In Sec. III, we apply the formalism to heavy quarkonia in the BLFQ approach. Sec. IV presents the numerical results of transition form factors. In Sec. V, we perform the calculations of vector meson decay constants with different magnetic projections. We conclude the paper in Sec. VI.

A. Transition form factor and decay width
The matrix element for the radiative transition between a vector meson (V) with four-momentum P and polarization m j and a pseudoscalar (P) with four-momentum P via emission of a photon, can be parametrized in terms of the transition form factor V(Q 2 ) as [10], where we define Q 2 ≡ −q 2 , with q µ = P µ − P µ representing the four-momentum of the photon. m P and m V are the masses of the pseudoscalar and the vector, respectively. e σ is the polarization vector of the vector meson. J µ (x) is the current operator.
In the physical process of V → P + γ, the photon is on shell (Q 2 = 0). The transition amplitude is: where µ,λ is the polarization vector of the final-state photon with its spin projection λ = ±1. The decay width of V → P + γ follows by averaging over the initial polarization and summing over the final polarization. In the rest frame of the initial particle, it reads, The momentum of the final photon is determined by the energy-momentum conservation, | q| = (m 2 V − m 2 P )/2m V . J V = 1 is the spin of the initial vector meson. To calculate the width of P → V + γ, exchange m V and m P , and replace J V with J P = 0 for the initial pseudoscalar in Eq. (3).

B. Light-front dynamics
In principle, the Lorentz invariant function V(Q 2 ) defined in Eq. (1) can be extracted from any of the four sets of hadron matrix elements, µ = +, −, x, y (v ± = v 0 ± v z , see definitions of the light-front variables in the Appendix). However, results from different current components may be different due to violations of Lorentz symmetry introduced by the Fock sector truncation as well as by the modeling of systems. These approximations have led to extensive discussions in the literature [11][12][13][14][15]. The "+" component, known as the "good current", is typically used, together with the Drell-Yan frame (q + = 0), to avoid contributions from pair production/annihilation in vacuum. The transverse components have been shown to be consistent with the "+" component in the limit of zero momentum transfer in certain theories, such as the φ 3 theory [13] and the spin-0 two-fermion systems [15]. Another option, the "-" component, is known as the "bad current", due to its association with the zero-mode contributions.
Here, we present formulae for the transition form factor, for both J + and J ⊥ , along with different magnetic projections (m j = 0, ±1) of the vector meson. Note that when the rotational symmetry on the transverse plane is preserved, which is usually the case, using J x or J y component or even combinations of the two are equivalent. In particular, we use J R ≡ J x + iJ y to carry out the calculation in the case of the transverse current. For any transverse vector k ⊥ , which is expressed as (k x , k y ) in the Cartesian coordinate or (k ⊥ , θ) in the polar coordinate, we will write its complex form as k R ≡ k x + ik y = k ⊥ e iθ and k L ≡ k x − ik y = k ⊥ e −iθ . From the vector decomposition in Eq. (1), where we have introduced two variables z ≡ P + /P + and ∆ ⊥ = P ⊥ − z P ⊥ , which are invariant under the transverse Lorentz boost specified by the velocity vector β ⊥ , This boost is kinematic and survives the Fock space truncation, whereas the full Lorentz transformation does not. The two sets of hadron matrix elements in Eqs. (4) and (5) can be related through such a boost, P(P + , P ⊥ + P + β ⊥ )| J ⊥ |V(P + , P ⊥ + P + β ⊥ , m j ) = P(P + , P ⊥ )| J ⊥ |V(P + , P ⊥ , m j ) + β ⊥ P(P + , P ⊥ )| J + |V(P + , P ⊥ , m j ) .
By applying the above relation to Eqs. (4) and (5), we find that for m j = ±1, J + and J R should give the same V(Q 2 ). For m j = 0, on the other hand, V(Q 2 ) cannot be extracted from J + , but can be extracted from transverse currents, such as J R .
Note that for the purpose of comparison, we label the transition form factors with their corresponding current components and the m j values of the vector wavefunctions. In practice, the different prescriptions of extracting the same transition form factor could provide a test of violation of the Lorentz symmetry in the calculation. In the covariant light-front dynamics, the transition form factor is extracted from combinations of several hadron matrix elements [11].

C. Impulse approximation
In the impulse approximation, the interaction of the external current with the meson is the summation of its coupling to the quark and to the antiquark, as illustrated in Fig. 1. The vertex dressing as well as pair creation/annihilation from higher order diagrams are neglected. The hadron matrix element can be written accordingly as a sum of the quark term and the antiquark term: The current operator is defined as The electric charge e = √ 4πα EM . For quarkonium, due to the charge conjugation symmetry, the antiquark gives the same contribution as the quark to the total hadronic current. So, for our purpose, we calculate the hadron matrix element for the quark part. As such, we computeV(Q 2 ) which is related to V(Q 2 ) by The hadron matrix element can be written explicitly in terms of the convolution of LFWFs. To begin with, the valence Fock space representation of quarkonium reads: where the color index c = 1, 2, . . . , N c , and the number of quark colors N c = 3. ψ is the LFWF written in relative coordinates. x ≡ p + /P + is the longitudinal light-front momentum fraction, k ⊥ = p ⊥ − x P ⊥ is the relative transverse momentum, where p is the single-particle 4-momentum of the quark. s represents the fermion spin projection in the x − direction.
The hadron matrix element in the Drell-Yan frame (q + = 0, i.e. z = 1) follows, We could then obtain the transition form factor from such hadron matrix elements according to Eqs. (8) and (9). Ideally,V(Q 2 ) is independent of the spin projection m j and the current components. Nevertheless, one needs to carefully choose the proper matrix elements to evaluate certain quantities, when approximations break the Lorentz symmetry [16]. For instance, there are different ways of choosing matrix elements to calculate the spin-1 electromagnetic form factors when the angular condition is violated [17]. Among those, some are preferred in the sense that unphysical terms could be partially or entirely suppressed [15,18].
In the case of the M1 transition form factorV(Q 2 ), using the combination of the transverse current J R with the m j = 0 polarization of the vector meson, according to the expression in Eq. (9), would give the LFWF representation as, where we define ψ ↑↓±↓↑ ≡ (ψ ↑↓ ± ψ ↓↑ )/ √ 2. Note that in deriving Eq. (13), we have taken advantage of symmetries in LFWFs, ψ (m j =0) ↑↓−↓↑/V = 0 and ψ ↑↓+↓↑/P = 0. The traditionally used good current J + is also worth looking at. As we have discussed, with this current component, the transition form factor can be extracted only from the m j = ±1 polarizations of the vector meson. We present the expression for m j = 1 according to Eq. (8), while the expression for m j = −1 is similar. It is evident from this expression that the overlapped spin components of the two wavefunctions indicate no spin-flip (between spin-triplet and spin-singlet), which may appear counter-intuitive for the M1 transition. Indeed, this calculation relies on subdominant terms and is less robust, as we will discuss in the following section for heavy quarkonia.
D. The nonrelativistic limit In the nonrelativistic limit, the M1 transition with the same radial or angular quantum numbers (e.g. nS → nS + γ), is often referred to as allowed, for which the transition amplitude is large andV(0) → 2 as a result of the similarity between the spatial wavefunctions of the vector and the pseudoscalar mesons with the same spatial quantum numbers; whereas the transition between states with different radial or angular excitations is referred to as hindered, for which the transition amplitude is zero andV(0) → 0 at leading order due to the orthogonality of the wavefunctions [1,[19][20][21]. The deviations of experimentally measured results from those nonrelativistic limits indicate relativistic effects [22]. For a heavy quarkonium system, which is close to the nonrelativistic domain, such deviations are expected to be small but nonzero.
The wavefunctions of heavy quarkonia, treated as relativistic bound states, are dominated by those components that are non-vanishing and reduce to the nonrelativistic wavefunction in the nonrelativistic limit. These wavefunction components are therefore referred to as the dominant components. It is necessary to emphasize that despite the correspondence between the dominant spin components and the nonrelativistic wavefunctions, the former carries relativistic contributions when solved in a relativistic formalism. There are also wavefunction components of purely relativistic origin, which vanish in the nonrelativistic limit and are therefore subdominant. In practice, the dominant components tend to be better constrained, while the subdominant ones are more sensitive to the model and numerical uncertainties. For the pseudoscalar states resembling S-waves (in particular n 1 S 0 ), η c (nS ) and η b (nS ), their dominant components are the spinsinglets ψ ↑↓−↓↑/P , while relativistic treatments would also allow them to have subdominant components, such as ψ ↑↑/P = ψ * ↓↓/P . Analogously, for the vector states of heavy quarkonia resembling S-waves (in particular n 3 S 1 ), ψ(nS ) and Υ(nS ), the dominant components are the spin-triplets, ψ For those vector states identified as D-waves, ψ(n 3 D 1 ) and Υ(n 3 D 1 ), where orbital excitations occur, all the spin-triplet components ψ ↑↑/V (m = 0), could also exist in the nonrelativistic limit. Moreover, the spin components with m = 0 admit the admixtures of S-waves, though the actual amount of such admixtures is small and sensitive to both the model parameters and the truncation. For example, the ψ(3770) [ψ(1D)] state, though primarily a 1 3 D 1 state, has contributions from n 3 S 1 states (notably 2 3 S 1 ) [23][24][25][26], and these S-wave admixtures are responsible for the nD → n S + γ transitions [21,[27][28][29]. In order to have a more intuitive view of the dominant and subdominant spin components for those states, we take the LFWFs from Ref. [4] to show in Fig. 2 the proportions of those dominant and subdominant components of heavy quarkonia. For all those pseudoscalar and vector states, the dominant terms could each occupy 88% ∼ 100% of the whole LFWF, suggesting that the heavy quarkonium indeed resembles a nonrelativistic system. The comparison between the same states of the charmonium and those of the bottomonium also reveals that the dominant component is more pronounced in the heavier, and less relativistic, system.
The nonrelativistic limit can be achieved forV m j =0 (Q 2 ) by adopting nonrelativistic wavefunctions where only the dominant spin components exist. However, withV m j =1 (Q 2 ), simply taking the nonrelativistic wavefunction would always lead to zero since the expression in Eq. (14) involves the vanishing subdominant terms. To be specific, we examine the transition form factors at Q 2 = 0, where they can be interpreted as the overlaps of wavefunctions in coordinate space [ψ (m j ) ss ( r ⊥ , x)], shown in Eqs. (15) and (16). Though equivalent to Eqs. (13) and (14) at Q 2 = 0 respectively, Eqs. (15) and (16)    of 1/q R , and are therefore more intuitive for the purpose of illustration.
Note that in the nonrelativistic limit, the wavefunctions of the respective pseudoscalar and vector states with the same radial and angular numbers are identical in their spatial dependence, and they only differ in their spin structures. For the allowed transition,V m j =0 (Q 2 = 0) → 2 due to the normalization of the spatial wavefunctions, which can be seen from Eq. (15) along with taking x → 1/2 + k z /(2m q ) [4] and small hyperfine splitting m P ≈ m V . For the hindered transition,V m j =0 (Q 2 = 0) → 0 due to the orthogonality of the two wavefunctions. Such a nonrelativistic reduction that takes advantage of the near orthonormality of wavefunctions is reminiscent of the nonrelativistic quark model (see Refs. [1,19,20]). However, forV m j =1 (Q 2 ), the realization of the nonrelativistic limits depends strongly on the details of the subdominant wavefunctions which are less constrained in the parameter fitting. For the hindered transition, where cancellation occurs, this leads to a strong sensitivity to the model parameter and potentially to the truncation. Fig. 3 presents the integrands (those inside {. . .}) ofV m j =0 (0) andV m j =1 (0) according to Eqs. (15) and (16) for an allowed (1S → 1S + γ) as well as a hindered (2S → 1S + γ) transition. In the left panel of Fig. 3, the integrands of the allowed transition have no nodes resulting from the coherent overlaps of the two wavefunctions. On the other hand, the right panel of Fig. 3 shows significant cancellations of contributions from the integrands which change sign due to nodes in the 2S wavefunctions. Based on these lines of reasoning, we takeV(Q 2 ) =V m j =0 (Q 2 ), using the transverse current, to evaluate the transition form factors for heavy quarkonia. The less robustV m j =1 (Q 2 ), using the plus current, which has strong sensitivity to the violation of rotational symmetry will also be presented for comparison.

III. CALCULATION IN BASIS LIGHT-FRONT QUANTIZATION
We adopt wavefunctions of heavy quarkonia from recent work [3,4] in the BLFQ approach [6]. The effective Hamiltonian extends the holographic QCD [30] by introducing the one-gluon exchange interaction with a running coupling [31]. The starting point for the model of Refs. [3,4] is transverse light-front holography, inspired by string theory, which approximates QCD at long distance. As a complementary part, the longitudinal confining potential is introduced to allow a more consistent treatment of the mass term and the longitudinal excitation. The one-gluon exchange implements the short-distance physics and determines the spin structure of the mesons. The mass spectrum and LFWFs are the direct solutions of the eigenvalue equation, and are obtained by diagonalizing the Hamiltonian in a basis representation. The spectrum agrees with the PDG data with an r.m.s mass deviation of 30 to 40 MeV for states below the open flavor thresholds. The LFWFs have been used to produce several observables and are in reasonable agreement with experiments and other theoretical approaches [32]. We now use these same LFWFs to calculate radiative transitions with the formalism described in Sec. II.
The LFWFs are solved in the valence Fock sector using a basis function representation: In the transverse direction, the 2D harmonic oscillator (HO) function φ nm is adopted as the basis. In the longitudinal direction, we use the modified Jacobi polynomial χ l as the basis. m is the orbital angular momentum projection, related to the total angular momentum projection as m j = m + s +s, which is conserved by the Hamiltonian. The basis space is truncated by their reference energies in dimensionless units: As such, the N max -truncation provides a natural pair of UV and IR cutoffs: Λ uv κ √ N max , λ ir κ/ √ N max , where κ is the oscillator basis energy scale parameter. L max represents the resolution of the basis in the longitudinal direction. See Ref. [4] for details on basis functions and parameter values. The LFWFs are calculated at N max = L max = 8, 16, 24 and 32. Transition form factors are computed at each of these basis truncations. Fig. 4 shows the convergence trends ofV(0) as functions of 1/N max . The left panel compares three different fitting functions to extrapolate our results obtained at finite basis sizes to the complete basis by taking the limit N max = L max → ∞. The extrapolations using these three functions agree to within 1% of each other. For the remainder of this paper we adopt the second-order polynomial in 1/N max for our extrapolations, fitted to the 4 basis sizes, with an extrapolation uncertainty given by the difference between the result in the largest basis (N max = L max = 32) and the extrapolated value. (Note that this uncertainty does not include any systematic uncertainty coming from the model for the interaction or from the Fock space truncation.) In the right panel of Fig. 4 we show our results for all allowed transitions at finite basis sizes, together with our extrapolation to the complete basis, including our extrapolation uncertainty estimate. Note that they are all close to the nonrelativistic limit 2, which is expected according to our discussion in Sec. II D. For the hindered transitions, the uncertainties from such basis extrapolations are comparatively larger, since their calculations are more sensitive to the details of wavefunctions.

IV. RESULTS
Here we present our results for selected pseudoscalar-vector transition form factors for charmonia and bottomonia below their respective open flavor thresholds. Fig. 5 shows our numerical results of the transition form factors in three groups, the allowed transition nS → nS + γ, the radial excited transition nS → n S + γ (n n ) and the angular excited transition nD → n S + γ, through a progression from upper to lower panels. As already discussed in the previous section, for the allowed transitions we findV(0) ≈ 2, whereas for the hindered transitions involving either radial or angular excitations we haveV(0) ≈ 0. Transitions involving higher radial excited states feature more wiggles in the curve, which is especially evident in the nS → nS + γ transitions as n increases. This is because the radial excited states have transverse nodes. As a result, the transition form factors, in the form of their convolutions [see Eq. (13)], are not monotonic. The comparison between charmonia and bottomonia is also of interest. For comparable transition modes, the transition form factors show similarity in their patterns as well as their behaviour as a function of N max . Furthermore, as illustrated in the second row of panels in Fig. 5, one observes that the P(nS ) → V(n S ) + γ transition form factors are very similar to the V(nS ) → P(n S ) + γ form factors for n > n .
Comparisons ofV(0) from this work, with experiments (compiled by PDG [22]) and with other models (Lattice QCD [33][34][35][36][37], Quark Model [38][39][40]) are collected in Table I and visualized in Fig. 6. Most calculations, as well as available experimental data, give a value of the the order of 2 for the allowed transitions nS → nS + γ: all such data in Table I are between 1.5 and 2.5 with only one exception, the relativistic quark model calculation of J/ψ → η c + γ. This is in agreement with the vector V(nS ) and the pseudoscalar P(nS ) mesons possessing very similar spatial wavefunctions, but different spin structures. On the other hand, there is a significant spread in the theoretical results of the hindered transitions. This is expected because the hindered transitions involve changes in radial quantum numbers and/or orbital angular motions and are sensitive to delicate cancellations as discussed above. Considering the fact that only two free parameters are employed by the model for quarkonia in Ref. [4] and the fact that we did not adjust any parameters in our calculation for the transitions, the agreement to within an order of magnitude is encouraging.
The results with the "+" current, in combination with the m j = 1 vector meson wavefunctions,V m j =1 (0), are presented as a ratio toV m j =0 (0) in Fig. 7. As already mentioned, we expect these calculations to be much less robust. This is because the calculation ofV m j =1 (Q 2 ) according to Eq. (14), depends on subdominant components of the wavefunctions, which is less constrained from the model. Indeed, the dependence of these calculations on the basis truncation is much larger, resulting in significantly larger extrapolation uncertainties. Furthermore, the hindered transitions have a much larger fluctuation than the allowed transition, due to their sensitivity to the subdominant components in one of the two spatial wavefunctions with different radial quantum numbers and/or different orbital motions. Our results with the "+" component of the currentV m j =1 (0), differ by up to 2 orders of magnitude from our more reliable results with the transverse component of the currentV m j =0 (0).

V. DECAY CONSTANTS
In the discussion of the transition form factors above, we saw that differences could arise when different magnetic projections of the vector mesons are used, in combination with different components of the electromagnetic current operator. These differences are linked with violations of rotational symmetry in the model. We argue that the violation of rotational symmetry is not a major factor by checking two representative observables. The first one is the meson masses. It has been shown in Ref. [4] that the mass eigenvalues associated with different magnetic projections are in reasonable agreement, with a mean mass spread of 17 MeV (8 MeV) for charmonia (bottomonia). The second observable is the decay constant, which we present in this section. The decay constant is defined with the same current operator as the transition form factor. For a vector meson, its value, too, could be extracted from different magnetic projections. On the other hand, it features the simplicity of involving only one LFWF instead of convoluting two LFWFs. Therefore, the decay constant provides a pathway for examining the rotational symmetry of LFWFs.
Decay constants for vector mesons are defined as the local vacuum-to-hadron matrix elements:    . The first row shows the allowed transitions, the second row shows transitions between different radial excitations, and the third row presents those involving angular excitations. The dashed and solid curves are calculated with LFWFs at N max = L max = 8 and N max = L max = 32 respectively. The shaded areas in between indicates the numerical uncertainty from basis truncation. As a consequence of the UV cutoff from the basis, the largest Q 2 ( Λ 2 uv ) at N max = 32 truncation is 31 GeV 2 (44 GeV 2 ) for charmonia (bottomonia).
In Ref. [4], the "+" current is used. Since e + (P, m j = ±1) = 0, with the "+" component of the current the decay constant can only be extracted from m j = 0, using However, in analogy to the transition form factor, we can also use the transverse current. In the case of m j = 0, we get exactly the same expression for the decay constant as with the "+" current, namely Eq. (20). TABLE I.V(0) for radiative decay between 0 −+ and 1 −− charmonia (bottomonia) below the DD (BB) threshold. Values from PDG [22] are converted from their decay widths according to Eq. (3). Note that the uncertainties of meson masses propagate into that ofV(0). The BLFQ results are from Eq. (13). For these results, all meson masses are taken from PDG [22], except that Υ(1 3 D 1 ), Υ(2 3 D 1 ) and η b (3S ) masses are taken from Ref. [4]. Extrapolations for BLFQ are made from N max = L max = 8, 16, 24, 32 using second-order polynomials in N −1 max . We use the difference between the extrapolated and the N max = 32 results to quantify numerical uncertainty, which does not include any systematic uncertainty. Uncertainties are quoted in parenthesis and apply to the least significant figures of the result. Some lattice results are quoted with more than one source of uncertainty. The lattice nonrelativistic QCD (NRQCD) [41] results are converted from their three-point matrix elements with meson masses from PDG [22]. Values from the relativistic quark model (rQM) [38] and the Godfrey-Isgur (GI) model [39,40] are converted from their decay widths according to Eq. (3) with their suggested meson masses, respectively. These results are plotted in Fig. 6.  (10) 0.025(13) 0.035 0.084 Υ(3S ) → η b (2S ) + γ < 0.167 0.145 (33) 0.056 0.65 This should not come as a surprise. Recall that the "+" and the transverse matrix elements with the same m j can be related through a transverse Lorentz boost [see Eq. (7)]. Furthermore, with the transverse current we can also calculate the decay constant from the m j = ±1 components of the vector meson. The expression for the decay constant with m j = 1 follows as Here m q is the quark mass which is one of the model parameters determined in Ref. [4]. Note that using m j = −1 would lead to an equivalent expression considering the symmetry between the m j = ±1 LFWFs.  Table. I. Extrapolations are made from N max = L max = 8, 16, 24, 32 using second-order polynomials in N −1 max . We use the difference between the extrapolated and the N max = 32 results to quantify numerical uncertainty which is indicated by the vertical error bars on the BLFQ results (sometimes smaller than the symbols). We do not include any systematic uncertainty. Quarkonia in the initial and final states are labeled on the top and bottom of the figure. Single (double) apostrophe stands for the radial excited 2S (3S) state. The D-wave states are identified as n 3 D 1 . The heavy quark limitV(0) = 2 of the allowed (nS → nS + γ) transition is shown in the dashed line. Results from PDG [22], Lattice QCD [33][34][35][36] and Lattice NRQCD [37,41], the relativistic quark model (rQM) [38] and the Godfrey-Isgur (GI) model [39,40] are also presented for comparison.
Ratio ofV m j =1 (0) toV m j =0 (0), calculated from Eq. (14) and Eq. (13) respectively. BLFQ results are extrapolated to N max = ∞ from N max = L max = 8, 16, 24, 32 using second-order polynomials in N −1 max . We use the difference between the extrapolated and the N max = 32 results to quantify the numerical uncertainty (indicated by vertical error bars), which is dominated by the uncertainty in the m j = 1 result obtained with the "+" component of the current. The nS → nS + γ transitions are shown as filled triangles, whereas the hindered transitions, involving radial/angular excitations, are shown as open diamonds.
As is the transition form factor, the decay constant is also Lorentz invariant and thus it should be independent of the polarization m j . Therefore in practice, the difference between f V (m j = 0) and f V (m j = 1) provides another measure of the violation of rotational symmetry by our model. For vector meson states identified as S-wave states, both f V (m j = 0) and f V (m j = 1) arise primarily from the dominant spin components of LFWFs, which relate to the nonrelativistic wavefunctions. Moreover, the two expressions reduce to the same form in the nonrelativistic limit, where x → 1/2 + k z /(2m q ) [4] and m q → m V /2. That is, For vector meson states identified as Dwave states, both f V (m j = 0) and f V (m j = 1) are calculated mainly from the dominant but less occupied spin components of LFWFs (see Sec. II D for discussions of spin components for D-wave states). The two resulting small values are sensitive to model and numerical uncertainties, and could reveal differences. Nevertheless, for the S-wave states of heavy quarkonia, we expect to find robust results when either m j = 0 or m j = 1 is used to calculate f V . The parameters m V and m q involved in f V (m j = 1) might result in additional uncertainty but the resulting fluctuation should be small for heavy systems.
We use both Eqs. (20) and (21) to calculate the decay constants for the lowest three charmonium and five bottomonium vector states below their open flavor thresholds. The results are presented in Fig. 8, where basis truncations are chosen to match the UV cutoffs Λ uv κ √ N max ≈ 1.7m q [4]. The two sets of results, f V (m j = 0) and f V (m j = 1), are all within each others' extrapolation uncertainties for the five S-wave states in Fig. 8. Such consistency implies that the rotational symmetry is reasonably preserved in LFWFs. For the D-wave states, the decay constants are small but differences between f V (m j = 0) and f V (m j = 1) can be noticed since each result depends on different small components. By implication, the transition form factor V(Q 2 ) m j =0 calculated using the transverse current, with its overlapping dominant components of both the initial and the final LFWFs, is further supported as a robust result.

VI. SUMMARY AND DISCUSSIONS
We have derived formulae for the radiative transitions between vector meson and pseudoscalar mesons (M1 transitions) in light-front dynamics. We then calculated the transition form factors for heavy quarkonia obtained in the BLFQ approach. The majority of the predictions from the BLFQ approach are in reasonable agreement with experimental data and other model calculations when the m j = 0 state of the vector meson and the transverse current J R = J x + iJ y are employed. We have also shown that at least in the context of heavy quarkonium, this choice is preferred to the traditional good current J + (cf. Ref. [11]). The latter lacks the access to the m j = 0 state of the vector mesons, and generates a less robust result. Moreover, in extracting the decay constants of vector mesons, J + gives the same result as J R when the m j = 0 state is used, but can not be used when the m j = 1 state is taken.
The comparison of results extracted from different m j states of the vector mesons also allowed us to investigate the rotational symmetry of the model. The consistency of decay constants, f V (m j = 0) and f V (m j = 1), suggests that the rotational symmetry is achieved at a reasonable level with the adopted LFWFs for states identified as S-waves. By contrast, the M1 transition form factors,V(Q 2 ) m j =0 andV(Q 2 ) m j =1 , as convolutions of the initial and final states LFWFs, reveal details of the violation of symmetry between m j = 0 and m j = 1 LFWFs, even between S-wave states. In such circumstance,V(Q 2 ) m j =0 is found to be more robust due to overlapping dominant parts of LFWFs, and it became our suggested approach in heavy systems. However, for states identified as D-waves, our results for both the transition form factors and for the decay constants show large difference between calculations using m j = 0 and calculations with m j = 1. This suggests that our m j = 0 and m j = 1 LFWF for these states have noticeably different contributions beyond the D-waves, such as those from the S-waves (which can be present even in the nonrelativistic calculations) and those with purely relativistic origin. In addition, they could also have different components that vanish when rotational symmetry is restored. The breaking of rotational symmetry in the wavefunction mainly results from the Fock sector truncation and basis truncation. We expect to improve rotational symmetry by including higher Fock sectors, as shown in some simpler theories [42].
In future work, we will also extend the current formalism to electric dipole transitions (E1 transitions) and higher multipole contributions such as M2. With a variety of transition modes, we hope to establish a more comprehensive understanding on the internal structure of quark-antiquark bound states and to further address the role of Lorentz symmetry. Radiative transitions within other meson sectors, such as the light meson system and the mixed flavor systems, are also of great interest. Furthermore, besides the currently used Drell-Yan frame, we could also adopt various other frames where the transition form factor in the timelike region may be accessible. The frame dependence study could also serve as a metric for the violation of the Lorentz symmetry, as shown for the elastic form factors [16].