$E1$ and M1 radiative transitions involving heavy-light axial, pseudoscalar and vector quarkonia in the framework of Bethe-Salpeter equation

This work is an extension of our previous work in \cite{bhatnagar20} to calculate M1 transitions, $0^{-+}\rightarrow 1^{--} \gamma$, and E1 transitions involving axial vector mesons such as, $1^{+-} \rightarrow 0^{-+}\gamma$, and $0^{-+}\rightarrow 1^{+-} \gamma $ for which very little data is available as of now. We make use of the general structure of the transition amplitude, $M_{fi}$ derived in our previous work \cite{bhatnagar20} as a linear superposition of terms involving all possible combinations of $++$, and $--$ components of Salpeter wave functions of final and initial hadrons. In the present work, we make use of leading Dirac structures in the hadronic Bethe-Salpeter wave functions of the involved hadrons, which makes the formulation more rigorous. We evaluate the decay widths for both the above mentioned $M1$ and $E1$ transitions. We have used algebraic forms of Salpeter wave functions obtained through analytic solutions of mass spectral equations for ground and excited states of $1^{--}$,$0^{-+}$ and $1^{+-}$ heavy-light quarkonia in approximate harmonic oscillator basis to do analytic calculations of their decay widths. We have compared our results with experimental data, where ever available, and other models.


Introduction
One of the challenging areas in hadronic physics is probing the inner structure of hadrons. Charmonium occupies an intermediate regime between the bb system and the light mesons. Radiative decays of charmonia are good testing grounds for various models, due to the fact that emitted photon can be directly detected, and electromagnetic interactions are well understood.
Radiative transitions characterized by ∆L = 0 are the magnetic dipole, M1 transitions, while those characterized by |∆L| = 1 are the electric dipole, E1 transitions. The M1 transition mode is sensitive to relativistic effects, specially between different spatial multiplets (where n > n ), while the E1 transitions are much stronger than M1 transitions, and involve transitions between excited states. In this work, besides the M1 transitions, we study the E1 transitions involving 1 +− mesons, which are the P-wave states.
The P wave cc states were first observed in 1976 by the SLAC-LBL experiment at SLAC/SPEAR [1,2], where they observed the decay, Ψ(3684) → γ + χ c . The P wave bb states were first observed by the Columbia-Stony Brook (CUSB) experiment at the Cornell CESR electron-positron storage ring [3,4] and confirmed by the CLEO experiment at CESR [5].
The 1 ++ mesons are seen in pp collisions. However, not many decays of these mesons are experimentally observed as can be checked from PDG tables [11,12].
We wish to mention that decay rates of M1 transitions are much weaker than the rates for E1 transitions. But M1 decay rates are interesting as they allow access to spin-singlet states, that are very difficult to produce. Thus, as regards M1 decays, we study decays, η c (2S) → J/Ψ(1S)γ, and B c (2S) → B * c (1S)γ. B c meson was discovered in pp collisions at √ s= 1.8 TeV. using CDF detector at Fermilab tevatron. It is the only heavy meson with two heavy quarks with different flavours that forbid their decays into two photons.
We calculate the radiative decay widths of these heavy-light quarkonia for the above mentioned processes in the framework of 4×4 Bethe-Salpeter equation (BSE), which is a fully relativistic approach that incorporates the relativistic effect of quark spins and can also describe internal motion of constituent quarks within the hadron in a relativistically consistent manner, due to its covariant structure [29,30]. Our wave functions satisfy the 3D BSE, which is in turn obtained from 3D reduction of the 4D BSE under Covariant Instantaneous Ansatz (which is a Lorentz-invariant generalization of Instantaneous Approximation), already has relativistic effects. Further, our transition amplitudes also have relativistically covariant form.
The present work, where we make use of two leading Dirac structures in the structure of BS wave functions of and A − (1 +− ) quarkonia involved in these radiative transitions is more rigorous than our previous work in [26], where we made use of only the most leading Dirac structure (γ 5 for P-mesons, iγ. for V-meson, and I for S-meson) in the BS wave functions of the hadrons involved in the processes.
Our previous studies on mass spectral calculations of heavy-light quarkonia [31,32] were used to fit the input parameters of our model as C 0 = 0.69, ω 0 = 0.22 GeV, Λ QCD = 0.250 GeV, and A 0 = 0.01, with input quark masses m u = 0.300 GeV, m s = 0.430 GeV,m c = 1.490 GeV, and m b = 4.690 GeV. In the present work on radiative decays, we use these same input parameters to calculate the single photon decay widths for the above processes.
Now, as mentioned in our previous works [30][31][32][33][34], we are not only interested in studying the mass spectrum of hadrons, which no doubt is an important element to study dynamics of hadrons, but also the hadronic wave functions that play an important role in the calculation of decay constants, form factors, structure functions etc. for QQ, and Qq hadrons. These hadronic Bethe-Salpeter wave functions were calculated algebraically by us in [30][31][32]. The plots of these wave functions [32] show that they can provide information not only about the long distance non-perturbative physics, but also act as a bridge between the long distance, and short distance physics, and are provide us information about the contribution of the short ranged coulomb interactions in the mass spectral calculation of heavy-light quarkonia. These wave functions can also lead to studies on a number of processes involving QQ, and Qq states, and provide a guide for future experiments.
This paper is organized as follows: In section 2, we give the general formulation of the process, H → H + γ in the framework of 4 × 4 Bethe-Salpeter equation under the covariant instantaneous ansatz. In sections 3, we calculate the single photon decay widths for the processes, P → V γ, In Section 4, we deal with the processs, A − → P γ. In Section 5, we deal with the process, P → A − γ, where, P, A − , and V are the pseudoscalar, axial vector, and vector heavy-light quarkonium states. In section 6, we give numerical results and discussions.

Instantaneous Ansatz
The Bethe-Salpeter equation that describes the quark-anti quark bound state of momenta p 1 and p 2 , relative momentum q, and meson momentum P is where K(q, q ) is the interaction kernel, and S −1 F (±p 1,2 ) = ±i p 1,2 + m 1,2 are the quark and antiquark propagators. We now make use of the Covariant Instantaneous Ansatz (which is a Lorentz-covariant 3D support), where, K(q, q ) = K( q, q ) on the BS kernel, where the BS kernel depends entirely on the variable, q µ = q µ − q.P P 2 P µ is the component of internal momentum of the hadron that is orthogonal to the total hadron momentum, i.e. q.P = 0, while σP µ = q.P P 2 P µ is the component of q longitudinal to P , where the 4-dimensional volume element is, d 4 q = d 3 qM dσ. Now working on the right side of Eq.(1), and making use of the fact that and the fact that the longitudinal component of M dσ of q does not appear in K(q,q ), carrying out integration over M dσ on right side of Eq.(1), we obtain, where, Γ(q) is the hadron-quark vertex function, and is directly related to the 4D wave function, Ψ(P, q), and one can express the 4D BS wave function Ψ(P, q) in terms of Γ(q) as, Further the 4D hadron-quark vertex, that enters into the definition of the 4D BS wave function in the previous equation, can be identified as, Following a sequence of steps outlined in [30], we get four Salpeter equations which are effective 3D forms of BSE (Salpeter equations) given below: Thus, in our framework, a crucial role is played by the component,q µ , which is always orthogonal to P µ and satisfies the unconstrained relation,q.P = 0, regardless of whether q.P = 0 (i.e. σ = 0), or q.P = 0 (i.e. σ = 0).
The Lorentz-invariant nature ofq 2 increases the applicability of this framework of Covariant Instantaneous Ansatz all the way from low energy spectra to high energy transition amplitudes. For details, see [26,33].
The 3D B.S. wave function can be expressed in terms of the projected wave functions as where and the projection operators with the relation In radiative transitions involving single photon decays, such as H → H + γ, the process requires calculation of triangle quark-loop diagram, which involves two hadron-quark vertices that we attempt in the 4 × 4 representation of BSE. The single photon decay of Qq quarkonia is described by the direct and exchange Feynman diagrams as in Figure   1.

Figure 1: Radiative decays of heavy-light quarkonia
These processes involve two Lorentz frames, the rest frame of the initial meson, H, and the rest frame of final meson, H . Let P , and q be the total momentum and the internal momentum of initial hadron, while P , and q the corresponding variables of the final hadron. And let k, and λ be momentum and polarization vectors of emitted photon, while λ is the polarization vector of final emitted meson. Thus if p 1,2 , and p 1,2 are the momenta of the two quarks in initial and final hadron respectively, then, we have, the momentum relations: for initial and final hadrons respectively. From the Feynman diagrams we see that conservation of momentum demands that, P = P + k. Now, for the first diagram, we have the kinematical relations, p 1 = p 1 + k, and −p 2 = −p 2 , where k = P − P is the momentum of the emitted photon. And for the second diagram, we have the corresponding relations, p 1 = p 1 , and −p 2 = −p 2 + k.
Making use of the above equations, the relationship between the internal momenta, q, and q of the initial and final hadrons in terms of the photon momentum, k can be expressed as, where, the first equation is for Diagram 1, and second equation is for Diagram 2. Here,m 1, the Wightman-Garding definitions [33] of masses of individual quarks, which act like momentum partitioning functions for the two quarks in a hadron.
For initial hadron, its internal momentum has already been decomposed as, q = (q, iM σ), whereq, and M σ are defined relative to its external momentum, P . Similarly for final meson, we again decompose its internal momentum, q into two components q = (q , iM σ ), withq = q − σ P transverse to initial hadron momentum, P , and σ = q .P P 2 , longitudinal to P . Thus, P.q = 0. The relationship between the transverse components of internal momenta of the two hadrons,q, andq is [26] where, the first equation of Eq.(13), holds for Diagram 1, and the second equation holds for Diagram 2, andP is the component of total momentum P of final hadron transverse to initial hadron momentum, P . Here,q .P = 0 due tô P .P = 0. Now, the kinematics gets simplified in the rest frame of the initial meson, where we have P = ( and since the photon momentum can be decomposed as, . Thus the energy of the emitted meson can be expressed as, The dot products of momenta of the initial and the emitted meson can be expressed as [26], Thus, it can be seen that, −E acts as the projection of P along the direction of initial hadron momentum, P .
Similarly, the relationship between the longitudinal components of internal momenta of the two hadrons in the two diagrams can be worked out as [26], which is again a consequence of the transversality ofq with initial hadron momentum, P . Thus, up to Eq.(15), the kinematics is the same for all the three processes (P → V γ, A − → P γ, and P → A − γ) studied in this work.

M1 Radiative decays, P → V γ
In this section we study the M1 radiative decay process, P → V γ, having studied the process, V → P γ in our previous work [26]. In the present study, we make use of two leading Dirac structures (as in Eq. (35)) in the BS wave functions of the two hadrons involved in the process, in contrast to a single most leading Dirac structure (γ 5 for P meson, and iγ. for V meson) used in our previous work [26]. This makes the calculations more involved than with use of the single most leading Dirac stricture in the 4D hadronic BS wave functions. It is to be noted that 4D BS wave function of initial pseudoscalar meson involved in the process is Ψ P (P, q) = S F (p 1 )Γ P (q)S F (−p 2 ). Since for transition amplitude calculation, we choose to do calculation in the rest frame of the initial pseudoscalar meson, we write the wave function of the emitted vector meson as, where as defined earlier,q = q − q .P P 2 P is transverse to initial hadron momentum, P , where the hadron-quark vertex function Γ V (q ) for the final meson is, Similarly for final meson, the expression for the 3D BS wave function, ψ(q ) is expressed in terms of the projection operators as, where with the relation ω 2 1,2 = m 2 1,2 +q 2 . In Eq.(12), due to relationship between q , and q, we can express the electromagnetic transition amplitude of the process as a single integral over d 4 q = d 3q M dσ as, Here, the first term corresponds to the first diagram, where the photon is emitted from the quark, while the second term corresponds to the second diagram where the photon is emitted from the antiquark, and e q , and e Q are the electric charge of quark, and antiquark respectively, while λ µ is the polarization vector of the emitted photon. We can express Eq. (20) as, Now, we reduce the above equation to the effective 3D form by integrating over the longitudinal component, M dσ over the poles of the propagators, S F (p i ), that are expressed as [26], ; where the rest of the terms are anticipated to be zero on account of 3D Salpeter equations. The contour integrations over M dσ are performed over each of the four terms taking into account the pole positions in the complex σ plane: The contour integral over each of the four terms can be performed by closing the contour either above or below the real axis in the complex σ-plane.
We now make use of the Salpeter equations in variable,q in Eq. (6), and the Salpeter equations in variableq given below. It is to be noted that the Salpeter equations inq involve −E = P.P M , which is the projection of P along the direction of initial momentum, P , and are given as [35], It can be verified that the results of each of these four integrals, Ω 1 , ..., Ω 4 whether we close the contour above or below the real σ-axis comes out to be the same, thereby validating the correctness of the formalism employed. These results of integrals over dσ in Ω 1 , ..., Ω 4 , are given as, α 1 , ..., α 4 in Eqs. (32).
This leads to the expression for effective 3D form of transition amplitude, M 1 f i under Covariant Instantaneous Ansatz for Diagram 1, as, Here, the results of contour integrals over dσ are given as α 1 , ..., α 4 : and the projected wave functions, ψ ±± being taken from the 3D Salpeter equations [32] derived earlier, which for initial meson in internal variableq are given in Eq. (6).
[It is to be noted that the factors (M ± ω 1 ± ω 2 ) that were also present in the numerators of α's in Eq. (29) Here, the results of contour integrals over dσ are given as α 1 , ..., α 4 : , Thus, we make use of the generalized method for handling quark-triangle diagrams with two hadron -quark vertices in the framework of 4 × 4 BSE under Covariant Instantaneous Ansatz described in [26], Now, to calculate the amplitude, M f i for the process, we need the ++, and −− components, ψ ±± V (q ) for vector and ψ ±± P (q) for pseudoscalar mesons, given as, to calculate which we need the 3D wave functions, ψ P (q), and ψ V (q ). To derive these, we start with the general 4D decomposition of BS wave functions [36,37]. Using 3D decomposition under Covariant Instantaneous Ansatz, the wave function of vector mesons of dimensionality, M can be written as [26,30,31], where λ is the vector meson polarization vector, while for a pseudoscalar meson, the 3D wave function with dimensionality, M can be written as [26,30,31], However, in calculation of hadronic observables, it was noticed that some Dirac covariants [37] in structure of hadronic BS wave function covariants contribute much more than others. And in accordance with a naive power counting rule in [33,38,39], for pseudoscalar mesons, one could classify Dirac structures, M γ 5 , and P γ 5 associated with amplitudes φ 1 , and φ 2 respectively as leading, while those with φ 3 , and φ 4 as sub-leading. A similar behaviour was observed in case of vector mesons [33,38], where Dirac structures, M and P associated with χ 1 and χ 2 respectively are leading [33,38], while those associated with χ 3 , ..., χ 6 are sub-leading. A similar observation about the most leading Dirac structures from all the Dirac structures was made by [40,41]. Thus to simplify the algebra, we take 3D wave functions with these two leading Dirac structures (the present work with two leading Dirac structures is more involved than our previous work [26], where we took only the most leading Dirac structure in hadronic BS wave functions of P, V, and S mesons [26]), with the 3D radial wave functions, φ P (q), and φ V (q), obtained as solutions of mass spectral equations [30][31][32] of pseudoscalar, and vector quarkonia respectively, that were obtained from 3D Salpeter equations being, where β P,V [32] are the inverse range parameters.
These radial wave functions were recently used for calculations of some M1 and E1 transitions in [26]. We have made use of the 3D Salpeter equations in Eq. (6), that depend on the variableq 2 , that is Lorentz-invariant, and is a four-scalar, whose validity extends over the entire 4D space, while also keeping contact with the surface, P.q = 0 (hadron rest frame). Now, our mass spectrum, and the 3D wave functions φ(q) in Eq. (36) and (49) (please see [30]) were calculated from Salpeter equations in Eq.(6) in the rest frame of the hadron.
The Bethe-Salpeter normalizers, N P , and N V are obtained through current conservation condition, The transition amplitude, M f i is expressed as, where the antisymmetric tensor, µναβ ensures its gauge invariance. Here, F P V is the transition form factor for P → V γ, with expression, where the expressions for F 1 P V and F 1 P V are given in Eqs.(80)-(81) in the Appendix, and correspond to contributions from the two diagrams, which is a more accurate treatment. Here [T R 1 ], and [T R 2 ] involve trace over gamma matrices.
The above expression corresponds to F P V (k 2 = 0), that corresponds to emission of a real photon.
The kinematical relation connectingq 2 , withq 2 , is given in Eq. (37). To calculate the decay widths, we need to calculate the spin averaged amplitude square,  ). Thus, |M f i | 2 can be expressed as, The decay width of the process (P → V γ) in the rest frame of the initial meson is expressed as where we make use of the fact that modulus of the momentum of the emitted vector meson can be expressed in terms of masses of particles as, where, ω k is the kinematically allowed energy of the emitted photon. Thus, Γ in turn can be expressed as: The M1 decay widths for the processes studied are given in Table 2, for the meson masses given in Table 1, that were evaluated in our previous works [28,32] We now calculate the radiative decay widths for the process, A − → P + γ in the next section.

Radiative decays of heavy-light quarkonia through A − → P γ
The E1 transitions are characterized by |∆L| = 1. Thus in these transitions, there is change in parity between the initial and final hadronic states. The scattering amplitude of the decay process A − → P γ corresponding to Diagram 1, can be written as,   Here, the results of contour integrals over dσ are given as α 1 , ..., α 4 in Eq. (29). Similarly the amplitude for Diagram 2 can be written as, And, the results of contour integrals over dσ are given as α 1 , ..., α 4 in Eq. (31). After the 3D reduction of the 4D BS wave function of scalar meson under CIA, we express the 3D BS wave function with dimensionality M as Making use of the fact that the leading order Dirac structures in axial meson BS wave function in accordance with the power counting rule proposed in [33,38,39] are γ 5 .q, and iγ 5 .q P M , and making use of [31], we express the 3D axial meson BS wave function of dimensionality, M as, where φ A − (q) is the spatial part of this wave function, whose analytic forms obtained by power series solutions of 3D mass spectral equations (derived from 3D Salpeter equations in Eq. (6)), in the variableq (which is in fact |q|) for P-wave meson in its own rest frame, calculated in [31] are These wave functions in Eq.(49) involve even powers ofq, along with odd power,q. Hereq = |q| as explained above is the length of the 3D vectorq, and is expressed as |q| = q 2 − (q.P ) 2 /P 2 , and is a Lorentz-invariant quantity [26], along with even powers ofq, such asq 2 ,q 4 ,... etc. which are again Lorentz-invariant. While for S-wave mesons, the wave functions are only functions of even powers ofq. Thus when P-wave or S-wave mesons are in the final state, their wave functions after Lorentz transformation, would involve the variables, |q |, and/or even powers ofq . We expressq 2 in terms ofq 2 directly through Eq. (37), that connectsq 2 withq 2 , while odd power, |q | is expressed as, |q | = q 2 , where we again make use of Eq. (37). Thus, the time component, σ of q will not appear in the wave functions of both P-wave and S-wave mesons in final state in the transition amplitude calculation. The BS normalizer of axial meson, N A − , can be obtained by solving the current conservation condition in Eq. (38).
We now obtain the ++ and −− components of the axial meson wave function with the corresponding adjoint wave functions as in case of P and V mesons done earlier, and calculate expressions for + + ++, + + −−, − − ++, and − − −− terms that appear in the scattering amplitude in Eq. (46). We then evaluate trace over the gamma matrices in Eq. (46). We make use of the fact thatq =q +m 2P , where, P = P − P .P P 2 P . We combine various terms, and further make use of the fact that, for initial axial meson, P. λ = 0. We can expressP . = P . . We can then express the invariant matrix element, M 1 f i corresponding to Diagram 1 as, Similarly for Diagram 2, we write the invariant matrix element, M 2 f i as, Total amplitude M f i for the process can be expressed as the sum of amplitudes from the two diagrams, where integrals over d 3q are performed over Θ i , and Θ i .
Thus, R 1 , R 2 , and R 3 are the three form factors in the above equation. Now, with a change of variables, and making use of the electromagnetic gauge invariance, it can be shown that these three form factors appearing in M f i , are not independent, and we can express the invariant amplitude, M f i in terms of a single form factor.
To show this, we introduce two new external variables,P , and k, which are defined as, whereP is the sum of momenta of initial and emitted mesons, and k is the emitted photon momentum. Thus, we can express the initial and final meson momenta P and P in terms of new variables, as P =P +k 2 , and P =P −k 2 . In terms of these new variables, we can express the amplitude, M f i as, Now, the transversality property of polarization vector of axial meson, P. = 0 leads to, We now introduce a new form factor, S = −S 2 + S 3 , in terms of which we can write M f i as, Now due to Eq.(55), the term with S 2 vanishes. Thus, amplitude, M f i can be expressed as, Now, electromagnetic gauge invariance demands, k µ M µν = 0. This leads to the relation, between the form factors, which is like an equation of constraint. Thus, due to the electromagnetic gauge invariance, S 1 , and S are no longer independent, and we can express the amplitude, M f i in terms of a single form factor, S 1 , whose expression is given in the next equations.
Now, to calculate the decay widths, we need to calculate the spin averaged amplitude square, the spin-averaged amplitude square of the process can be written as We can write the decay width, where we make use of the fact that modulus of the momentum of the emitted pseudoscalar meson can be expressed in terms of masses of particles as, | − → P | = 1 2M (M 2 − M 2 ).

Radiative decays of heavy-light quarkonia through P → A − γ
We proceed to evaluate the process, P → A − γ using Fig.1. Here, the initial pseudoscalar (0 −+ ) meson decays into an axial vector (1 +− ) meson and a photon. We can then express the effective 3D form of invariant matrix element, M 1 f i corresponding to Diagram 1 as, Similarly the amplitude for Diagram 2 can be written as, Taking trace over the gamma matrices in the previous equations, we get, the invariant matrix element, M 1 f i for Diagram 1 as, Similarly for Diagram 2, we write the invariant matrix element, M 2 f i as, Total amplitude M f i for the process can be expressed as the sum of amplitudes from the two diagrams, The structure of M f i above is similar to [25].
To calculate the decay widths, we again need to calculate the spin averaged amplitude square, where we sum over the final polarization states, λ of photon, and λ of V-meson. Following a similar procedure as in A − − > P γ, it is seen that the contributions to spin averaged amplitude modulus square arises only from S 1 , while the contributions from S 2 , and S 3 vanish after doing the averaging over the ploarization states.

Results and Discussion
The present work is an extension of our work in [26] to study of radiative M1 decays, P → V γ, and E1 decays, A − → P γ, In our previous work [26], we had simplified the calculation by considering only the most leading Dirac structures in the wave functions of P, V and A mesons, that contribute maximum to calculation of all meson observables in accordance with our power counting rule [33,38]. However in the present calculation, we consider the two leading order Dirac structures in the BS wave functions of P, and A mesons given in Eqs. (35), and (48) in accordance with the power counting scheme we proposed in [33,38], which makes this calculation more rigorous.
Using this generalized expression for M f i , in Eq.(28-31), we have evaluated the decay widths for M 1 transitions, We wish to mention that as seen from Tables 2 and 3, the decay rates of M1 transitions are much weaker than the rates for E1 transitions. But M1 decay rates are interesting as they allow access to spin-singlet states, that are very difficult to produce. However, the known M1 decay rates show a serious disagreement between theory and experiment, as can be seen from Table 2.
Similarly we again see a wide range of variations in different models for both M 1, and E1 transitions, particularly for decays of η c , and h c mesons. Further, our decay widths for nS− > n S transitions in M1 decays, and nP − > n S transitions in E1 decays show a marked decrease as we go from ground to higher excited states, which is in conformity with data and other models. We have also given our predictions for radiative decays, h c (2P )− > η c (2S)γ, h c (2P )− > η c (1S)γ, and η c (2S) → h c (1P ) for which data is not yet available. As regards M1 transitions, we have given our prediction for the decay width of η c (2S) → J/Ψ(1S)γ, for which the PDG tables [12] give only the upper limit on the decay width. Also we calculated the decay width B c (2S) → B * c (1S)γ for which data is not available. The aim of doing this work was mainly to study the processes, P → V γ, A − → P γ, and P → A − γ for which very little data is available. This study was also to test the algebraic forms of wave functions of A − mesons [28] that we have recently derived from the mass spectral equations of these axial mesons, along with the wave functions of P mesons and V mesons by studying their transitions. These wave functions were obtained as solutions of their mass spectral equations in an approximate harmonic oscillator basis obtained analytically from 4 × 4 BSE as a starting point, that has so far given good predictions [30][31][32] not only of the mass spectrum of heavy-light quarkonia, but also their leptonic decays, two-photon, and two gluon decays. The present work would in turn lead to the validation of our approach, which provides a much deeper insight than the purely numerical calculations in 4 × 4 BSE approach that are prevalent in the literature.
A more detailed study on the transition form factors of both M 1, and E1 transitions, but also the "static" form factors describing meson-photon interactions through the vertex M γM for various mesons will be relegated to a separate paper.