Probing the ATOMKI X17 vector boson using Dalitz decays V → Pe⁺e⁻

In 2015, the ATOMKI collaboration reported an anomaly in the M1 transition of the 18.15 MeV isoscalar excited state ($ J^\pi = 1^+, T = 0 $) to the ground state ($ J^\pi = 0^+, T = 0 $) in $ ^8 {{\rm{Be}}} $ [1, 2]. The excited state of $ ^8 {{\rm{Be}}} $ was created by firing proton beam at thin strip foils of $ ^7 {{\rm{Li}}} $. The first then decayed to the ground state and produced pairs of electrons and positrons. The $ e^+e^- $ emission is referred to as the internal pair creation (IPC). The angular correlation of the IPC process (IPCC) was measured at various opening angles θ in the laboratory rest frame. Quantum electrodynamics (QED) predicts the IPCC to drop rapidly and monotonically with the increase in the opening angle [3, 4]. However, a deviation was observed at the proton beam energy $ E_p = 1.10\; {\rm{MeV}} $ and at the angle $ \theta \approx 140^\circ $ with a $ 6.8\sigma $ significance. The deviation appeared as a “bump” in the angular correlation spectrum, suggesting an additional contribution to the $ e^+e^- $ production from a hypothetical bosonic particle of mass $ \approx 17\; {\rm{MeV}} $, namely X17, apart from the main photon contribution.

After the first observation of the anomaly in the $ ^7 {{\rm{Li}}}( {\rm{p}},e^+e^-)^8 {{\rm{Be}}} $ nuclear reaction [1], the ATOMKI collaboration further studied the reactions $ ^3 {\rm{H}}( {\rm{p}},e^+e^-)^4 {\rm{He}} $ [5] and $ ^{11} {{\rm{B}}}( {\rm{p}},e^+e^-)^{12} {\rm{C}} $ [6] and reported the observation of similar anomalies. The $ e^+e^- $ excesses in the decays of excited $ ^4 {\rm{He}} $ and $ ^{12} {\rm{C}} $ were observed at different opening angles compared to the case of excited $ ^8 {{\rm{Be}}} $, but the inferred mass of the unknown particle remained consistent with the $ ^8 {{\rm{Be}}} $ result. A brief summary of the ATOMKI's results is given in Table 1. Independent experimental attempts from other collaborations to test the ATOMKI anomaly have been done. Based on data samples with total statistics corresponding to $ 8.4\times 10^{10} $ electrons collected in 2017 and 2018, the NA64 collaboration at CERN searched for the X17 produced in the bremsstrahlung reaction $ e^-Z\to e^-ZX $. No evidence was found and an improved limit on the $ X-e^- $ coupling was set to be $ 1.2\times 10^{-4}\lesssim \varepsilon_e \lesssim 6.8\times 10^{-4} $ [7−9]. In 2023, the MEG II collaboration at PSI reported no significant signal above the expected background in the $ e^+e^- $ angular correlation spectrum that would indicate the presence of the X17 boson [10]. However, their results were found to be compatible with the ATOMKI's within $ 1.5\sigma $. Meanwhile, a group at VNU University of Science reported the observation of the anomaly in the $ ^7 {{\rm{Li}}}( {\rm{p}},e^+e^-)^8 {{\rm{Be}}} $ nuclear reaction at proton beam energy 1.225 MeV and opening angle around $ 135^\circ $. The mass was found to be $ m_X = 16.66\pm 0.47( {{\rm{stat}}})\pm 0.35 ( {{\rm{sys}}}) $ MeV with a confidence above $ 4\sigma $ [11]. Recently, the PADME collaboration has performed a dedicated search for the X17 boson via $ e^+e^- $ annihilation and found an excess of events over the predicted background expectation with a local significance of $ 2.5\sigma $ for an X17 mass of 16.90 MeV [12, 13].

Transition	Opening angle θ	$ E_p $ (MeV)	Best fit $ m_X $ (MeV)	Confidence	Ref.
$ ^7 {{\rm{Li}}}( {\rm{p}},e^+e^-)^8 {{\rm{Be}}} $	$ \approx 140^\circ $	1.10	$ 16.70\pm 0.35 ( {{\rm{stat}}})\pm 0.50 ( {{\rm{sys}}}) $	$ 6.8\sigma $	[1]
$ ^3 {\rm{H}}( {\rm{p}},e^+e^-)^4 {\rm{He}} $	$ \approx 115^\circ $	0.51	$ 17.01 \pm 0.12 ( {{\rm{stat}}})\pm 0.21 ( {{\rm{sys}}}) $	$ 7.3\sigma $	[5]
		0.61	$ 16.88 \pm 0.16 ( {{\rm{stat}}})\pm 0.21 ( {{\rm{sys}}}) $	$ 6.6\sigma $
		0.90	$ 16.68 \pm 0.30 ( {{\rm{stat}}})\pm 0.21 ( {{\rm{sys}}}) $	$ 8.9\sigma $
$ ^{11} {{\rm{B}}}( {\rm{p}},e^+e^-)^{12} {\rm{C}} $	$ \approx 155^\circ - 160^\circ $	1.50	$ 16.81 \pm 0.15 ( {{\rm{stat}}})\pm 0.20 ( {{\rm{sys}}}) $	$ 3.0\sigma $	[6]
		1.70	$ 16.93 \pm 0.08 ( {{\rm{stat}}})\pm 0.20 ( {{\rm{sys}}}) $	$ 7.0\sigma $
		1.88	$ 17.13 \pm 0.10 ( {{\rm{stat}}})\pm 0.20 ( {{\rm{sys}}}) $	$ 8.0\sigma $
		2.10	$ 17.06 \pm 0.10 ( {{\rm{stat}}})\pm 0.20 ( {{\rm{sys}}}) $	$ 3.0\sigma $

Table 1.  ATOMKI anomaly summary.

The ATOMKI anomaly has attracted a great deal of attention in the particle physics community and has led to a large number of theoretical studies searching for possible explanations. One promising solution to the anomaly is the “protophobic” vector boson explanation proposed by Feng et al. [14, 15] (see, however, Ref. [16]). They also pointed out that the X17 boson could not be a scalar due to angular momentum conservation in the decay $ ^8 {{\rm{Be}}}^*\to\, ^8 {{\rm{Be}}}+e^+e^- $. Kozaczuk et al. considered the possibility of the X17 being a light axial vector boson and found that such a scenario is consistent with experimental data [17] (see also Ref. [18]). Ellwanger and Moretti studied the pseudoscalar boson explanation in Ref. [19]. In Ref. [20], Kirpichnikov, Lyubovitskij, and Zhevlakov performed a global analysis of possible impacts of hidden sub-GeV bosons with different spin-parity quantum numbers (including the X17 state) on the ATOMKI anomaly, $ (g-2)_\mu $ anomaly, proton-charge-radius puzzle, and electric dipole moments of fermions. A large variety of models following the above mentioned directions has been proposed and studied [21−25]. A detailed review of the X17 anomaly can be found in Ref. [26].

If the X17 exists, it would affect similar processes in particle physics, including the Dalitz decays of vector mesons. Recently, the BESIII collaboration measured the Dalitz decay $ D^{*0}\to D^0e^+e^- $ for the first time [27] and reported a $ 3.5\sigma $ excess over the theoretical prediction based on the vector meson dominance (VMD) model [28, 29]. This excess may be another signal of the X17. The idea of using the Dalitz decays to probe for new physics beyond the Standard Model (SM) is not new. For instance, the Dalitz decays $ J/\psi\to P\ell^+\ell^- $ have been proposed [30] as a probe for dark photon searches at the BESIII experiment, where a large $ J/\psi $ data sample has been accumulated [31]. In the light of the ATOMKI anomaly, Castro and Quintero proposed tests of the anomaly using the Dalitz decays $ D^*_{(s)}\to D_{(s)}e^+e^- $ and $ B^*_{(s)}\to B_{(s)}e^+e^- $ [32], assuming that X17 is a vector boson. Ban et al. [33] studied various strategies for searching the X17 boson in the Dalitz decay $ J/\psi\to \eta_c\ell^+\ell^- $ at BESIII and Belle II. Recently, Lee et al. [34] performed a fit for the couplings between the X17 and quarks using experimental data from several Dalitz channels, including $ \psi(2S)\to \eta_c\ell^+\ell^- $, $ \phi\to \eta\ell^+\ell^- $, $ D_s^*\to D_s\ell^+\ell^- $, and $ D^{*0}\to D^0\ell^+\ell^- $. Note that in Refs. [32−34] the authors relied on the VMD model and the heavy quark effective theory (HQET) to calculate the hadronic form factors. Also, there were some disagreements on the prediction for $ D_s^*\to D_s\ell^+\ell^- $ between Ref. [32] and Ref. [34].

In this study, we aim at providing independent predictions for the Dalitz decays $ D^*_{(s)}\to D_{(s)}e^+e^- $, $ B^*_{(s)}\to B_{(s)}e^+e^- $, and $ J/\psi\to \eta_ce^+e^- $ in the SM, as well as in the presence of the X17. We follow the authors of Refs. [32, 34] and assume that X17 is a vector boson. We use several favored combinations of the coupling parameters between the X17 and quarks to calculate the X17's effects on the widths of these decays and compare our results with those given in Refs. [32, 34]. The hadronic form factors required for the calculation of the Dalitz decays are obtained in the Covariant Confined Quark Model (CCQM) which has been developed by our group. These form factors are calculated directly in our model for the entire range of momentum transfer without any extrapolations, heavy-quark limit or VMD assumptions. Therefore, our results would shed more light on the understanding of the X17 effects and help probe for the X17 boson in the Dalitz decays of vector mesons.

The rest of the paper is organized as follows. In Section II we briefly introduce the Covariant Confined Quark Model as a tool for hadronic calculation. Section III is devoted for the description of the Dalitz decays in terms of the Feynman diagrams, invariant matrix elements, and form factors. Here we consider both the contributions from the photon (within the SM) and from the X17 vector boson (beyond the SM). Numerical results are presented in Section IV. A short description of the form-factor evaluation within the CCQM is also provided. Finally, a brief summary is given in Section V.

The covariant confined quark model is a framework based on quantum field theory that describes hadronic bound states as quantum fields interacting with their quark constituents. This interaction is defined by a Lagrangian. For a meson M, the Lagrangian is written as

where $ g_M $ is the meson-quark coupling constant, $ \Gamma_M $ is the Dirac matrix with quantum numbers specific to the meson M, and $ F_M(x;x_1,x_2) $ is the vertex function that effectively describes the meson's size and the relative positioning of the quarks and the meson as a whole. The vertex function has the form

where $ \omega_{i}=m_{q_i}/(m_{q_1}+m_{q_2}) $ is the mass fraction of the quarks. The function $ \Phi_M[(x_1-x_2)^2] $ is assumed to be Gaussian for simplicity, and takes the following form in momentum space

Here, $ \Lambda_M $ is a free parameter of the model, referred to as the size parameter of the meson M. It is important to note that the specific mathematical form chosen for $ \widetilde{\Phi}_M(-p^2) $ is not critical, as long as it decreases rapidly enough at high momentum values (in the Euclidean space) to ensure the ultraviolet finite of Feynman diagrams.

The normalization of particle-quark vertices is provided by the compositeness condition [35, 36]

where $ Z_M $ is the meson's wave function renormalization constant and $ \Pi'_M $ is the derivative of the meson's mass function. The meson mass function $ \Pi_M(p) $ is calculated using a one-loop Feynman diagram representing the meson's self-energy (see Fig. 1). The specific forms for pseudoscalar $ (\Pi_P) $ and vector $ (\Pi_V) $ mesons read

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

where the free quark propagator has been used

The CCQM has several adjustable parameters: the constituent quark masses $ m_q $, the size parameters for hadrons $ \Lambda_H $, and a universal cutoff parameter λ that accounts for the confinement of quarks within hadrons at low energies. These parameters are determined by fitting the model's predictions to experimental data. Once these parameters are fixed, the CCQM can be used to perform hadronic calculations across the entire physical range of momentum transfer without needing to extrapolate. A key advantage of the model is its ability to consistently describe not only mesons [37−42] but also baryons [43−45], tetraquarks [46, 47], and other multiquark states [48, 49]. The values of the model parameters used in this work are provided in Tables 2 and 3. Other physical quantities and model-independent parameters are taken from the Particle Data Group [50].

$ \Lambda_{D} $	$ \Lambda_{D_s} $	$ \Lambda_B $	$ \Lambda_{B_s} $	$ \Lambda_{\eta_c} $	$ \Lambda_{D^*} $	$ \Lambda_{D_s^*} $	$ \Lambda_{B^*} $	$ \Lambda_{B_s^*} $	$ \Lambda_{J/\psi} $
1.600	1.748	1.963	2.050	3.777	1.529	1.558	1.805	1.794	1.738

Table 2.  Meson size parameters (in GeV).

$ m_{u/d} $	$ m_s $	$ m_c $	$ m_b $	λ
0.241	0.428	1.672	5.046	0.181

Table 3.  Quark masses and infrared cutoff parameter (in GeV).

To estimate theoretical uncertainties within the CCQM, we note that we employed the MINUIT program, which finds the best parameter values by minimizing the $ \chi^2 $ function when fitting the model to data. Instead of a complex error propagation analysis, we adopt a simpler, though less precise, method for estimating errors in our physical predictions. We proceeded with calculations using only these best-fit parameter values. Upon comparison, we noted that the model's fitted values typically deviated from experimental data by about 5%-10%. Since our calculations of hadronic quantities are analogous to those for leptonic and electromagnetic decay constants used in the fit, we estimate the error in the CCQM predictions to be approximately 10%.

Within the SM, the Dalitz decay has been calculated previously in many theoretical studies (see, e.g., Refs. [51−53]). We therefore do not repeat the calculation and just write down the well known expression for the decay width distribution as

The physical range for the momentum transfer is given by $ 4m^2_\ell \leq q^2 \leq (m_V-m_P)^2 $.

In the CCQM, the transition form factor $ F^\gamma_{VP}(q^2) $ is calculated based on the Feynman diagrams in Fig. 2, which describe direct photon emission from the constituent quarks. Note that within the CCQM, resonance diagrams corresponding to photon emission via intermediate $ V^\prime\to\gamma^* $ transitions are also allowed (see, e.g., Ref. [51]). However, since the mass splittings $ \Delta m(VP)\equiv m_V-m_P $ for the decays studied in this paper are too small ($ \Delta m(D^*D)\approx 140\; {\rm{MeV}} $, $ \Delta m(B^*B)\approx 50\; {\rm{MeV}} $, and $ \Delta m(J/\psi \eta_c)\approx 110\; {\rm{MeV}} $) compared to masses of the vector resonances $ \rho(770) $, $ \omega(782) $, $ \phi(1020) $, and $ J/\psi(3096) $, the contributions from resonance diagrams are negligibly small, and are well below the estimated inherent error of the model predictions ($ \sim 10 $%). We therefore only consider the triangle diagrams in this study. Diagrams describing the contribution of the X17 boson to the Dalitz decays are obtained from those in Fig. 2 by replacing the photon by the X17 boson.

Figure 2.  Feynman diagrams for Dalitz decays $ V\to P\ell^+\ell^- $.

The quark-photon coupling is described by the interaction Lagrangian

where $ e_Q $ and $ e_q $ are the charges of heavy and light quarks in units of e. Similarly, for the interaction between the quarks and the X17 boson one has

where $ \varepsilon_{Q,q} $ are the coupling constants between the X17 boson and quarks of various flavors.

The invariant matrix elements of the Dalitz decay $ V(p,\epsilon_V)\to P(p^\prime)e^+(q_+)e^-(q_-) $ are given by [34]

where $ e\varepsilon_e $ is the coupling constant between the X17 boson and electron (positron), and $ T_\mu^{\gamma,X} $ are the hadronic amplitudes of the transitions $ V\to P\gamma^* $ and $ V\to PX^* $, respectively. The hadronic transition amplitudes can be parametrized in terms of the transition form factors $ g_{VP\gamma}(q^2) $ and $ g_{VPX}(q^2) $ as follows:

The differential decay rates induced by photon and X17 boson are given by

where $ \lambda(x,y,z) \equiv x^2+y^2+z^2-2(xy+yz+zx) $ is the Källén function.

By integrating the Eqs. (13) and (14) one obtains the decay width of the Dalitz decay. Note that the widths of the decays considered in this paper have not been measured directly. Instead, the CLEO [54] and BESIII [27] collaborations measured the ratios of the Dalitz decay with respective to the corresponding radiative decay, namely,

where the radiative decay width is given by [55]

Note that the interference of the amplitudes in Eq. (11) is negligible since the X17 boson's total width is very narrow [32]. Therefore, one can simply write

The differential decay rate for the case of the X17 mediator can be simplified using the narrow width approximation [34]

Additionally, by assuming that the X17 vector boson decays dominantly into $ e^+e^- $, one has

By substituting Eqs. (18) and (19) into Eq. (14), the coupling constant $ e\varepsilon_e $ cancels out and is therefore not required in the calculation.

The form factors $ g_{VP\gamma}(q^2) $ and $ g_{VPX}(q^2) $ defined in Eq. (12) are calculated within the CCQM by evaluating the hadronic amplitudes $ T_\mu^\gamma $ and $ T_\mu^X $. We show below the calculation steps for the case of $ g_{VP\gamma}(q^2) $. The calculation of $ g_{VPX}(q^2) $ is performed in a similar manner by exchanging $ e_{Q,q} \leftrightarrow \varepsilon_{Q,q} $. In the CCQM, the hadronic amplitude $ T_\mu^\gamma $ is written as

The ratios of quark masses now read $ \omega_1=m_{Q}/(m_{Q}+m_{q}) $ and $ \omega_2=m_{q}/(m_{Q}+m_{q}) $ with $ Q=b,c $, and $ q=u,d,s $.

Next, we substitute the Gaussian form for the vertex functions in Eq. (3). Note that Feynman diagrams are calculated in the Euclidean region where $ p^2=-p^2_E $. The vertex functions fall off in the Euclidean region, and therefore guarantee ultraviolet convergence. The loop integration is performed by using the Fock-Schwinger representation for the quark propagator

The integrals over the Fock-Schwinger parameters $ 0\le \alpha_i<\infty $ are treated by introducing an additional integration which converts the set of these parameters into a simplex as follows

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.

Finally, the form factor $ g_{VP\gamma}(q^2) $ can be written as

where $ I_{Q(q)}(m^2_{V},m^2_P,q^2) $ are twofold integrals which are calculated numerically. The expression for $ I_{Q}(m^2_{V},m^2_P,q^2) $ reads

The expression for $ I_{q}(m^2_{V},m^2_P,q^2) $ can be obtained by simply exchanging $ m_Q \leftrightarrow m_q $ and $ \omega_1 \leftrightarrow \omega_2 $, i.e. $ I_{q}(m^2_{V},m^2_P,m_q, m_Q, \omega_1,\omega_2)=I_{Q}(m^2_{V},m^2_P,m_Q, m_q, \omega_2,\omega_1) $.

For convenience in the calculation, we choose a double-pole parametrization to interpolate the calculated form factors as follows

This interpolation form well represents the calculated values of all form factors in this study. The parameters of the double-pole interpolation are listed in Table 4.

Transition	$ g_{VP\gamma}(q^2) $				$ g_{VPX}(q^2) $
	$ g_{VP\gamma}(0) $	a	b		$ g_{VPX}(0) $	a	b
$ D^{*+}\to D^+e^+e^- $	$ -0.32 $	4.02	9.33		$ -5.7\times 10^{-3} $	2.03	1.07
$ D^{*0}\to D^0e^+e^- $	1.74	1.55	0.13		$ 1.3\times 10^{-3} $	1.55	0.13
$ D^{*+}_s\to D^+_se^+e^- $	$ -0.18 $	3.90	9.63		$ -4.4\times 10^{-3} $	1.45	0.40
$ B^{*+}\to B^+e^+e^- $	1.27	2.25	1.14		$ 4.5\times 10^{-4} $	4.80	14.49
$ B^{*0}\to B^0e^+e^- $	$ -0.73 $	1.96	0.43		$ -6.4\times 10^{-3} $	1.96	0.43
$ B^{*0}_s\to B^0_se^+e^- $	$ -0.58 $	1.29	$ -0.05 $		$ -5.1\times 10^{-3} $	1.29	$ -0.05 $
$ J/\psi\to \eta_c e^+e^- $	0.64	0.14	0.003		$ 4.8\times 10^{-4} $	0.14	0.003

Table 4.  Parameters of the double-pole approximation for $ g_{VP\gamma}(q^2) $ and $ g_{VPX}(q^2) $ (with $ \varepsilon_u = \pm 5.0\times 10^{-4} $ and $ \varepsilon_d = \mp 2.9\times 10^{-3} $).

Usually, the form factor $ g_{VP\gamma}(q^2) $ is normalized to the corresponding radiative decay constant $ g_{VP\gamma}(0) $

This normalized form factor appears in Eq. (8). We can also define a similar normalized form factor for the transition induced by the X17 boson

In Figs. 3, 4, and 6 we plot the normalized form factors calculated in the CCQM and compare them with the form factors obtained using the VMD model.

Figure 3.  (color online) $ F^\gamma_{VP}(q^2) $ in our model (red, solid) and the VMD model (orange, dashed) [32, 34].

Figure 4.  (color online) $ F^X_{VP}(q^2) $ in our model (red, solid) and the VMD model (orange, dashed) [34]. The values $ \varepsilon_u = \pm 5.0\times 10^{-4} $ and $ \varepsilon_d = \mp 2.9\times 10^{-3} $ were used.

Figure 6.  (color online) $ J/\psi\to \eta_c e^+ e^- $ form factors in our model (red, solid) and the VMD model (orange, dashed) [33].

From an experimental point of view, we can directly measure the squared form factor $ |F^\gamma_{VP}(M)|^2 $, where $ M \equiv \sqrt{q^2} $ is the dilepton mass, by looking at the invariant mass distribution of the lepton pairs produced in the Dalitz decays. We then compare this measured spectrum to the prediction from QED for a point-like interaction [56]. We therefore also plot $ |F^\gamma_{VP}(M)|^2 $ in Figs. 5 and 6 for possible comparison with future experiments.

Figure 5.  (color online) Normalized form factors squared as functions of the dilepton mass $ M = \sqrt{q^2} $.

Firstly, we present our predictions for the decay widths of the Dalitz decays within the SM. The results are listed in Table 5. The largest decay width belongs to the channel $ D^{*0}\to D^0e^+e^- $, which is larger than others by 1–2 orders of magnitude. Note that these rates still have not been measured directly to date. Most often, the ratios $ R^\gamma_{ee}(V) =\Gamma(V\to Pe^+e^-)/\Gamma(V\to P\gamma) $ are of the main interest. However, the predictions given here will be important for future experiments, for instance, those at BESIII and the up-coming Super Charm-Tau factory [57].

$ D^{*+}\to D^+ $	$ D^{*0}\to D^0 $	$ D^{*+}_s\to D^+_s $	$ B^{*+}\to B^+ $	$ B^{*0}\to B^0 $	$ B^{*0}_s\to B^0_s $	$ J/\psi\to \eta_c $
4.03	122	1.32	1.69	0.58	0.53	8.20

Table 5.  CCQM predictions for the Dalitz decay widths $ \Gamma(V\to Pe^+e^-) $ in the SM (all in eV).

Next, we provide the predictions for the ratios $ R^\gamma_{ee}(V) $ and $ R^X_{ee}(V) $. Note that the ratio $ R^X_{ee}(V) $ depends on the coupling $ \varepsilon_{Q,q} $ between the X17 boson and the quarks. We will therefore briefly discuss several favored scenarios regarding these couplings. Assuming that the ATOMKI anomaly [1] is induced by the X17 boson via the process $ ^8 {{\rm{Be}}}^*+\to ^8 {{\rm{Be}}}+X(\to e^+e^-) $, the couplings between X17 and the quarks of the first generation are constrained by the condition [14, 58, 59]

The null result from the NA48/2 experimental search for $ \pi^0\to X\gamma $ requires the X17 boson to be “protophobic”, i.e. to have a suppressed coupling to the proton [60]

By combining the two conditions, one obtains the allowed ranges (i) $ \varepsilon_u \in (-4.5,-2.9) $ & $ \varepsilon_d =3.7-\varepsilon_u $ or (ii) $ \varepsilon_u \in (2.9,4.5) $ & $ \varepsilon_d =-3.7-\varepsilon_u $. Taking the central values in each range of $ \varepsilon_u $, one obtains the favored combination $ \varepsilon_u\approx \pm 3.7\times 10^{-3} $ and $ \varepsilon_d\approx \mp 7.4\times 10^{-3} $. This combination was used for the calculation in Refs. [32, 33]. In Table 6 we present our predictions for the ratios $ R_{ee}^\gamma(V) $, $ R_{ee}^X(V) $, and $ R_{ee}^{\rm tot}(V)\equiv R_{ee}^\gamma(V)+R_{ee}^X(V) $ using $ \varepsilon_u\approx \pm 3.7\times 10^{-3} $ and $ \varepsilon_d\approx \mp 7.4\times 10^{-3} $. Note that we assume universal couplings between the X17 boson with up-type and down-type quarks, i.e. $ \varepsilon_c = \varepsilon_u $ and $ \varepsilon_s =\varepsilon_b = \varepsilon_d $. We also compare with those obtained using the VMD model [32, 33]. Several comments should be addressed:

Transition	$ R_{ee}^\gamma(V) $	$ R_{ee}^X(V) $	$ R_{ee}^{ {\rm{tot}}}(V) $	Ref.	Experiment
$ D^{*+}\to D^+e^+e^- $	$ 6.47 $	$ 1.67 $	$ 8.14 $	CCQM
	$ 6.67 $	$ 1.05\pm0.07 $	$ 7.72\pm0.07 $	VMD [32]
$ D^{*0}\to D^0e^+e^- $	$ 6.45 $	$ 3.02 \times 10^{-2} $	$ 6.48 $	CCQM	$ 11.08\pm 0.90 $ [27]
	$ 6.67 $	$ 3.02\times 10^{-2} $	$ 6.70 $	VMD [32]
$ D^{*+}_s\to D^+_se^+e^- $	$ 6.50 $	$ 3.12 $	$ 9.62 $	CCQM	$ 7.2^{+1.8}_{-1.6} $ [54]
	$ 6.72 $	$ 3.10\pm0.60 $	$ 9.82\pm0.60 $	VMD [32]
$ B^{*+}\to B^+e^+e^- $	$ 4.67 $	$ 1.77\times 10^{-2} $	$ 4.69 $	CCQM
	$ 4.88 $	$ (1.91\pm0.03)\times 10^{-2} $	$ 4.90 $	VMD [32]
$ B^{*0}\to B^0e^+e^- $	$ 4.69 $	$ 3.96\times 10^{-1} $	$ 5.09 $	CCQM
	$ 4.88 $	$ 3.96\times 10^{-1} $	$ 5.40 $	VMD [32]
$ B^{*0}_s\to B^0_se^+e^- $	$ 5.63 $	$ 4.05\times 10^{-1} $	$ 6.04 $	CCQM
	$ 4.99 $	$ 4.08\times 10^{-1} $	$ 5.40 $	VMD [32]
$ J/\psi\to \eta_c e^+e^- $	$ 6.08 $	$ 2.98\times 10^{-2} $	$ 6.11 $	CCQM
	$ 6.06 $	$ 1.32\times 10^{-2} $	$ 6.07 $	VMD [33]

Table 6.  Ratios $ R_{ee}^\gamma (V) $ and $ R_{ee}^X(V) $ (in units of $ 10^{-3} $) in the CCQM and VMD model with the couplings $ \varepsilon_u =\pm 3.7\times 10^{-3} $ and $ \varepsilon_d =\mp 7.4\times 10^{-3} $.

● Our predictions for $ R_{ee}^\gamma(V) $ agree well with those in the VMD model [32, 33] within uncertainty. The central values in the two models differ by less than 5%, except for the case of $ B^{*0}_s\to B^0_se^+e^- $, for which our value is larger by approximately 11%. Note that in a recent study of charm-meson Dalitz decays [53], the authors also made use of the VMD model and obtained the values $ R_{ee}^\gamma(D^{*+})=6.44\times 10^{-3} $, $ R_{ee}^\gamma(D^{*0})=6.45\times 10^{-3} $, and $ R_{ee}^\gamma(D^{*+}_s)=6.46\times 10^{-3} $, which are almost the same as our predictions.

● Our predictions for $ R_{ee}^X(V) $ agree well with those in Ref. [32], except for the case of $ D^{*+}\to D^+e^+e^- $, for which our value is larger by 37%.

● Comparing with Ref. [33], our result for $ R_{ee}^\gamma(J/\psi) $ is almost the same as the one in [33]. However, our value for $ R_{ee}^X(J/\psi) $ is larger by a factor of two.

● The two channels $ D^{*+}\to D^+e^+e^- $ and $ D^{*+}_{s}\to D^+_{s}e^+e^- $ are most sensitive to the X17 boson contribution. However, the predicted ratio $ R_{ee}^{ {\rm{tot}}}(D^{*+}_s) =9.62\times 10^{-3} $ is slightly larger than the value $ 7.2^{+1.8}_{-1.6}\times 10^{-3} $ measured by CLEO [54].

● The X17 boson contribution to the channel $ D^{*0}\to D^0e^+e^- $ is smaller than the SM contribution by two orders of magnitude and therefore cannot account for the excess observed recently by BESIII [27].

Recently, Denton and Gehrlein performed a detailed analysis of the X17 vector boson hypothesis [61], taking into account isospin mixing and breaking effects on the rate for $ ^8 {{\rm{Be}}}^*\to \, ^8 {{\rm{Be}}}X $ [15]. By using constraints from ATOMKI anomalies in the nuclear decays of excited $ ^8 {{\rm{Be}}} $ [1], $ ^4 {\rm{He}} $ [5], and $ ^{12} {\rm{C}} $ [6], and from the NA48/2 experiment [60], they obtained the allowed regions $ |\varepsilon_u|\approx (0.5 - 0.9)\times 10^{-3} $ and $ |\varepsilon_d|\approx (2.5 - 2.9)\times 10^{-3} $ with $ \varepsilon_u \varepsilon_d <0 $. If isospin effects are ignored, the favored parameters read $ |\varepsilon_u| = \pm 5.0\times 10^{-4} $ and $ |\varepsilon_d| = \mp 2.9\times 10^{-3} $. If isospin mixing and breaking effects are taken into account, the favored values now become $ |\varepsilon_u| = \pm 9.0\times 10^{-4} $ and $ |\varepsilon_d| = \mp 2.5\times 10^{-3} $. These two parameter sets were used in Ref. [34] to calculate the X17 contributions to the $ D^{*}_{(s)}\to D_{(s)}e^+e^- $ decays within the VMD model. For comparison, we provide our predictions using these parameter sets in Tables 7 and 8. We have the following observations regarding the results in these tables:

Transition	$ R_{ee}^\gamma(V) $	$ R_{ee}^X(V) $	$ R_{ee}^{ {\rm{tot}}}(V) $	Ref.	Experiment
$ D^{*+}\to D^+e^+e^- $	$ 6.47 $	$ 0.31 $	$ 6.78 $	CCQM
	$ 6.6 $	$ 1.2\pm0.4 $	$ 7.8\pm0.4 $	VMD [34]
$ D^{*0}\to D^0e^+e^- $	$ 6.45 $	$ 5.51 \times 10^{-4} $	$ 6.45 $	CCQM	$ 11.08\pm 0.90 $ [27]
	$ 6.7 $	$ 5.6\times 10^{-4} $	$ 6.7 $	VMD [34]
$ D^{*+}_s\to D^+_se^+e^- $	$ 6.50 $	$ 0.61 $	$ 7.11 $	CCQM	$ 7.2^{+1.8}_{-1.6} $ [54]
	$ 6.7 $	$ 4.2\pm 3.1 $	$ 11\pm 3 $	VMD [34]
$ B^{*+}\to B^+e^+e^- $	$ 4.67 $	$ 0.98\times 10^{-4} $	$ 4.67 $	CCQM
$ B^{*0}\to B^0e^+e^- $	$ 4.69 $	$ 6.08\times 10^{-2} $	$ 4.75 $	CCQM
$ B^{*0}_s\to B^0_se^+e^- $	$ 5.63 $	$ 6.22\times 10^{-2} $	$ 5.69 $	CCQM
$ J/\psi\to \eta_c e^+e^- $	$ 6.08 $	$ 5.43\times 10^{-4} $	$ 6.08 $	CCQM

Table 7.  Ratios $ R_{ee}^\gamma (V) $ and $ R_{ee}^X(V) $ (in units of $ 10^{-3} $) in the CCQM and VMD model with the couplings $ \varepsilon_u =\pm 5.0\times 10^{-4} $ and $ \varepsilon_d =\mp 2.9\times 10^{-3} $.

Transition	$ R_{ee}^\gamma(V) $	$ R_{ee}^X(V) $	$ R_{ee}^{ {\rm{tot}}}(V) $	Ref.	Experiment
$ D^{*+}\to D^+e^+e^- $	$ 6.47 $	$ 0.21 $	$ 6.68 $	CCQM
	$ 6.6 $	$ 0.75\pm0.26 $	$ 7.4\pm0.3 $	VMD [34]
$ D^{*0}\to D^0e^+e^- $	$ 6.45 $	$ 1.8 \times 10^{-3} $	$ 6.45 $	CCQM	$ 11.08\pm 0.90 $ [27]
	$ 6.7 $	$ 1.8\times 10^{-3} $	$ 6.7 $	VMD [34]
$ D^{*+}_s\to D^+_se^+e^- $	$ 6.50 $	$ 0.40 $	$ 6.90 $	CCQM	$ 7.2^{+1.8}_{-1.6} $ [54]
	$ 6.7 $	$ 2.6\pm 1.9 $	$ 9.3\pm 1.9 $	VMD [34]
$ B^{*+}\to B^+e^+e^- $	$ 4.67 $	$ 8.65\times 10^{-4} $	$ 4.67 $	CCQM
$ B^{*0}\to B^0e^+e^- $	$ 4.69 $	$ 4.52\times 10^{-2} $	$ 4.74 $	CCQM
$ B^{*0}_s\to B^0_se^+e^- $	$ 5.63 $	$ 4.62\times 10^{-2} $	$ 5.68 $	CCQM
$ J/\psi\to \eta_c e^+e^- $	$ 6.08 $	$ 1.76\times 10^{-3} $	$ 6.08 $	CCQM

Table 8.  Ratios $ R_{ee}^\gamma (V) $ and $ R_{ee}^X(V) $ (in units of $ 10^{-3} $) in the CCQM and VMD model with the couplings $ \varepsilon_u =\pm 9.0\times 10^{-4} $ and $ \varepsilon_d =\mp 2.5\times 10^{-3} $.

● Our predictions for $ R_{ee}^\gamma(D^*_{(s)}) $ fully agree with those calculated in Ref. [34]. The central values in the two studies differ by less than 4%. Besides, the two studies predict the same value for $ R_{ee}^X(D^{*0}) $.

● Regarding the decays $ D^{*+}\to D^+e^+e^- $ and $ D^{*+}_{s}\to D^+_{s}e^+e^- $, our predictions for $ R_{ee}^X(D^{*+}_{(s)}) $ disagree with Ref. [34]. The predicted $ R_{ee}^X(D^{*+}_{(s)}) $ in Ref. [34] are larger than ours by approximately an order of magnitude. As a result, their total ratios $ R_{ee}^{\rm tot}(D^{*+}_{s}) $ are slightly larger than the CLEO result [54]. Meanwhile, our values for $ R_{ee}^{ {\rm{tot}}}(D^{*+}_s) $ are consistent with the last.

● The two new parameter sets lead to very small contributions of the X17 boson to the Dalitz decays of $ B^*\to B $ and $ J/\psi\to \eta_c $. And finally, the X17 boson contribution to the channel $ D^{*0}\to D^0e^+e^- $ is now even smaller than in Table 6 and cannot account for the BESIII excess [27].

Finally, in Table 9, we present our predictions for the decays $ D^{*0}\to D^0e^+e^- $, $ D^{*+}_s\to D^+_se^+e^- $, and $ J/\psi\to \eta_c e^+e^- $ using the best-fit couplings $ \varepsilon_u =6.0\times 10^{-2} $, $ \varepsilon_c =6.4\times 10^{-3} $, and $ \varepsilon_s =-2.0\times 10^{-3} $ recently obtained in Ref. [34]. Note that these values were determined by fitting the VMD predictions for the Dalitz decays of $ D_s^* $, $ D^{*0} $, $ \psi(2S) $, and ϕ to available experimental data, assuming that $ \varepsilon_u\neq\varepsilon_d\neq\varepsilon_c $ in general. With these coupling parameters, the X17 contributions to $ D^{*+}_s\to D^+_se^+e^- $ and $ J/\psi\to \eta_c e^+e^- $ are negligibly small. On the other hand, the X17 contributes largely to the decay $ D^{*0}\to D^0e^+e^- $ and can explain the BESIII excess. However, the X17 contribution to $ D^{*0}\to D^0e^+e^- $ is now approximately equal to the SM one, which seems unlikely. In other words, as already discussed in Ref. [34], there is a tension between the value of the coupling parameter $ \varepsilon_u $ determined from the ATOMKI experiments and the one from the BESIII measurement [27].

Transition	$ R_{ee}^\gamma(V) $	$ R_{ee}^X(V) $	$ R_{ee}^{ {\rm{tot}}}(V) $	Ref.	Experiment
$ D^{*0}\to D^0e^+e^- $	$ 6.45 $	5.23	$ 11.68 $	CCQM	$ 11.08\pm 0.90 $ [27]
$ D^{*+}_s\to D^+_se^+e^- $	$ 6.50 $	$ 1.87\times 10^{-3} $	$ 6.50 $	CCQM	$ 7.2^{+1.8}_{-1.6} $ [54]
$ J/\psi\to \eta_c e^+e^- $	$ 6.08 $	$ 8.90\times 10^{-5} $	$ 6.08 $	CCQM

Table 9.  Ratios $ R_{ee}^\gamma (V) $ and $ R_{ee}^X(V) $ (in units of $ 10^{-3} $) with the best-fit couplings $ \varepsilon_u =6.0\times 10^{-2} $, $ \varepsilon_c =6.4\times 10^{-3} $, and $ \varepsilon_s =-2.0\times 10^{-3} $ [34].

Electromagnetic Dalitz decays of the vector mesons $ D^{*}_{(s)} $, $ B^{*}_{(s)} $, and $ J/\psi $ were studied in the framework of the Covariant Confined Quark Model in light of new experimental data from the ATOMKI and BESIII collaborations. Predictions for the form factors and the decay widths (normalized to the corresponding radiative decays) were provided, both in the SM and in the presence of the hypothetical X17 vector boson. We considered four sets of favored coupling constants between the X17 boson and the quarks based on experimental constraints. Within the SM, we found a full agreement between our predictions and those in other theoretical studies, in which the VMD model was used to obtain the form factors. A detailed comparison of the form factors in the CCQM and in the VMD model was also presented. However, the X17 contributions to the Dalitz decays calculated in the CCQM sometimes disagree with those in the VMD model. The discrepancy can be as large as an order of magnitude. This suggests more independent theoretical calculations are needed.

Finally, among seven decay modes considered in this study, the two decays $ D^{*0}\to D^0e^+e^- $ and $ D^{*+}_s\to D^+_se^+e^- $ are most sensitive to the X17 contribution. The first has not been observed, but is expected to be measured in the near future by BESIII. Meanwhile, the second was measured by CLEO and the experimental data agrees well with our prediction, especially when the X17 contribution is included. The most interesting decay mode is $ D^{*0}\to D^0e^+e^- $, for which theoretical calculations using the CCQM and the VMD model suggest a negligibly small contribution from the X17 (of the order of $ 10^{-3}\sim 10^{-2} $ compared to the SM contribution). However, the recent BESIII data for this decay showed a $ 3.5\sigma $ excess over the SM prediction, which cannot be accounted for by the X17 alone. Independent theoretical studies and more experimental data from BESIII for this decay are therefore needed to shed more light on this curious case.

C.T. Tran and A.T.T. Nguyen thank HCMC University of Technology and Education for support in their work and scientific collaboration. C.T. Tran thanks Nestor Quintero for confirming the value $ \varepsilon_e $ used in Ref. [32]. C.T. Tran thanks Thuc Uyen T.L. and her collaborators for fruitful discussion and for providing us with the VMD values given in Table 8 (Ref. [34]). The authors thank Mauro Raggi and the unknown referee for pointing out relevant references to NA64 and PADME searches for the X17 boson.