Probing axion-like particles in leptonic decays of heavy mesons

Peccei and Quinn introduced a new global chiral symmetry in 1977 [1], known as $ U(1)_{\rm PQ} $ symmetry, to address the strong CP problem in QCD, that is, the absence of CP violation in the strong interactions and the electric dipole moment (EDM) of the neutron being forbidden. At energy scale $ f_a $, such a symmetry is assumed to break spontaneously, resulting in the appearance of a pseudo-Nambu-Goldstone boson, namely the axion [1−5], whose mass $ m_a $ is constrained by the relation $ m_a \sim \frac{m_\pi f_\pi}{f_a} $, where $ m_\pi $ and $ f_\pi $ are the mass and decay constants of pion, respectively. This means that if the energy scale $ f_a $ is very high, the axion should be extremely light. The condition can be relaxed if we are not limited to such a QCD axion but consider more general axion-like particles (ALPs) [6−8]. In such cases, both the PQ symmetry breaking scale $ f_a $ and the ALP mass $ m_a $ can be considered as independent parameters.

People are interested in ALPs in many aspects. Theoretically, several ALPs are predicted by string theories [9, 10] and supersymmerty [11, 12]. These particles may be crucial to the evolution of cosmology, most importantly, as a candidate of dark matter [13]. The cosmological and astronomical observations can set strict constraints on very light ALPs [14−20]. The phenomenology of ALPs have also been extensively studied at the Large Hadron Collider (LHC) [21−29]. For example, they can be produced by the decays of on-shell Higgs/Z boson [6, 29], or participate as off shell mediators in the scattering processes [30]. At high-luminosity electron-positron colliders, ALPs with masses in the MeV-GeV range can also be explored non-resonantly or resonantly [31−37], which means they can be produced directly by coupling to the charged leptons or the gauge bosons, or produced through the decays of final mesons. For example, in Ref. [35] the production of ALPs at B factories is considered, including the contributions of $ e^{+}e^{-} \to \gamma a $ and $ e^+e^-\to\Upsilon\to\gamma a $. Recently, the BESIII experiment [38] used a data sample of $ \psi(3686) $ to obtain the upper limits of the branching fraction of the decay $ J/\psi \to \gamma a $ and set constrains on the coupling $ g_{a\gamma\gamma} $ in the mass range of $0.165\; {\rm GeV} \leq m_a \leq 2.84\; {\rm GeV}$. The fixed target experiments are also promising methods for searching for ALPs [38−43], as the detectors can extend tens of meters and are suitable for detecting long-lived particles. ALP production through K meson decays are usually investigated in such cases [38−40].

Except for the methods mentioned above, the decays of heavy-flavored mesons also provide a way to probe ALPs. On one hand, a large number of bottomed and charmed mesons have been produced at B factories and the τ-charm factory, respectively, and on the other, a larger range of $ m_a $ can be looked into in heavy meson decays compared with that of K meson. We will focus on the decay processes $ h^- \to \ell^- \not {E} $, where $ h^-=D^-,\; D_s^-, \; B^- $, or $ B_c^- $; $ \ell^-=e^-,\; \mu^- $, or $ \tau^- $; $ \not {E} $ represents the missing energy. In the Standard Model, $ \not {E} $ is transported by the antineutrino (see Fig. 1), and the corresponding partial width is

Figure 1.  Diagram of two-body leptonic decays of a charged heavy meson.

where $ G_F $ is the Fermi coupling constant, $ f_h $ is the decay constant of h meson, $ V_{Qq} $ is the CKM matrix element, M and $ m_\ell $ are the masses of the meson and the lepton, respectively, and the neutrino mass is assumed to be zero. If an ALP can also be produced and has sufficient lifetime to escape the detector, we will have $\Gamma(h^-\to \ell^-\not {E})=\Gamma(h^-\to \ell^-\bar{\nu}_{\ell})+\Gamma(h^-\to \ell^-\bar{\nu}_{\ell}a)$. Then the experimental results of $h^-\to \ell^- \not {E}$, which are presented in Table 1, can be used to constrain the coupling between the ALP and SM particles.

Channel	Experiment values
$ {\mathcal B}(B^-\to e^- \not {E}) $	$ \le 9.8 \times 10^{-7} $ [44]
$ {\mathcal B}(B^-\to \mu^- \not {E}) $	$ (5.3 \pm 2.0 \pm 0.9) \times 10^{-7} $ [45]
$ {\mathcal B}(B^-\to \tau^- \not {E}) $	$ (7.2 \pm 2.7 \pm 1.1 ) \times 10^{-5} $ [46]
$ {\mathcal B}(D_{s}^-\to e^- \not {E}) $	$ \le 0.83 \times 10^{-4} $ [47]
$ {\mathcal B}(D_{s}^-\to \mu^- \not {E}) $	$ (0.5294 \pm 0.0108 \pm 0.0085) \times 10^{-2} $ [48]
${\mathcal B}(D_{s}^-\to \tau^- \not {E})$	$ (5.44 \pm 0.17 \pm 0.13 ) \times 10^{-2} $ [49]
${\mathcal B}(D^-\to e^- \not {E})$	$ \le 8.8 \times 10^{-6} $ [50]
${\mathcal B}(D^-\to \mu^- \not {E})$	$ (3.71 \pm 2.7 \pm 1.1) \times 10^{-3} $ [51]
${\mathcal B}(D^-\to \tau^- \not {E})$	$ (1.20 \pm 0.24 \pm 0.12 ) \times 10^{-5} $ [52]

Table 1.  Experiment results of $ {\mathcal B}(B^- \to \ell^-\not {E}) $, $ {\mathcal B}(D^-\to $$ \ell^-\not {E}) $, and $ {\mathcal B}(D_s^-\to \ell^-\not {E}) $.

This study is inspired by Refs. [53−56] where the processes $ h^- \to \ell^- \bar{\nu}_{\ell}a $ with a being an ALP were studied. Beyond the considerations for the case of a pseudoscalar, Ref. [53] also considered the possibility of a being a vector particle. Here, we will also probe such processes, but with three main differences. First, the Hadron transition matrix elements will be calculated using the Bethe-Salpeter (BS) method. This method has been widely used to study the weak decays [57−59], strong decays [60], and electromagnetic decays [61] of heavy mesons, and the theoretical results are consistent with the experimental data. Second, for the coupling of ALP with quarks, we consider the interference effects of different diagrams and extract the allowed range of parameters in a more general way. Third, we will explore ALP production from the decays of the double heavy meson $ B_c $, that is $ B_c^- \to \ell^- \bar{\nu}_{\ell}a $, which could give the constraint in a larger ALP mass range.

The remainder of the paper is organized as follows: In Sec. II, we outline the effective Lagrangian, which describes the interaction between ALPs and charged fermions. Then, we calculate the Hadronic transition matrix elements of the $ h^-\to \ell^-\nu_{\ell} a $ processes using the BS method. In Sec. III, we apply the experimental results to extract the allowed regions of the parameters, which are then used to calculate the upper limits of the branching fractions of the $ B_c $ decay processes. Finally, we present the summary in Sec. IV.

Generally, the ALPs can interact with the fermions and gauge bosons. Herein, in the case of leptonic decays of heavy mesons, we only consider the tree level diagrams shown in Fig. 2, which represent the interaction of the ALPs with SM fermions. The contributions of gauge bosons are of a higher order and assumed to be neglected. The effective Lagrangian, which describes flavor conserving processes, comprises dimension-five operators [53−56]

Figure 2.  Feynman diagrams of $ h \rightarrow \ell\bar{\nu}_\ell a $. (a), (b), (c), and (d) represent the ALP couples with the quark, antiquark, charged lepton, and antineutrino, respectively.

where c_i is the effective coupling constant with the index i extending over all the fermions. Similar with Refs. [54−56], in the last step, we have used the equation of motion to rewrite the expression. This corresponds to the use of another operator basis, which will lead to different upper limits of c_i. However, when we calculate the decay channels of $ B_c $ meson, two methods give quite similar results. The dependence of the fermion mass $ m_i $ indicates that the ALP does not couple directly with the neutrinos and Fig. 2(d) does not need to be considered, which will simplify the calculation.

The Feynman amplitude corresponding to each diagram can be written using the meson wave function. Here, the instantaneous BS wave function is appropriate for the calculation, as we are considering the heavy flavored meson, which can be seen as a two-body bound state including at least one heavy quark (or antiquark). This kind of wave function has been extensively used to study the decay processes of heavy mesons. For the pseudoscalar meson, its wave function has the form [62]

where P is the meson momentum and q is the relative momentum between the quark and antiquark; $q^\mu_\perp \equiv q^\mu- \dfrac{P \cdot q }{M^2} P^\mu = (0,\; \vec q)$; $ f_i $ is the radial wave function, which can be achieved by solving the eigenvalue equations numerically (e.g., the detailed results for the $ B_c $ meson can be found in Ref. [63]).

The amplitude corresponding to Fig. 2(a) can be written as

where $ c_Q $ is the coupling constant; $ m_Q $ is the mass of quark Q; $ N_c=3 $ is the color factor; $ k_a $, $ k_\ell $, and $ k_{\bar\nu} $ are the momenta of ALP, charged lepton, and antineutrino, respectively; $ p_Q $ is the momentum of quark Q, which is related to P and q, as $p_Q = \dfrac{m_Q}{m_Q+m_q} P+q$, with $ m_q $ being the mass of the antiquark. To obtain Eq. (4), we have made an approximation, namely $p_Q \approx \dfrac{m_Q}{m_Q+m_q} P+q_\perp$, so that the denominator is independent of the time component of the relative momentum q. After integrating out $ q_\perp $, $ \mathcal{M}_a $ becomes

where we have used $A_Q = \dfrac{c_Q m_Q G_F V_{qQ}}{f_a \sqrt{2}}$ for short; $ F_1 $ and $ F_2 $ are the form factors of the hadronic transition matrix element, which are functions of the squared momentum transition $ (P-k_a)^2 $.

Similarly, the Feynman amplitude corresponding to Fig. 2(b) can be written as

where $ c_q $ is the coupling constant; $ m_q $ and $ p_q $ are the mass and momentum of the antiquark, respectively; and $ p_q $ has the form $p_q\approx \dfrac{m_q}{m_Q+m_q} P-q_\perp$. After integrating out $ q_\perp $, we get

where $A_q=\dfrac{c_q m_q G_FV_{qQ}}{f_a \sqrt{2}}$; $ G_1 $ and $ G_2 $ are form factors.

When calculating the integral, one must be careful, since as the relative momentum $ \vec q $ changes, the denominator of the propagator may be zero, which will lead to a nonvanishing imaginary part of the form factors [64]. Specifically, we can write the propagator as

where θ is the angle between $ \vec q $ and $ \vec k_a $. a and b are expressed as

where $ E_a $ is the energy of the ALP. We will first integrate out $ \cos\theta $ analytically, and then integrate out $ |\vec q| $ numerically. As $ \vec q $ changes, $ a+b $ and $ a-b $ may have opposite signs, which means the pole can exist with a specific value of θ.

The Feynman amplitude corresponding to Fig. 2(c) can be written as

where $ c_\ell $ is the coupling constant. The hadronic part of the amplitude is proportional to the meson momentum [65], with the proportionality being the decay constant $ f_h $

Using Eq. (10) and defining $A_\ell = \dfrac{c_\ell m_\ell G_FV_{qQ}}{f_a \sqrt{2}}$ for short, $ \mathcal{M}_c $ is simplified to

The total Feynman amplitude is $ \mathcal M=\mathcal M_a +\mathcal M_b +\mathcal M_c $. The partial width of such a decay channel is achieved by finishing the three-body phase space integral

where the upper and lower limits have the following forms:

In this section, we will first use the formulas obtained above to calculate the limits of the coupling constants, and then the results will be applied to investigate the similar decays of the $ B_c $ meson. We will focus on three different scenarios: First, the ALP couples only with a single fermion, namely the quark, the antiquark, or the charged lepton; second, the ALP couples with quarks and antiquarks but not with leptons; third, the ALP couples with all the charged fermions with the same coupling constant.

During the calculation process, the following numerical values of three groups of physical quantities are adopted: (1) The masses of constituent quarks are from Ref. [66]: $ m_u = 0.305 $ GeV, $ m_d = 0.311 $ GeV, $ m_s = 0.5 $ GeV, $ m_c = 1.62 $ GeV, and $ m_b = 4.96 $ GeV; (2) the relevant CKM matrix elements are from Particle Data Group (PDG) [67]: $ |V_{ub}| = 3.82 \times 10^{-3} $, $ |V_{cd}| = 0.211 $, $ |V_{cs}| = 0.975 $, and $ |V_{cb}| = 41.1 \times 10^{-3} $; (3) the decay constant of D, $ D_s $, and B mesons are the Lattice QCD results [68]: $ f_D = 212.6 $ MeV, $ f_{D_s}=249.9 $ MeV, and $ f_B = 188 $ MeV.

We separately consider the coupling between the ALP and each charged fermion. For example, by setting $ c_q $ and $ c_\ell $ to zero, we get the branching fraction of $ h^-\to\ell^-\bar\nu a $ with $ c_Q/f_a $ and $ m_a $ being unknown parameters. Then, using the experimental results in Table 1, we get the upper limits of $ c_Q/f_a $ as functions of $ m_a $. The same is true for the other two cases. The results are shown in Fig. 3. One notices that the curves raise rapidly when $ m_a $ is sufficiently, which is mainly owing to the reduced phase space. The most stringent upper limits for different coupling constants result from different decay channels. For example, Fig. 3(a) shows that the most strict constraints of $ c_{b}/f_a $ and $ c_{c}/f_a $, which are of the order of $ 10^2\; {\rm TeV}^{-1} $ when $ m_a $ is less than $ 1 $ GeV, result from $ B^-\to\mu^-\bar\nu_\mu a $ and $ D_s^+\to e^+\nu_e a $, respectively. One can also notice that the upper limit of $ c_e/f_a $ in Fig. 3(c) is much higher than those of $ c_{Q(q)}/f_a $ set up by the same channels, which is because of the small $ m_e $.

Figure 3.  (color online) The upper limits of the coupling constants (a) $ c_Q / f_a $, (b) $ c_q / f_a $ and (c) $ c_{\ell} / f_a $ derived from the leptonic decays of D, $ D_s $, and B mesons.

Using the above results, we can impose restrictions on $ \mathcal B(B_c^-\to \ell^- \bar{\nu}_\ell a) $, which are presented in Fig. 4. We can see that when $ m_a $ is less than 1 GeV, the branching ratio upper limits for ALP-quark coupling cases are of the order of $ 10^{-4} $. For the ALP-tau coupling case, the order of magnitudes is $ 10^{-2} $, which is much larger than those of the ALP-$ e/\mu $ coupling cases. We also notice that all the curves, except the one corresponding to $ B_c^-\to\tau^-\bar\nu_\tau a $ in Fig. 4(a), show an ascending trend. This is because the maximum value of $ m_a $ permitted in the $ B_c $ decay is larger than that in the same decays of the B and $ D_{(s)} $ mesons. For example, when we use $ c_{b}/f_a $ from $ B^-\to\mu^-\bar\nu_\mu a $ to constrain the branching ratios of $ B_c^-\to\mu^-\bar\nu_\mu a $, the result blows up before $ m_a $ reaches its maximum. The appearance of kinks in Fig. 4(c) is because we used $ c_{e/\mu}/f_a $ from different channels when $ m_a $ takes different values.

Figure 4.  (color online) The upper limits of $ \mathcal{B}(B_c^- \to \ell^- \bar{\nu}_\ell a) $ in scenario 1.

In this scenario, we consider the case that the ALP couples with quarks but not with leptons. Two assumptions are made. First, all the coupling constants are assumed to be real numbers, and no additional relative phase is introduced. Second, the coupling constants of the ALP and light quarks, denoted by $ c_q $, are assumed to be equal, while for those of heavy quarks, $ c_Q $ ($ Q=c,\; b $) may have different values. We let $ c_q $ be a free parameter, then we use the experimental results of B and $ D_{(s)} $ mesons to provide constraints on $ c_b $ and $ c_c $. Specifically, we first write the following expression:

where the three terms on the right side of the last equation represent the contributions of Fig. 2(a), 2(b), and their interference, respectively, and $ \tilde\Gamma $ denotes the decay width devided by the coupling constant sqaured, which can be calculated theoretically.

The experimentally allowed region for $ c_Q $ and $ c_q $ is determined by Eq. (13). We investigate all the decay channels, and find that the decays of $ B^- $ and $ D^- $ mesons give the most stringent constraint. The area of the parameter space depends on $ m_a $. In Fig. 5, we show the results when $ m_a $= 0, and 1.6 GeV, respectively. One can see that when $ m=0 $ GeV, the two areas are comparable; however, when $ m=1.6 $ GeV, the area for $ Q=c $ (just partially plotted) is larger than that for $ Q=b $.

Figure 5.  (color online) The experimentally allowed region (inside the ellipse) of the coupling constants. The green and blue areas correspond to $ Q=c $ and $ Q=b $, respectively.

Next, we scan the experimentally allowed parameter area in Fig. 5 and calculate the limits of $ \mathcal{B}(B_c^- \to \ell^- \bar{\nu}_\ell a) $. Because both $ c_Q $ and $ c_q $ changing the sign does not affect the results, we only need to consider the region with $ c_q\ge 0 $. Our strategy is as follows: First, by solving Eq. (13), we obtain the boundary values of $ c_Q $ for a fixed $ c_q $,

where we have defined $ C_{Q(q)}=c_{Q(q)}/f_a $. The condition $C_q^2 \le \dfrac{4 \tilde\Gamma_Q \delta\Gamma}{4 \tilde\Gamma_Q \tilde\Gamma_{q}-\tilde\Gamma_{Qq}^2}$ should be satisfied to make sure the quantity under the square root nonnegative. Second, we scan $ C_b $ and $ C_c $ together for a fixed $ C_q $ whose allowed region is determined by the small ellipse in Fig. 5.

Finally, we calculate the upper limits of $ \mathcal{B}(B_c^- \to \ell^- \bar{\nu}_\ell a) $ as functions of $ C_q $ with a specific $ m_a $. The results are shown in Fig. 6(a) and (c). Here, we only present the result for the positive $ C_q $, and that of the negative $ C_q $ is symmetrical with it about the vertical axis. As $ m_a $ increases, the experimental allowed region of $ C_q $ shrinks, and the upper limit of $ \mathcal{B}(B_c^- \to \ell^- \bar{\nu}_\ell a) $ roughly increases with $ m_a $. This can be understood from Fig. 5, where as $ m_a $ changes from 0 GeV to 1.6 GeV, the allowed range of $ c_c/f_a $ becomes quite large, resulting in a larger upper limit. The interesting thing is that for some range of $ c_q $, a nonvanishing lower limit of the branching fraction exists (see Fig. 6(b) and (d)). From Fig. 5, we can see that the two ellipses are tilted; therefore, at some $ c_q $, $ c_b $ and $ c_c $ cannot be zero simultaneously, which leads to the nonvanishing $ \mathcal{B}(B_c^- \to \ell^- \bar{\nu}_\ell a) $.

Figure 6.  (color online) The upper and lower limits of $ \mathcal B(B_c \rightarrow \ell \bar{\nu}_\ell a) $ in scenario 2.

In this scenario, we assume that the ALP couples with all the charged fermions with the same coupling constant $ c_q=c_Q=c_l=c $, the upper limits of which, as functions of $ m_a $, are shown in Fig. 7(a). We can see the $ B^-\to\mu\bar\nu_\mu a $ channel provides the most stringent restriction. We also notice that the curves have undulation around 1 GeV. We take the $ B^-\to\mu\bar\nu_\mu a $ channel as an example to understand why this happened, and plot the contributions of three Feynman diagrams in Fig. 2 and the corresponding interference terms. The effect of the lepton can be neglected owing to its small mass, while the coupling between the ALP and quarks provides the main contribution, especially the interference term, which has a sizeable negative value when $ m_a<2 $ GeV. This causes the total result to be undulant around $ m_a=1 $ GeV, which is also transferred to the coupling constant.

Figure 7.  (color online) (a) shows the upper limits of the coupling constants in scenario 3. (b) Shows the contributions to $ \mathcal{B}(B^- \to e^- \bar{\nu}_e a) f_a^2 /c_{i}^{2} $ by each Feynman diagram and the interference terms.

Using the upper limit of the coupling constant given by the $ B^-\to\mu\bar\nu_\mu a $ channel, we obtain the constraints of the branching fraction of $ B_c^- \to \ell^- \bar{\nu}_{\ell} a $, which are shown in Fig. 8(a). We can see, for $ \mathcal{B}(B_c^- \to e^- \bar{\nu}_e a) $ and $ \mathcal{B}(B_c^- \to \mu^- \bar{\nu}_{\mu} a) $, that the upper limit is approximately $ 10^{-6} $ when $ m_a $ is less than 1 GeV, and then it keeps increasing until it is larger than one. This is because the phase space of the B decay is small compared with that of the $ B_c^- $ decay when $ \ell=e $ or μ. For $ \mathcal{B}(B_c^- \to \tau^- \bar{\nu}_e a) $, the upper limit is approximately $ 10^{-5} $ when $ m_a<1 $ GeV, and then fluctuates with the peak value of $ 5.4\times 10^{-4} $ at $ m_a=3.4 $ GeV. Finally, the branching fractions reach zero owing to the vanishing phase space. Comparing these results with the branching fractions of $ B_c^-\to\ell^-\bar\nu_\ell $ in the SM is noteworthy. Using Eqs. (1) and (10), we get $\mathcal B(B_c^-\to \ell\bar\nu_\ell)= 2.2\times 10^{-9},\; 9.5\times 10^{-5}$, and $ 2.2\times 10^{-2} $ for $ \ell=e,\; \mu $, and τ, respectively. One can see that $ \mathcal B(B_c^-\to e^-\bar\nu_e) $ is much smaller than the upper limit of $ \mathcal B(B_c^-\to e^-\bar\nu_e a) $ for the chiral suppression; conversely, $ \mathcal B(B_c^-\to \tau^-\bar\nu_\tau) $ is much larger than the upper limit of $ \mathcal B(B_c^-\to \tau^-\bar\nu_\tau a) $.

Figure 8.  (color online) (a) shows the upper limit of $ \mathcal{B}(B_c^- \to \ell^- \bar{\nu}_{\ell} a) $ in scenario 3. (b), (c), and (d) show the upper limits of the differential branching fraction with $ \ell=e $, μ, and τ, respectively.

In Fig. 8(b), (c), and (d), we present the upper limits of the differential branching fractions as functions of the lepton energy. For $ \ell $ being e or μ, the results are similar owing to the small lepton mass; however, for $ \ell=\tau $, the distribution shape is quite different. We can see that the peak values shrink as $ m_a $ increases. Experimentally, the detection of the unmonoenergetic charged lepton in the decay $ h^-\to\ell^-\not {E} $ may indicate the existence of ALP, as it has the largest experimentally allowed probability to find the charged lepton around the peak.

In conclusion, we have studied ALP production through the processes $ h^- \to \ell^- \bar{\nu}_\ell a $. The instantaneous BS wave functions of the heavy mesons are applied to compute the branching fractions of such decay channels. We adopt three scenarios, that is, the ALP coupling only to one charged fermion, the ALP coupling only to quarks, and the ALP coupling to all the charged fermions with the same coupling constant. In each scenario, by comparing the theoretical and experimental results, we obtain the upper limits of the coupling constants, which are then used to calculate the upper limits of the branching fractions of the $ B_c^- \to \ell^- \bar{\nu}_{\ell} a $ channel. For the second scenario, we also obtain the nonzero lower limit of the branching fraction at some range of $ c_q $. This study is expected to be helpful for the future detection of ALPs via heavy meson decays.