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 a proton beam at thin strip foils of $ ^7 {{\rm{Li}}} $ and 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 an 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 $ 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}}} $; however, the inferred mass of the unknown particle remained consistent with the $ ^8 {{\rm{Be}}} $ result. A brief summary of ATOMKI's results is given in Table 1. Independent experimental attempts have been made by other collaborations to test the ATOMKI anomaly. 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 $ 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 those of ATOMKI within $ 1.5\sigma $. Meanwhile, a group at the 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 a proton beam energy of 1.225 MeV and an opening angle of 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 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 considerable attention in the particle physics community and 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 Ref. [16]). They pointed out that the X17 boson could not be a scalar because of 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 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 the 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 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 X17. The idea of using the Dalitz decays to probe for new physics beyond the standard model (SM) is not new. For example, the Dalitz decay $ J/\psi\to P\ell^+\ell^- $ has been proposed [30] as a probe for dark photon searches in the BESIII experiment, accumulating a large $ J/\psi $ data sample [31]. In the light of the ATOMKI anomaly, Castro and Quintero proposed testing 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 to search 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 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 for calculating hadronic form factors. Furthermore, disagreements exist regarding 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 X17. We follow the authors of Refs. [32, 34] and assume that X17 is a vector boson. We use several favored combinations of coupling parameters between X17 and quarks to calculate the effects of X17 on the widths of these decays and compare our results with those given in Refs. [32, 34]. Hadronic form factors required for calculating the Dalitz decays are obtained in the covariant confined quark model (CCQM) 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 understanding X17 effects and help probe 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 CCQM as a tool for hadronic calculation. Section III describes Dalitz decays in terms of Feynman diagrams, invariant matrix elements, and form factors. Here, we consider both contributions from the photon (within the SM) and 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 CCQM is a framework based on the 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 $ represents the meson-quark coupling constant, $ \Gamma_M $ represents the Dirac matrix with quantum numbers specific to the meson M, and $ F_M(x;x_1,x_2) $ represents the vertex function that effectively describes the size of the meson 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}) $ represents the mass fraction of 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 $ represents a free parameter of the model, referred to as the size parameter of the meson M. 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 $ represents the wave function renormalization constant of the meson and $ \Pi'_M $ represents the derivative of the mass function of the meson. The meson mass function $ \Pi_M(p) $ is calculated using a one-loop Feynman diagram representing the self-energy of the meson (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 considers the confinement of quarks within hadrons at low energies. These parameters are determined by fitting the predictions of the model 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 cut-off parameter (in GeV).

To estimate theoretical uncertainties within the CCQM, we employed the MINUIT program, which finds the best parameter values by minimizing the $ \chi^2 $function when fitting the model to the data. Instead of a complex error propagation analysis, we adopted a simpler, less precise method to estimate errors in our physical predictions. We proceeded with calculations using only these best-fit parameter values. Upon comparison, we noted that the fitted values of the model deviate from experimental data by approximately 5%−10%. Because our calculations of hadronic quantities were analogous to those for leptonic and electromagnetic decay constants used in the fit, we estimated 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]). Therefore, we do not repeat the calculation and only utilize the well known expression for the decay width distribution given 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 the direct photon emission from constituent quarks. Note that within the CCQM, resonance diagrams corresponding to photon emission via intermediate $ V^\prime\to\gamma^* $ transitions are also allowed (see Ref. [51]). However, because 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 the masses of the vector resonances $ \rho(770) $, $ \omega(782) $, $ \phi(1020) $, and $ J/\psi(3096) $, the contributions from resonance diagrams are negligibly small and well below the estimated inherent error of model predictions ($ \sim 10 $%). Therefore, we consider only 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.  (color online) 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 $ represent the charges of heavy and light quarks in units of e, respectively. Similarly, for the interaction between the quarks and the X17 boson,

where $ \varepsilon_{Q,q} $ represents 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 $ represents the coupling constant between the X17 boson and electron (positron), and $ T_\mu^{\gamma,X} $ represents the hadronic amplitudes of transitions $ V\to P\gamma^* $ and $ V\to PX^* $, respectively. Hadronic transition amplitudes can be parametrized in terms of the transition form factors $ g_{VP\gamma}(q^2) $ and $ g_{VPX}(q^2) $ as

Differential decay rates induced by the photon and the 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 Eqs. (13) and (14), one obtains the decay width of the Dalitz decay. The widths of 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 respect 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 because the total width of the X17 boson 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]

In addition, by assuming that the X17 vector boson decays dominantly into $ e^+e^- $,

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.

Form factors $ g_{VP\gamma}(q^2) $ and $ g_{VPX}(q^2) $ defined in Eq. (12) are calculated within the CCQM by evaluating hadronic amplitudes $ T_\mu^\gamma $ and $ T_\mu^X $. We show the calculation steps for the case of $ g_{VP\gamma}(q^2) $ below. 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 they guarantee ultraviolet convergence. The loop integration is performed using the Fock-Schwinger representation for the quark propagator

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

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

The infrared cut-off 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 calculated numerically. The expression for $ I_{Q}(m^2_{V},m^2_P,q^2) $ is

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

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} $).

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 those 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 perspective, we can directly measure the squared form factor $ |F^\gamma_{VP}(M)|^2 $, where $ M \equiv \sqrt{q^2} $ is the dilepton mass observed by looking at the invariant mass distribution of lepton pairs produced in the Dalitz decays. Then, we compare this measured spectrum to the prediction from QED for a point-like interaction [56]. Therefore, we plot $ |F^\gamma_{VP}(M)|^2 $ in Figs. 5 and 6 for possible comparisons with future experiments.

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

First, we present our predictions for the decay widths of 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 that of the others by 1–2 orders of magnitude. 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 primary interest. However, the predictions given here will be important for future experiments such as 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 Dalitz decay widths $ \Gamma(V\to Pe^+e^-) $ in the SM (all in eV).

Next, we provide predictions for 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 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 $. By 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 has been 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 $. Furthermore, we also compare them 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%. In a recent study of charm-meson Dalitz decays [53], the authors used the VMD model and obtained the values of $ 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%.

● Compared 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 the 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] considering 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 considered, 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 made 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}) $.

● For 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 $. Finally, the X17 boson contribution to the channel $ D^{*0}\to D^0e^+e^- $ is now even smaller than that 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]. These values were determined by fitting 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 $. 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. By contrast, 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 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 CCQM in light of new experimental data from the ATOMKI and BESIII collaborations. Predictions for the form factors and 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 quarks based on experimental constraints. Within the SM, we found a full agreement between our predictions and those in other theoretical studies, wherein the VMD model was used to obtain the form factors. A detailed comparison of the form factors in the CCQM and VMD model was 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 that more independent theoretical calculations are required.

Finally, among the seven decay modes considered in this study, the two decays $ D^{*0}\to D^0e^+e^- $ and $ D^{*+}_s\to D^+_se^+e^- $are the most sensitive to the X17 contribution. The first has not been observed; however, it is expected to be measured in the near future by BESIII. Meanwhile, the second was measured by CLEO, and the experimental data agree 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 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 show a $ 3.5\sigma $ excess over the SM prediction, which cannot be accounted for by X17 alone. Independent theoretical studies and more experimental data from BESIII for this decay are therefore required 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 providing 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.