Radiative decays D ∗ ( s ) → D ( s ) γ in covariant conﬁned quark model

Radiative decays D ∗ ( s ) → D ( s ) γ are revisited in light of new experimental data from the BABAR and BESIII collaborations. The radiative couplings g D ∗ Dγ encoding nonperturbative QCD eﬀects are calculated in the framework of the covariant conﬁned quark model developed by us. We compare our results with other theoretical studies and experimental data. The couplings (in GeV − 1 ) | g D ∗ + D + γ | = 0 . 45(9) and | g D ∗ 0 D 0 γ | = 1 . 72(34) calculated in our model agree with the corresponding experimental data | g D ∗ + D + γ | = 0 . 47(7) and | g D ∗ 0 D 0 γ | = 1 . 77(16). The most interesting case is the decay D ∗ s → D s γ , for which a recent prediction based on light-cone sum rules at next-to-leading order | g D ∗ s D s γ | = 0 . 60(19) deviates from the ﬁrst (and only to date) lattice QCD result | g D ∗ s D s γ | = 0 . 11(2) at nearly 3 σ . Our calculation yields | g D ∗ s D s γ | = 0 . 29(6), which falls somehow between the two mentioned results, although it is larger than those predicted in other studies using quark models or QCD sum rules.


I Introduction
Decays of D * mesons are dominated by the radiative D * → Dγ (γ emission) and strong D * → Dπ (π emission) processes.The hadronic effects in these decays are encoded in the magnetic moment coupling g D * Dγ and axial coupling g D * Dπ , respectively.Owing to experimental difficulties, data concerning their width are rare.The first measurement of the width of a vector D meson was achieved by the CLEO collaboration for the charged D * + meson.They observed Γ(D * + ) = (96 ± 4 ± 22) keV [1], whose uncertainties are statistical and systematic.In 2013, the BABAR collaboration obtained the width with much higher precision: Γ(D * + ) = (83.3± 1.2 ± 1.4) keV [2].In 2014, the BESIII collaboration reported the measurements of the branching fractions of the D * 0 decays to be B(D * 0 → D 0 π 0 ) = (65.5 ± 0.8 ± 0.5)% and B(D * 0 → D 0 γ) = (34.5 ± 0.8 ± 0.5)%, assuming the dominance of these two channels [3].These are the most precise measurements of the branching fractions to date.Recently, the BESIII collaboration has announced the measurement of the branching of D * + s → D + s π 0 relative to D * + s → D + s γ to be (6.16 ± 0.43 ± 0.19)% [4].This measurement confirmed the previous BABAR's result and made it more accurate [5].It should be noted that the decay D * + s → D + s π 0 violates isospin symmetry and is highly suppressed with respect to the dominating mode D * + s → D + s γ.The mechanism of D * + s → D + s π 0 is still an open question.However, it is evidently not the same π emission as for D * + → (Dπ) + and D * 0 → D 0 π 0 .It can be concluded that, in light of these new experimental data, as well as data to be reported in the near future from BESIII, the study of D * -meson radiative decays is fruitful and promising.
Radiative decays of D * meson can shed light on the dynamics of many hadronic processes.
In this paper, we revisit the radiative decays of the D * mesons in the framework of the covariant confined quark model (CCQM) developed previously by our group.This is a model based on quantum field theory in which the interaction between the constituent quarks and corresponding hadronic bound state is described by a Lagrangian in the form of a hadron field coupled to a quark current.Our model has the advantages of treating multiquark states in a similar manner and obtaining observables in the full physical range without any extrapolation.We aimed to provide predictions for both D * and D * s decays; nevertheless, we further focused on the decay D * s → D s γ, for which theoretical predictions still vary significantly and experimental data on its width are expected to be available soon.
The rest of the paper is organized as follows.In Sec.II, we briefly introduce the covariant confined quark model.The calculation of the radiative decays D * (s) → Dγ using our model is demonstrated in Sec.III.Numerical results are presented in Sec.IV, where we also provide a detailed comparison of our results with experimental data and other theoretical studies.
Finally, a brief summary is given in Sec.V.

II Covariant Confined Quark Model
The covariant confined quark model (CCQM) is a model based on quantum field theory.
In this model, hadronic bound states are described by quantum fields that interact with their constituent quarks via an interaction Lagrangian.In the case of a meson M(q 1 q2 ), this Lagrangian has the form where Γ M is the relevant Dirac matrix and g M is the meson-quark coupling constant.The vertex function F M (x; x 1 , x 2 ) effectively represents the finite size of the meson and relative quark-hadron position.The vertex function has to obey the identity F M (x+a; x 1 +a, x 2 +a) = F M (x; x 1 , x 2 ) for any given four-vector a so as to satisfy the translational invariance.Its form is chosen as where ω i = m q i /(m q 1 + m q 2 ) with m q i being the constituent quark mass; therefore, 1.This chosen form satisfies the translational invariance mentioned above.Besides, it corresponds to an effective description of a meson.The Dirac delta function effectively represents the barycenter of the two-quark system.The function Φ M [(x 1 − x 2 ) 2 ] depends on the distance between the two quarks, which can be understood as the effective size of the meson.The Fourier transform of the function Φ M [(x 1 − x 2 ) 2 ], which we denote as Φ M (−k 2 ), can be calculated, in principle, from the solutions of the Bethe-Salpeter equation for meson bound states.In a series of previous papers of us [26][27][28], after trying various forms for the vertex function, we found that the hadronic observables are insensitive to the details of the functional form of the quark-hadron vertex function.This observation is used as a guiding principle, and the function Φ M [(x 1 − x 2 ) 2 ] is assumed to be Gaussian for simplicity; it is expressed in terms of the momentum representation as The size parameter Λ M of the meson M is one of the free parameters of the model.The minus sign in the argument of the function Φ M (−k 2 ) emphasizes that we are working in Minkowski space.Given that k 2 turns into −k 2 in the Euclidean space, the Gaussian function exp(k 2 /Λ 2 M ) falls off appropriately in the Euclidean region.We stress again that any choice for Φ M is applicable provided it has the sufficiently fast fall-off behavior in the ultraviolet region of the Euclidean space to render the ultraviolet finite of Feynman diagrams.
The normalization of particle-quark vertices is provided by the compositeness condition [29,30] where Z M is the wave function renormalization constant of the meson M and Π ′ M is the derivative of the meson mass function taken on the mass-shell p 2 = m 2 M .The square root of the renormalization constant Z 1/2 M can be interpreted as the matrix element between the physical state and corresponding bare state.When Z M = 0, the physical state does not contain the bare state; therefore, it is described as a bound state.The Lagrangian given in Eq. (1) describes the interaction of the constituents (quarks) and corresponding physical state (meson), which is the bound state of the constituents.As a result of the interaction, the physical state is dressed, and its mass and wave function have to be renormalized.
Note that, even before the rise of QCD, it was well understood that describing composite particles such as hadrons within QFT is not easy.In QFT, free fields are quantized by imposing (anti)commutator relations between creation and annihilation operators, which act on vacuum and help construct the asymptotic in-and out-states.Physical processes are described by the S-matrix elements, which are convolutions of propagators.The Lagrangian that describes free fields and their interactions must be renormalized; or in other words, bare or unrenormalized quantities such as wave functions, coupling constants, and mass must be renormalized so that they represent physical (renormalized) ones.When describing a bound state, e.g., the decays of a hadron in a physical process, physical (renormalized) quantities such as its mass must be used.The relation between the bare field and the dressed one can be expressed using the wave function renormalization constant as φ 0 = Z 1/2 φ r .Note that the bare field φ 0 can be eliminated from the Lagrangian by setting its wave function renormalization constant to be zero, i.e., Z = 0. Jouvet [31] suggested that the equation Z = 0 can be used as a compositeness condition.In particular, he proved the equivalence between a four-fermion theory and a Yukawa-type theory if the renormalization constant of the boson field is set to zero.Then, the renormalized boson mass and Yukawa coupling can be expressed in terms of the Fermi constant via the compositeness condition.A more detailed discussion of the compositeness condition is given in the Appendix.
The compositeness condition effectively excludes the constituent degrees of freedom from the physical state space and avoids double counting.In other words, the constituents exist only in virtual states.This point is illustrated in a more intuitive manner in Fig. 1.
The meson mass function in Eq. ( 4) is defined by the Feynman diagram shown in Fig.  and has the following form for a pseudoscalar meson and vector meson, respectively: where we use the free quark propagator It is convenient to calculate the derivatives of the meson mass functions by using the following identities: Then, the derivatives of the meson mass functions have the form The Fock-Schwinger representation of quark propagators is then used to calculate the loop integrations in Eqs. ( 9) and ( 10): As will be seen later, the Fock-Schwinger representation constitutes an efficient approach for the calculation of tensor loop integrals because loop momenta can be converted into derivatives of the exponent function.
Loop integrations are performed in the Euclidean space.The transition from the Minkowski space to the Euclidean space is done by the Wick rotation Simultaneously, all external momenta are rotated, i.e., p 0 → ip 4 , so that p 2 = −p 2 E ≤ 0. Note that the quadratic form in Eq. ( 11) now becomes positive-definite, and the integral over α is absolutely convergent.The Minkowski notation will be kept to avoid excessive relabeling; keep in mind that k 2 ≤ 0 and p 2 ≤ 0.
Using the representations of the vertex functions and quark propagators given by Eqs. ( 3) and (11), respectively, we calculate the Gaussian integration in the derivatives of the mass functions in Eqs. ( 9) and (10).The exponent has the form ak 2 + 2kr + z 0 , where r = b p.
Next, with the help of the properties (k is the loop momentum) we replace k by ∂ r = γ µ ∂ ∂rµ , which allows us to exchange the tensor integrations for a differentiation of the Gaussian exponent.For example, Eq. ( 6) now has the form The r-dependent Gaussian exponent e −r 2 /a is moved to the left through the differential operator ∂ r by using the properties Finally, the derivatives are moved to the right by using the commutation relation The last step is done by using a form code that works for any numbers of loops and propagators.
The remaining integrals over the Fock-Schwinger parameters 0 ≤ α i < ∞ are treated by introducing an additional integration that converts the set of these parameters into a simplex as follows: As an example, we show below the expression for the derivative of the mass function of a pseudoscalar meson: The function f P (t, α) arises from trace evaluation.
At this stage, an infrared cutoff is introduced to avoid any possible thresholds in the The infrared cutoff parameter λ effectively guarantees the confinement of quarks within hadrons.This method is generic and can be used for diagrams with an arbitrary number of loops and propagators.In the CCQM, the parameter λ is assumed to be universal.The numerical evaluation of the final integral in Eq. ( 19) was done by FORTRAN code using the NAG library.
Let us enumerate the number of free parameters in the CCQM and describe the fitting process to determine their values.For a given meson M i , five parameters are needed: the size in the CCQM as where N c = 3 is the number of colors and O µ = γ µ (1 − γ 5 ) is the weak Dirac matrix with left chirality.The mesons are taken on their mass-shells.The matrix elements in Eq. ( 20) are calculated in a similar manner as that for the case of the mass functions described above.
In addition, we decided to use eight fundamental electromagnetic decays to establish further constraints on the model parameters.These decays are listed in Table II.The results of (overconstrained) least-squares fitting for the leptonic decay constants and electromagnetic decay widths are reported in Tables I and II.The results of the fitting for the values of quark masses, size parameters, and infrared cutoff parameter are reported in Tables III and   IV.
Once the free parameters are fixed, the CCQM can be employed as a frame-independent tool for hadronic calculation.Note that this data fitting was performed in 2011 (see Ref. [38]).One of the advantages of the CCQM is the possibility to calculate hadronic quantities in the full physical range of momentum transfer without any extrapolation.Another advantage of the model is that it can be used to treat not only mesons [39][40][41][42], but also baryons [43][44][45], tetraquarks [46,47], and other multiquark states [48,49] in a consistent way.
Concerning the estimation of theoretical errors in the CCQM, we use MINUIT for fitting, which is based on χ 2 minimalization.When we performed the fitting, we obtained the bestfitting values for model parameters, which are listed in Tables III and IV.We did not consider the errors of these parameters because it is complicated to take into account the error propagations in the final physical predictions.Instead, we followed an easier but less 1.58 ± 0.37 rigorous approach for error estimation.We only used the best-fitting parameters for further calculation.Then, we observed that the fitted values deviated from experimental data by approximately 5% − 10%.In our model, hadronic quantities are calculated in a similar manner to the calculations of the leptonic and electromagnetic decay constants.Therefore, we estimated the errors of the hadronic quantities to be approximately 5% − 10%.When these hadronic quantities are used for further calculations, for example, the decay widths, the errors accumulate.Consequently, we estimated the error for the decay widths to be approximately 10% − 20%.In particular, in this study, we estimated the error for the couplings g D * Dγ to be 10% and the error for the decay widths to be 20%.
The radiative decays D * (s) → D (s) γ are described by the Feynman diagrams shown in Fig. 4. The coupling between the quarks and the photon is described by the interaction Lagrangian, where e c and e q are the quark charges in units of e.
The transition amplitude is expressed as where The ratios of quark masses now read ω 1 = m c /(m c + m q ) and ω 2 = m q /(m c + m q ) with q = u, d, s.
To calculate the matrix element M(D * → Dγ), the Gaussian form is substituted for the vertex functions in Eq. (3) and Fock-Schwinger representation of quark propagators in Eg. (11), and then the techniques described by Eqs. ( 12)-( 19) are applied.
Next, using the transversality conditions ǫ γ µ (q)q µ = 0 and ǫ D * ν (p)p ν = 0, the matrix element can be expressed in the form where The expression for I q (m 2 D * , m 2 D ) can be obtained by simply exchanging m c ↔ m q and ω 1 ↔ ω 2 , i.e.

IV Numerical results and discussion
Our results for the radiative decay constants are listed in Table V together with results obtained in other theoretical approaches.It is worth mentioning that the sign difference in the predictions in Table V is simply due to the choice when one defines the matrix element and does not affect the physical observables, i.e., the decay width, given that the last is proportional to the square of the matrix element.In other words, the sign of the radiative decay constants cannot be measured experimentally.The experimental value g D * + D + γ = −0.47(7)GeV −1 was extracted from the total width of D * + and measured branching fraction [25].Meanwhile, the total width of D * 0 is still unknown experimentally.
Note that our predictions g D * + D + γ = −0.45(9)GeV −1 and g D * 0 D 0 γ = 1.72(34)GeV −1 are in full agreement with the current experimental data.Note also that the decay width of the channel D * + → D + γ was included in the model parameter fitting (see Ref. [38] and Table II).However, channels D * 0 → D 0 γ and D * + s → D + s γ were not included in the fit.In the case of g D * + s D + s γ , the theoretical predictions vary significantly with respect to the case of g D * + D + γ and g D * 0 D 0 γ .The central value of our result, |g D * + s D + s γ | = 0.29(6) GeV −1 , is larger than all other theoretical predictions, except for a recent value, |g D * + s D + s γ | = 0.60(19) GeV −1 , obtained by using LCSR(NLO) [25].Our result deviates from others within 2σ owing to large uncertainties in the predictions.In Ref. [25], a discussion on D-spin symmetry for the decays D * + → D + γ was also provided.Here, the D-spin symmetry is understood as the symmetry under the exchange of down and strange quarks d ↔ s.As pointed out in Ref. [25], this is a good approximate symmetry in QED.Given that the decay D * + s → D + s γ can be obtained from the decay D * + → D + γ by exchanging d ↔ s, the order of this symmetry breaking is expected not to be too far from that in QED.It was roughly estimated that the deviation between g D * + D + γ and g D * + s D + s γ is approximately 20-30%.Our results show a deviation of approximately 35% which is reasonable.Note also from Table V the difference between the predictions of Refs.[25] and [10], despite the similarity in the approach, which is LCSR at NLO.The authors of Ref. [25] explained this issue as the result of computational differences between both studies.In particular, the twist-1 O(α s )-corrections were calculated in Ref. [25] while they were neglected in Ref. [10], and the twist-4 corrections were dismissed in Ref. [25] while they were considered in Ref. [10].Experiment [54] 1.33(36) 19.9(5.5) - Table VI compares the predictions for the width of the D * (s) radiative decays.Available theoretical predictions for Γ(D * + → D + γ) vary in the range (0.2 − 1.7) keV.Our result Γ(D * + → D + γ) = 1.21 (48) agrees well with experimental data, suggesting the width of the decay D * + → D + γ to be at the higher end of the predicted range.Note that the experimental value Γ(D * + → D + γ) = 1.33 (36) keV was obtained from the world average of the total width Γ(D * + ) = 83.4(1.8)keV and branching fraction B(D * + → D + γ) = 0.016(4) reported by Particle Data Group (PDG) [54].We stress again that the decay D * + → D + γ was included in our model parameter fitting (see Table II).However, this fitting was performed in 2011 [38], and the input value Γ(D * + → D + γ) = 1.5 (5) keV was taken from PDG at that time [33].The slight difference for the central values of Γ(D * + → D + γ) given in Tables VI and II, which are 1.21 and 1.22, respectively, is because we used the most updated values for meson masses in this study.
It follows from Eq. ( 28) that At the final step, we have neglected the mass differences between charged and neutral D ( * ) mesons.Note also that we have assumed from the beginning that m d = m u .In the heavy quark limit (HQL) m c → ∞, the integral I c is suppressed as 1/m c and the width ratio yields 0.25.However, from Table VI, one obtains Γ(D * + → D + γ)/Γ(D * 0 → D 0 γ) = 0.066.
Therefore, the use of HQL is not suitable for the decays D * → Dγ.This is opposed to the case of the decays B * → Bγ, in which the use of HQL for b quark is reliable, as pointed out in a previous paper of us [55].
In the case of the decay D * s → D s γ, there is still much speculation about its width.Predictions for Γ(D * s → D s γ) range from 0.07 keV to 2.4 keV.Our result Γ(D * s → D s γ) = 0.55 (22) keV deviates from the recent prediction based on LCSR at NLO Γ(D * s → D s γ) = 2.36 +1. 49 −1.41 keV [25] and the first LQCD result Γ(D * s → D s γ) = 0.066 (26) keV [24] by approximately 2σ.Note also that the central value of our result Γ(D * s → D s γ) = 0.55 (22) keV is larger than that of most predictions; in particular, it exceeds the LQCD central value Γ(D * s → D s γ) = 0.066 keV [24] by roughly an order of magnitude.By contrast, the recent prediction based on LCSR at NLO Γ(D * s → D s γ) = 2.36 s → e + ν e ), likely to be obtained in the near future, will help clarify the value of the decay width Γ(D * + s → D + s γ).Given that our predictions for the case of D * + and D * 0 fully agree with the experimental values, we wanted to check if the discrepancy in the case of D * s is due to the value of the constituent s-quark mass.We varied the mass of s-quark by approximately ±20% with respect to its fitted value m s = 0.424 GeV while keeping other parameters unchanged.The dependence of the width Γ(D * s → D s γ) on m s is shown in Fig. 5.Note that our prediction is not consistent with the experimental data within the errors.More accurate experimental data are needed in the case of D * s radiative decay.If this discrepancy still remains, we will have to refit our model parameters in the future.
Finally, we compared the contribution from the charm quark and the other one (up, down, and strange quark) into the radiative coupling constants.As mentioned above, we have that where the expression of I c(q) is given in Eq. ( 27).The electromagnetic aspect of the decays is mostly visible in the proportionality of the quark contribution to its electric charge.Given that e c = e u = +2/3 and e s = e d = −1/3, there is a large cancellation between the contributions from the charm and down quarks in g D * + D + γ , as well as from the charm and strange quarks in g D * + s Dsγ .Meanwhile, for g D * 0 D 0 γ , the contributions from the charm and up quarks add up (see Table V).The strong effects in the decays are captured in the quantities I c(q) .Table VII compares I c and I q , (q = u, d, s) by taking the ratios I q /I c .This table also includes the numerical values for the ratios e q I q /e c I c .

TABLE VII. Comparison of quark contributions to radiative decay constant
Channel Quark q I q /I c e q I q /e c I c g D * Dγ (GeV s γ, our result for the decay width 0.55 (22) keV is smaller than the recent LCSR prediction 2.36 +1. 49 −1.41 keV by more than 1σ but exceeds the first LQCD prediction 0.066 (26) keV by more than 2σ.Using the branching fraction of the purely leptonic decay s → e + ν e recently measured by the BESIII collaboration, we estimate Γ(D * + s → D + s γ) to be approximately 0.05(3) keV.Meanwhile, theoretical predictions for this width range from 0.07 keV to 2.4 keV, and most of them are around 0.2 keV.In particular, this estimated value is smaller than our prediction by more that 2σ.Therefore, the decay D * + s → D + s γ is an interesting case that requires further theoretical and experimental studies to give a final Comparing the renormalized generating functionals of the Yukawa theory in Eq. (A9) and the Fermi theory in Eq. (A16), it can be concluded that the condition for their equality is Therefore, the vanishing of the wave function renormalization constant in the Yukawa theory can be used as the condition for the bare (unrenormalized) field φ 0 = Z 1/2 φ r to vanish for a composite boson.
D ) is the radiative decay constant.The quantities I c(q) (m 2 D * , m 2 D ) are threefold integrals which are calculated numerically.The expression for I c (m 2 D * , m 2 D ) reads FIG. 5. Γ(D * s → D s γ) vs. strange-quark mass.
s → D + s γ s 2.97 −1.49−0.29(6) V Summary Radiative decays of the vector mesons D * + , D * 0 , and D * + s are studied in the framework of the Covariant Confined Quark Model in light of new experimental data from the BABAR and BESIII collaborations.Predictions for the radiative couplings g D * + D + γ , g D * 0 D 0 γ , and g D * + s D + s γ , as well as the decay widths are reported.A full agreement between these predictions and experimental data was found for the decays D * + → D + γ and D * 0 → D 0 γ.In the case of D * + s → D + answer.If the condition G ΠS=0 (m 2 ) = 1 is required and the boson field is rescaled as φ → φΠ′ S=0 (m 2 ), one obtains the free Lagrangian of the boson field with mass m and the correct residue of the Green function.The renormalized generating functional of the Fermi theory finally readsZ ren F = Dφ exp i 2 dx φ(x)(✷ − m 2 )φ(x)

TABLE I .
Input values for the leptonic decay constants f M (in MeV) and our least-squares fitting [38]EII.Input values for some basic electromagnetic decay widths and our least-squares fitting values (in keV)[38].

TABLE V .
Radiative decay constants in the CCQM and other approaches (all in GeV −1 )

TABLE VI .
Decay widths Γ(D * (s) → D (s) γ) in the CCQM and other approaches (all in keV) [56]9−1.41keV[25]has a central value approximately four times as large as our value. Otheoretical studies seem to agree onthe central value Γ(D * s → D s γ) ≈ 0.2 keV.It is instructive to estimate the width of the decay D * s → D s γ using the recent measurement of the leptonic decay D * + s → e + ν e performed by the BESIII collaboration[56].This is the first measurement of the decay branching fraction; it reads B(D * + s → e + ν e ) = (2.1 +1.2 −0.9stat.± 0.2 syst. ) × 10 −5 .By taking the ratio of branching fractions B(D * + s → e + ν e )/B(D + s → µ + ν µ ) one obtains = 46(30) eV, which is much smaller than most theoretical predictions listed in Table VI, including ours.It should be noted that this estimation is based on the first measured branching fraction B(D * + s → e + ν e ) = (2.1 +1.2 −0.9stat.± 0.2 syst. ) × 10 −5 which still suffers from large uncertainties.A more precise measurement of B(D * +