Radiative decays ${{\boldsymbol{D^{*}}}_{\boldsymbol{(s)}} {{\bf\to}\boldsymbol{ D}_{\boldsymbol{(s)}}\boldsymbol\gamma}}$ in covariant confined quark model

Decays of $ D^* $ mesons are dominated by the radiative $ D^*\to D\gamma $ (γ emission) and strong $ D^*\to D\pi $ (π emission) processes. The hadronic effects in these decays are encoded in the magnetic moment coupling $ g_{D^*D\gamma} $ and axial coupling $ g_{D^*D\pi} $, 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 $\Gamma(D^{*+})=(96\pm 4\pm 22)\,\; \rm{keV}$ [1], whose uncertainties are statistical and systematic. In 2013, the BaBar Collaboration obtained the width with much higher precision: $\Gamma(D^{*+})=(83.3\pm 1.2\pm 1.4)\, \;\rm{keV}$ [2]. In 2014, the BESIII Collaboration reported the measurements of the branching fractions of the $ D^{*0} $ decays to be $ \mathcal{B}(D^{*0}\to D^0\pi^0)=(65.5\pm 0.8\pm 0.5){\text{%}} $ and $\mathcal{B}(D^{*0}\to D^0\gamma) = (34.5\pm 0.8\pm 0.5){\text{%}}$, 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^{*+}\to D_s^+\pi^0 $ relative to $ D_s^{*+}\to D_s^+\gamma $ to be $ (6.16\pm 0.43\pm 0.19){\text{%}} $ [4]. This measurement confirmed the previous BaBar Collaboration's result and made it more accurate [5]. It should be noted that the decay $ D_s^{*+}\to D_s^+\pi^0 $ violates isospin symmetry and is highly suppressed with respect to the dominating mode $ D_s^{*+}\to D_s^+\gamma $. The mechanism of $ D_s^{*+}\to D_s^+\pi^0 $ is still an open question. However, it is evidently not the same π emission as for $ D^{*+}\to (D\pi)^+ $ and $ D^{*0}\to D^0\pi^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. They allow testing various non-perturbative approaches in QCD. This explains the large amount of theoretical attempts to calculate these decays, especially prior to CLEO's first measurement of $ \Gamma(D^*) $. Radiative decays $ D^*\to D\gamma $ have been studied in the framework of QCD sum rules (QCDSR) [6−9], light cone sum rules (LCSR) [10, 11], and various quark models [12−22]. Lattice QCD (LQCD) calculation has been carried out for decays $ D^*\to D\gamma $ [23] and $ D^*_s\to D_s\gamma $ [24]. However, the theoretical uncertainties in these calculations are still significant. Recently, the first complete next-to-leading order (NLO) computation of the decays $ D^*\to D\gamma $ using LCSR has been reported [25]. In Ref. [25], the authors pointed out some differences between the first LQCD prediction [24] of $g_{D_s^*D_s\gamma}= 0.11(2) ~ \rm{GeV}^{-1}$ and the current experimental value $g_{D^*D\gamma}= -0.47(7)\,\; \rm{GeV}^{-1}$, assuming the level of D-spin breaking to be no more than 20%−30%. Moreover, the LQCD result [24] predicted the decay width to be $ \Gamma(D_s^*\to D_s\gamma)= $ 0.066(26) keV. This value disagrees with the result obtained from the light front quark model $ \Gamma(D_s^*\to D_s\gamma)= $ 0.18(1) keV [21] by $ 4\sigma $, and with the result from QCDSR $ \Gamma(D_s^*\to D_s\gamma)=0.25(8) $ keV [7] by $ 2\sigma $. The central value of this prediction is one order of magnitude smaller than most of other predictions available in the literature. Given that $ D^*_s\to D_s\gamma $ dominates the decay of the $ D^*_s $ meson, it is interesting to examine this case more carefully.

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\to D_s\gamma $, 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)}\to D\gamma $ 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.

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\bar{q}_2) $, this Lagrangian has the form

where $ \Gamma_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 $ \omega_{i}=m_{q_i}/(m_{q_1}+m_{q_2}) $ with $ m_{q_i} $ being the constituent quark mass; therefore, $ \omega_1+\omega_2=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 $ \Phi_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 $ \Phi_M[(x_1-x_2)^2] $, which we denote as $ \widetilde{\Phi}_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−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 $ \Phi_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 $ \Lambda_M $ of the meson M is one of the free parameters of the model. The minus sign in the argument of the function $ \widetilde{\Phi}_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/\Lambda^2_M) $ falls off appropriately in the Euclidean region. We stress again that any choice for $ \widetilde{\Phi}_M $ is applicable provided that 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 $ \widetilde{\Pi}'_M $ is the derivative of the meson mass function taken on the mass-shell $ p^2=m_M^2 $. The square root of the renormalization constant $ Z_M^{1/2} $ 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; 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 $ \phi_0 = Z^{1/2}\phi_r $. Note that the bare field $ \phi_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.

Figure 1.  (color online) Compositeness condition: how it works.

The meson mass function in Eq. (4) is defined by the Feynman diagram shown in Fig. 2 and has the following form for a pseudoscalar meson and vector meson, respectively:

Figure 2.  One-loop self-energy diagram for a meson.

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

so that $ k^2=k_0^2-\vec{k}^2=-k_4^2-\vec{k}^2=-k_E^2 \leq 0. $ Simultaneously, all external momenta are rotated, i.e., $p_0 \to {\rm i}p_4$, so that $ p^2=-p_E^2 \leq 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 \leq 0 $ and $ p^2 \leq 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 $ \not k $ by $ {\not \partial}_r = \gamma^\mu\dfrac{\partial}{\partial r_\mu} $, 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 ${\rm e}^{-r^2/a}$ is moved to the left through the differential operator $ \not \partial_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\le \alpha_i<\infty $ 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,\alpha) $ arises from trace evaluation.

At this stage, an infrared cutoff is introduced to avoid any possible thresholds in the Feynman diagram:

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 parameter $ \Lambda_{M_i} $, meson-quark coupling $ g_{M_i} $, two of the four constituent quark masses $ m_{q_j} $, ($ m_u=m_d, m_c, m_s, m_b $), and universal cutoff parameter λ. If $ n_M $ mesons are considered, one has $ 2n_M+5 $ model parameters. However, the compositeness condition in Eq. (4) provides $ n_M $ constraints. These constraints are used to exclude the coupling $ g_{M_i} $ from the parameter set. The remaining $ n_M+5 $ parameters are obtained by fitting the experimental data. Naturally, we first chose experimental values for the leptonic decay constants of $ n_M $ mesons to fit the model parameters. Whenever experimental data were not available, or were available but with large errors, we decided to use LQCD results instead.

The matrix elements of the leptonic decays are described by the Feynman diagram shown in Fig. 3. The leptonic decay constants of the pseudoscalar and vector mesons are defined in the CCQM as

Figure 3.  Quark model diagram for meson leptonic decay.

where $ N_c=3 $ is the number of colors and $ O^\mu=\gamma^\mu(1-\gamma_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 2. The results of (overconstrained) least-squares fitting for the leptonic decay constants and electromagnetic decay widths are reported in Tables 1 and 2. The results of the fitting for the values of quark masses, size parameters, and infrared cutoff parameter are reported in Tables 3 and 4.

	Fit Values	Input	Ref.			Fit Values	Input	Ref.
$ f_\pi $	128.7(6.4)	$ 130.4\pm 0.2 $	[33, 34]		$ f_\omega $	198.5(9.9)	$ 198\pm 2 $	[33]
$ f_K $	156.1(7.8)	$ 156.1\pm 0.8 $	[33, 34]		$ f_\phi $	228.2(11.4)	$ 227\pm 2 $	[33]
$ f_D $	205.9(10.3)	$ 206.7\pm 8.9 $	[33, 34]		$ f_{J/\psi} $	415.0(20.8)	$ 415\pm 7 $	[33]
$ f_{D_s} $	257.5(12.9)	$ 257.5\pm 6.1 $	[33, 34]		$ f_{K^*} $	213.7(10.7)	$ 217\pm 7 $	[33]
$ f_B $	191.1(9.6)	$ 192.8\pm 9.9 $	[35]		$ f_{D^*} $	243.3(12.2)	$ 245\pm 20 $	[37]
$ f_{B_s} $	234.9(11.7)	$ 238.8\pm 9.5 $	[35]		$ f_{D_s^*} $	272.0(13.6)	$ 272\pm 26 $	[37]
$ f_{B_c} $	489.0(24.5)	$ 489\pm 5 $	[36]		$ f_{B^*} $	196.0(9.8)	$ 196\pm 44 $	[37]
$ f_\rho $	221.1(11.1)	$ 221\pm 1 $	[33]		$ f_{B_s^*} $	229.0(11.5)	$ 229\pm 46 $	[37]

Table 1.  Input values for the leptonic decay constants $ f_M $ (in MeV) and our least-squares fitting values [38].

Process	Fit Values	Input [33]		Process	Fit Values	Input [33]
$ \pi^0\to \gamma\gamma $	$ 5.06(25)\times 10^{-3} $	$ (7.7\pm 0.4)\times 10^{-3} $		$ K^{*+}\to K^+\gamma $	55.1(2.8)	$ 50\pm 5 $
$ \eta_c\to \gamma\gamma $	1.61(8)	$ 1.8\pm 0.8 $		$ K^{*0}\to K^0\gamma $	116(6)	$ 116\pm 10 $
$ \rho^+\to \pi^+\gamma $	76.0(3.8)	$ 67\pm 7 $		$ D^{*+}\to D^+\gamma $	1.22(6)	$ 1.5\pm 0.5 $
$ \omega\to \pi^0\gamma $	672(34)	$ 703\pm 25 $		$ J/\psi\to \eta_c\gamma $	1.43(7)	$ 1.58\pm 0.37 $

Table 2.  Input values for some basic electromagnetic decay widths and our least-squares fitting values (in keV) [38].

$ \Lambda_{\pi} $	$ \Lambda_K $	$ \Lambda_{D} $	$ \Lambda_{D_s} $	$ \Lambda_B $	$ \Lambda_{B_s} $	$ \Lambda_{B_c} $	$ \Lambda_{\rho} $
0.87(10)	1.04(26)	1.47(39)	1.57(42)	1.88(26)	1.95(27)	2.42(10)	0.61(7)
$ \Lambda_{\omega} $	$ \Lambda_{\phi} $	$ \Lambda_{J/\psi} $	$ \Lambda_{K^*} $	$ \Lambda_{D^*} $	$ \Lambda_{D^*_s} $	$ \Lambda_{B^*} $	$ \Lambda_{B_s^*} $
0.47(5)	0.88(10)	1.48(17)	0.72(18)	1.16(22)	1.17(22)	1.72(27)	1.71(27)

Table 3.  Meson size parameters (in GeV) [38].

$ m_{u/d} $	$ m_s $	$ m_c $	$ m_b $	λ
0.235(12)	0.424(19)	2.160(71)	5.090(21)	0.181(9)

Table 4.  Quark masses and infrared cutoff parameter (in GeV) [38].

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−42], but also baryons [43−45], tetraquarks [46, 47], and other multiquark states [48, 49] in a consistent way.

Concerning the estimation of theoretical errors in the CCQM, we used MINUIT for fitting, which is based on $ \chi^2 $ minimization. When we performed the fitting, we obtained the best-fitting values for model parameters, which are listed in Tables 3 and 4. 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 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\gamma} $ to be 10% and the error for the decay widths to be 20%.

The radiative decays $ D^*_{(s)}\to D_{(s)}\gamma $ are described by the Feynman diagrams shown in Fig. 4. The coupling between the quarks and the photon is described by the interaction Lagrangian,

Figure 4.  Feynman diagrams for radiative decays $ D^*_{(s)}\to D_{(s)}\gamma $.

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 $\omega_1=m_{c}/ (m_{c}+m_{q})$ and $ \omega_2=m_{q}/(m_{c}+m_{q}) $ with $ q=u,d,s $.

To calculate the matrix element $ \mathcal{M}(D^*\to D\gamma) $, the Gaussian form is substituted for the vertex functions in Eq. (3) and Fock-Schwinger representation of quark propagators in Eq. (11), and then the techniques described by Eqs. (12)–(19) are applied.

Next, using the transversality conditions $ \epsilon_\mu^\gamma(q) q^\mu = 0 $ and $ \epsilon_\nu^{D^*}(p) p^\nu = 0 $, the matrix element can be expressed in the form

where $ g_{D^*D\gamma}=e_c\,I_c(m^2_{D^*},m^2_D)+e_q\,I_q(m^2_{D^*},m^2_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

The expression for $ I_{q}(m^2_{D^*},m^2_D) $ can be obtained by simply exchanging $ m_c \leftrightarrow m_q $ and $ \omega_1 \leftrightarrow \omega_2 $, i.e. $ I_{q}(m^2_{D^*},m^2_D,m_q, m_c, \omega_1,\omega_2)=I_{c}(m^2_{D^*},m^2_D,m_c, m_q, \omega_2,\omega_1) $.

Finally, the decay width for $ D^*\to D\gamma $ is given by

Our results for the radiative decay constants are listed in Table 5 together with results obtained in other theoretical approaches. It is worth mentioning that the sign difference in the predictions in Table 5 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^+\gamma}=-0.47(7)\,\; \rm{GeV}^{-1}$ was extracted from the total width of $ D^{*+} $ and measured branching fraction $ \mathcal{B}(D^{*+}\to D^+\gamma) $ [25]. Meanwhile, the total width of $ D^{*0} $ is still unknown experimentally. However, it can be calculated using the isospin symmetry related to the strong couplings $ g_{D^{*+}D\pi} $ and $ g_{D^{*0}D\pi} $ [23, 51], which yields $ \Gamma(D^{*0})=56.5(14.0)\, \rm{keV} $ [25]. Then, the value $ g_{D^{*0}D^0\gamma}= 1.77(16)\, \rm{GeV}^{-1} $ was obtained with the help of the measured $ \mathcal{B}(D^{*0}\to D^{0}\gamma) $ [25]. Note that our predictions $ g_{D^{*+}D^+\gamma}=-0.45(9)\, \rm{GeV}^{-1} $ and $g_{D^{*0}D^0\gamma}=1.72(34)\, \rm{GeV}^{-1}$ are in full agreement with the current experimental data. Note also that the decay width of the channel $ D^{*+}\to D^+\gamma $ was included in the model parameter fitting (see Ref. [38] and Table 2). However, channels $ D^{*0}\to D^0\gamma $ and $ D^{*+}_s\to D_s^+\gamma $ were not included in the fitting.

Ref.	$ g_{D^{*+}D^+\gamma} $	$ g_{D^{*0}D^0\gamma} $	$ g_{D^{*+}_sD^+_s\gamma} $
This work	$ -0.45(9) $	$ 1.72(34) $	$ -0.29(6) $
LCSR (NLO) [25]	$ 0.40^{+0.12}_{-0.13} $	$ -2.11^{+0.35}_{-0.34} $	$ 0.60^{+0.19}_{-0.18} $
LCSR (NLO)[10]	$ -0.15^{+0.11}_{-0.10} $	$ 1.48^{+0.29}_{-0.27} $	$ -0.079^{+0.086}_{-0.078} $
HQET$ + $CQM [22]	$ -0.38^{+0.05}_{-0.04} $	$ 1.91^{+0.09}_{-0.09} $	–
LQCD [23]$ ^\dagger $, [24]$ ^* $, [50]$^{**}$	$ -0.2(3)^\dagger $	$ 2.0(6)^\dagger $	$0.11(2)^*, 0.04(1)^{**}$
Bag Model [16]	0.5	1.1	–
RQM [18]	$ -0.30 $	1.85	–
RQM [19]	$ -0.44(6) $	$ 2.15(11) $	$ -0.19(3) $
LCSR [11]	$ -0.50(12) $	$ 1.52(25) $	–
QCDSR [7]	$ -0.19^{+0.03}_{-0.02} $	$ 0.62^{+0.03}_{-0.03} $	$ -0.20^{+0.03}_{-0.03} $
HQET$ + $VMD [52]	$ -0.29^{+0.19}_{-0.11} $	$ 1.60^{+0.35}_{-0.45} $	$ -0.19^{+0.19}_{-0.08} $
HHχPT [53]	$ -0.27(5) $	$ 2.19(11) $	$ 0.041(56) $
Experiment [54]	$ -0.47(7) $	$ 1.77(16) $	–

Table 5.  Radiative decay constants in the CCQM and other approaches (all in $ \rm{GeV}^{-1} $).

In the case of $ g_{D^{*+}_sD^+_s\gamma} $, the theoretical predictions vary significantly with respect to the case of $ g_{D^{*+}D^+\gamma} $ and $ g_{D^{*0}D^0\gamma} $. The central value of our result, $|g_{D^{*+}_sD^+_s\gamma}| = 0.29(6)~ \rm{GeV}^{-1}$, is larger than all other theoretical predictions, except for a recent value, $|g_{D^{*+}_sD^+_s\gamma}|= 0.60(19)\,\; \rm{GeV}^{-1}$, obtained by using LCSR(NLO) [25]. Our result deviates from others within $ 2\sigma $ owing to large uncertainties in the predictions. In Ref. [25], a discussion on D-spin symmetry for the decays $ D^{*+}\to D^+\gamma $ was also provided. Here, the D-spin symmetry is understood as the symmetry under the exchange of down and strange quarks $ d\leftrightarrow s $. As pointed out in Ref. [25], this is a good approximate symmetry in QED. Given that the decay $ D_s^{*+}\to D_s^+\gamma $ can be obtained from the decay $ D^{*+}\to D^+\gamma $ by exchanging $ d\leftrightarrow 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^+\gamma} $ and $ g_{D^{*+}_sD^+_s\gamma} $ is approximately $ 20{\text{%}} $–$ 30{\text{%}} $. Our results show a deviation of approximately $ 35{\text{%}} $, which is reasonable. Note also from Table 5 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 $ \mathcal{O}(\alpha_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].

Table 6 compares the predictions for the width of the$ D^*_{(s)} $ radiative decays. Available theoretical predictions for $ \Gamma(D^{*+}\to D^+\gamma) $ vary in the range $ (0.2 - 1.7) $ keV. Our result $ \Gamma(D^{*+}\to D^+\gamma)=1.21(48) $ agrees well with experimental data, suggesting the width of the decay $ D^{*+}\to D^+\gamma $ to be at the higher end of the predicted range. Note that the experimental value $ \Gamma(D^{*+}\to D^+\gamma)= 1.33(36)~ \rm{keV} $ was obtained from the world average of the total width $ \Gamma(D^{*+})=83.4(1.8)\, \rm{keV} $ and branching fraction $ \mathcal{B}(D^{*+}\to D^+\gamma)=0.016(4) $ reported by Particle Data Group (PDG) [54]. We stress again that the decay $ D^{*+}\to D^+\gamma $ was included in our model parameter fitting (see Table 2). However, this fitting was performed in 2011 [38], and the input value $ \Gamma(D^{*+}\to D^+\gamma)=1.5(5)\, \rm{keV} $ was taken from PDG at that time [33]. The slight difference for the central values of $ \Gamma(D^{*+}\to D^+\gamma) $ given in Tables 6 and 2, which are 1.21 and 1.22, respectively, is because we used the most updated values for meson masses in this study.

Ref.	$ \Gamma(D^{*+}\to D^+\gamma) $	$ \Gamma(D^{*0}\to D^0\gamma) $	$ \Gamma(D^{*+}_s\to D^+_s\gamma) $
This work	$ 1.21(48) $	$ 18.4(7.4) $	$ 0.55(22) $
LCSR (NLO) [25]	$ 0.96^{+0.58}_{-0.62} $	$ 27.83^{+9.23}_{-9.50} $	$ 2.36^{+1.49}_{-1.41} $
LQCD [23]$ ^\dagger $, [24]$ ^* $	$ 0.8(7)^\dagger $	$ 27(14)^\dagger $	$ 0.066(26)^* $
HQET$ + $CQM [22]	$ 0.9\pm 0.3 $	$ 22.7\pm 2.2 $	–
Bag Model [16]	1.73	7.18	–
NJL Model [17]	0.7	19.4	0.09
RQM [14]	0.46	23.05	0.38
RQM [18]	0.56	21.69	–
RQM [19]	(0.94–1.42)	(26.0–32.0)	(0.2–0.3)
RQM [20]	1.04	11.5	0.19
LFQM [21]	$ 0.90(2) $	$ 20.0(3) $	$ 0.18(1) $
QCDSR [9]	$ 0.23(10) $	$ 12.9(2.0) $	$ 0.13(5) $
LCSR [11]	1.50	14.40	–
QCDSR [7]	$ 0.22(6) $	$ 2.43(21) $	$ 0.25(8) $
HQET$ + $VMD [52]	$ 0.51(57) $	$ 16.0(10.8) $	$ 0.24(24) $
Experiment [54]	$ 1.33(36) $	$ 19.9(5.5) $	–

Table 6.  Decay widths $ \Gamma(D^*_{(s)}\to D_{(s)}\gamma) $ in the CCQM and other approaches (all in $ \rm{keV} $).

The experimental value $ \Gamma(D^{*0}\to D^0\gamma)=19.9(5.5)\, \rm{keV} $ in Table 6 was calculated by using the estimation $ \Gamma(D^{*0})=56.5(14.0)\, \rm{keV} $ mentioned above and the branching fraction $ \mathcal{B}(D^{*0}\to D^0\gamma)=0.353(9) $ fitted by PDG [54]. Our prediction $ \Gamma(D^{*0}\to D^0\gamma)=18.4(7.4)\, \rm{keV} $ also agrees well with this experimental value.

At this point, the ratio of decay widths $ \Gamma(D^{*+}\to D^+\gamma)/\Gamma(D^{*0}\to D^0\gamma) $can be checked. 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\to\infty $, the integral $ I_c $ is suppressed as $ 1/m_c $ and the width ratio yields $ 0.25 $. However, from Table 6, one obtains $ \Gamma(D^{*+}\to D^+\gamma)/ \Gamma(D^{*0}\to D^0\gamma)= 0.066 $. Therefore, the use of HQL is not suitable for the decays $ D^{*}\to D\gamma $. This is opposed to the case of the decays $ B^{*}\to B\gamma $, 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\to D_s\gamma $, there is still much speculation about its width. Predictions for $ \Gamma(D^*_s\to D_s\gamma) $ range from 0.07 keV to 2.4 keV. Our result $\Gamma(D^*_s\to D_s\gamma)=0.55(22)$ keV deviates from the recent prediction based on LCSR at NLO $ \Gamma(D^*_s\to D_s\gamma)= 2.36^{+1.49}_{-1.41} $ keV [25] and the first LQCD result $ \Gamma(D^*_s\to D_s\gamma)=0.066(26) $ keV [24] by approximately $ 2\sigma $. Note also that the central value of our result $ \Gamma(D^*_s\to D_s\gamma)=0.55(22) $ keV is larger than that of most predictions; in particular, it exceeds the LQCD central value $ \Gamma(D^*_s\to D_s\gamma)=0.066 $ keV [24] by roughly an order of magnitude. By contrast, the recent prediction based on LCSR at NLO $\Gamma(D^*_s\to D_s\gamma)= 2.36^{+1.49}_{-1.41}$ keV [25] has a central value approximately four times as large as our value. Other theoretical studies seem to agree on the central value $ \Gamma(D^*_s\to D_s\gamma)\approx 0.2 $ keV.

It is instructive to estimate the width of the decay $ D^*_s\to D_s\gamma $ using the recent measurement of the leptonic decay $ D_s^{*+}\to e^+\nu_e $ performed by the BESIII collaboration [56]. This is the first measurement of the decay branching fraction; it reads $\mathcal{B}(D_s^{*+}\to e^+\nu_e) = (2.1^{+1.2}_{-0.9_{ \rm{stat.}}}\pm 0.2_{ \rm{syst.}})\times 10^{-5}$. By taking the ratio of branching fractions $ \mathcal{B}(D_s^{*+}\to e^+\nu_e)/\mathcal{B}(D_s^{+}\to \mu^+\nu_\mu) $, one obtains

Using the world average $ \mathcal{B}(D_s^{+}\to \mu^+\nu_\mu) $ and average value $ \dfrac{f_{D_s^{*+}}}{f_{D^+_s}}=1.12(1) $ given by LQCD calculations, the BESIII collaboration found the total width of $ D_s^{*+} $ to be $ \Gamma^{ \rm{total}}_{D_s^{*+}}=121.9^{+69.6}_{-52.2}\pm 11.8\, \rm{eV} $ [56], in agreement with the LQCD prediction $ 70(28)\, \rm{eV} $ [24] within $ 1\sigma $. The world average value $ \mathcal{B}(D^{*+}_s\to D^+_s\gamma)=0.381(29) $ [54] can now be used to obtain $\Gamma(D^{*+}_s\to D^+_s\gamma)=46(30)\,\; \rm{eV}$, which is much smaller than most theoretical predictions listed in Table 6, including ours. It should be noted that this estimation is based on the first measured branching fraction $ \mathcal{B}(D_s^{*+}\to e^+\nu_e)=(2.1^{+1.2}_{-0.9_{ \rm{stat.}}}\pm 0.2_{ \rm{syst.}})\times 10^{-5} $, which still suffers from large uncertainties. A more precise measurement of $ \mathcal{B}(D_s^{*+}\to e^+\nu_e) $, likely to be obtained in the near future, will help clarify the value of the decay width $ \Gamma(D^{*+}_s\to D^+_s\gamma) $.

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 $ \pm 20\% $ with respect to its fitted value $ m_s = 0.424 $ GeV while keeping other parameters unchanged. The dependence of the width $ \Gamma(D^*_s\to D_s\gamma) $ 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.

Figure 5.  (color online) $ \Gamma(D^*_s\to D_s\gamma) $ vs. strange-quark mass.

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 cancelation between the contributions from the charm and down quarks in $ g_{D^{*+}D^+\gamma} $, as well as from the charm and strange quarks in $ g_{D_s^{*+}D_s\gamma} $. Meanwhile, for $ g_{D^{*0}D^0\gamma} $, the contributions from the charm and up quarks add up (see Table 5). The strong effects in the decays are captured in the quantities $ I_{c(q)} $. Table 7 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 $.

Channel	Quark q	$ I_q/I_c $	$ e_q I_q/e_c I_c $	$ g_{D^{*}D\gamma} $ (GeV$ ^{-1} $)
$ D^{*+}\to D^+\gamma $	d	3.74	$ -1.87 $	$ -0.45(9) $
$ D^{*0}\to D^0\gamma $	u	3.74	$ +3.74 $	$ +1.72(34) $
$ D^{*+}_s\to D^+_s\gamma $	s	2.97	$ -1.49 $	$ -0.29(6) $

Table 7.  Comparison of quark contributions to radiative decay constant.

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^{+}\gamma} $, $ g_{D^{*0}D^{0}\gamma} $, and $ g_{D^{*+}_s D^{+}_s\gamma} $, as well as the decay widths are reported. A full agreement between these predictions and experimental data was found for the decays $ D^{*+}\to D^{+}\gamma $ and $ D^{*0}\to D^{0}\gamma $. In the case of $ D^{*+}_s\to D^{+}_s\gamma $, our result for the decay width $ 0.55(22)\, \rm{keV} $ is smaller than the recent LCSR prediction $ 2.36^{+1.49}_{-1.41}\, \rm{keV} $ by more than $ 1\sigma $ but exceeds the first LQCD prediction $ 0.066(26)\, \rm{keV} $ by more than $ 2\sigma $. Using the branching fraction of the purely leptonic decay $ D_s^{*+}\to e^+\nu_e $ recently measured by the BESIII collaboration, we estimated $ \Gamma(D^{*+}_s\to D^{+}_s\gamma) $ to be approximately $ 0.05(3)\, \rm{keV} $. Meanwhile, theoretical predictions for this width range from $ 0.07\, \rm{keV} $ to $ 2.4\, \rm{keV} $, and most of them are around $ 0.2\, \rm{keV} $. In particular, this estimated value is smaller than our prediction by more that $ 2\sigma $. Therefore, the decay $ D^{*+}_s\to D^{+}_s\gamma $ is an interesting case that requires further theoretical and experimental studies to give a final answer.

CTT deeply thanks Dr. Anh-Vu Phan-Gia and Dr. Hai-Cat Tran for their help and fruitful discussion. This paper is dedicated to the memory of our dear friend and colleague Jürgen Körner.

Let us consider the simplest case of a zero-spin bound state whose constituents are fermions. The Yukawa-type "bare" Lagrangian describing the coupling of the boson field $ \phi_0 $ to the fermion field q has the form

The vacuum generating functional of the Yukawa theory is given by

Integrating out quark fields and dismissing all irrelevant normalization factors, we have

where $S_q(x-y)=\int\dfrac{{\rm d}^4k}{(2\pi)^4i}\dfrac{{\rm e}^{-{\rm i}k(x-y)}}{m_q-\not k}$ is the quark Green function (propagator).

The bi-linear terms in boson field read

where $\Pi_{S=0}(x_1-x_2) = {\rm i} \langle T[(\bar{q}\Gamma q)_x(\bar{q}\Gamma q)_y]\rangle_0 = -{\rm i}\, \rm{tr}[\Gamma S_q (x_1-x_2)\Gamma S_q(x_2-x_1)]$ is the mass function of the spin-0 boson.

Expanding the Fourier transform of the mass function at the physical boson mass up to second order

we obtain

The wave function, boson mass, and Yukawa coupling are renormalized as follows:

Note that the wave function renormalization constant Z is related to the renormalized coupling constant $ g_r $ via

The renormalized generating functional of the Yukawa theory is finally expressed as

We dismiss the linear boson term because it is absent for pseudoscalar mesons; for scalar mesons, it can be removed by a shift of the field.

Next, we move on to the Fermi theory for which the generating functional is expressed as

where the Fermi-type Lagrangian $ \mathcal{L}_F(x) $ is given by

For the exponential of the four-fermion interaction, the Gaussian functional representation can be used as follows:

We then have

Integrating out quark fields and dismissing all irrelevant normalization factors, we obtain

All bi-linear terms in the boson field can be collected, and the renormalized mass function is introduced in a similar manner to that of the Yukawa case as follows:

If the condition $ G\tilde{\Pi}_{S=0}(m^2)=1 $ is required and the boson field is rescaled as $ \phi\to\dfrac{\phi}{\sqrt{\tilde{\Pi}^\prime_{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 reads

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

Taking into account the relation (A8), this condition is equivalent to

Therefore, the vanishing of the wave function renormalization constant in the Yukawa theory can be used as the condition for the bare (unrenormalized) field $ \phi_0 = Z^{1/2}\phi_r $ to vanish for a composite boson.