The study of light invisible particles in $B_c$ decays

In this paper, we study the light scalar and pseudoscalar invisible particles in the flavor changing neutral current processes of the $B_c$ meson. Effective operators are introduced to describe the couplings between quarks and light invisible particles. The Wilson coefficients are extracted from the experimental results of the $B$ and $D$ mesons, which are used to predict the upper limits of the branching fractions of the similar decay processes for the $B_c$ meson. The hadronic transition matrix element is calculated with the instantaneously approximated Bethe-Salpeter method. The upper limits of the branching fractions when $m_\chi$ taking different values are presented. It is found that at some region of $m_\chi$, the channel $B_c\to D_s^{(\ast)}\chi\chi$ has the largest upper limit which is of the order of $10^{-6}$, and for $B_c\to D_s^\ast\chi\chi^\dagger$, the largest value of the upper limits can achieve the order of $10^{-5}$. Other decay modes, such as $B_c\to D^{(*)}\chi\chi^{(\dagger)}$ and $B_c\to B^{(*)}\chi\chi^{(\dagger)}$, are also considered.


I. INTRODUCTION
The standard model (SM) is extremely successful. However, it is considered to be an effective field theory which is valid only up to certain energy scale. For example, it will be invalid at the Planck scale, with gravity giving large contribution. Far below that, there are many arguments supporting that new physics (NP) will appear at the TeV scale. The NP can show itself as the missing energy in the collision at the pp or e + e − colliders. If we assume the possible new particle to be the candidate for the dark matter (DM), the high energy collision will provide a powerful way to detect such particles. Among the DM candidates, the weakly interacting massive particle (WIMP), which appears in many theoretical models, has attracted extensive attention (see [1] for reviews). The WIMP annihilation cross section is constrained by the observed dark matter density, which sets the lower bound of the WIMP mass to a few GeV (the so-called Lee-Winberg limit [2]). However, this result is modeldependent. If the DM is nonfermionic and the weak mass scales or weak interactions are not assumed [3], this constraint can be relaxed, and more lower mass, such as a few keV, will be possible. Theoretically, this kind of light dark matter (LDM) can have different spins, for example, it can be a scalar particle [4], sterile neutrino [5], or hidden vector DM [6]. The MeV-scale LDM is proposed [7,8] to explain the unexpected emission of 511 keV photons from the galaxy center. Experimentally, the parameter space for the WIMP with mass larger than several GeV has been severely constrained by the recent experiment [9], which also provides a strong motivation for the study of the sub-GeV dark matter.
The LDM emission from the heavy meson decays is an interesting approach for such studies. Phenomenologically, the LDM of some hidden sector can weakly interact with the SM fermions through different ways. For example, it can couple directly to the Higgs boson [10,11]. Or there are some connectors with quantum numbers of both SM and hidden sectors. Such connector can be a chiral fermion [12] or a dark gauge boson [13].
At the energy level of heavy mesons, these processes will be greatly suppressed by the large mass in the propagator of the connector or by the small coupling constant between the connector and the SM fermions. By a model-independent way, we can introduce an effective Lagrangian to describe phenomenologically the interaction between the LDM and SM fermions. This method has been extensively used in Refs. [14][15][16][17][18][19] to study the flavorchanging neutral current (FCNC) processes of K, D, and B mesons. The SM background comes from the decays with νν in the final states, which has small branching fraction and makes the detection of NP possible. The difference between the experimental results for The rest of the paper is organized as follows: In Sec. II, we first construct the effective Lagrangian which describes the coupling between quarks and light dark matter particles.
Then by comparing the theoretical and experimental results, we extract the upper limits of the Wilson coefficients. In Set. III, these upper limits are used to constrain the branching fractions of the decay channel B c → hχχ with h and χ being the final meson and dark matter particle, respectively. Finally, we give the summary and perspective in Sec. IV.

II. EFFECTIVE OPERATORS
A. χ is a scalar At the quark level, the LDM emission processes of the heavy meson can be described by the effective Lagrangian [16], where q and q f are the Dirac spinor fields of the initial and final quarks, respectively; g s1 and g s2 are the phenomenological coupling constants. This Lagrangian is model-independent.
And for specific models, the four-particle vertex may be generated at the tree or loop level [14][15][16][17][18] by introducing other new particles. In this work, we will not focus on any specific model, but consider the FCNC processes of B c meson induced by such effective operators. Theoretically, there are many studies [23][24][25][26] of the FCNC processes of the B c meson, while the corresponding detection is still missing. So we cannot use the experimental data of B c meson to set constraints on the coupling constants. Our strategy is in the opposite direction. That is, the allowed-region of the coupling constants from other processes are used to constraint the branching ratios of the B c decays. Experimentally, there are data for such decays of B and D mesons. The corresponding channels are bounds for their branching ratios are listed in Table I. Within the SM, the missing energy / E represents the νν pair, and the branching fractions are calculated in Refs. [27][28][29][30]. The difference between theoretical predictions and experimental bound allows the existence of NP. Here the NP processes are described by the Feynman diagrams in Fig. 1.
BR(B ± → π ± / E) < 14 BR(B ± → π ± νν) = 9.7 ± 2.1 BR(B ± → π ± χχ) < 6.4 For the B − → π − (K − )χχ processes, only the scalar current gives contribution to the transition amplitude, which can be written as where P and P f are the momenta of the initial or final mesons, respectively; m q and m q f are the masses of quarks; s is defined as (P − P f ) 2 . In the second step, the equation of motion is used. The hadronic transition matrix is parameterized as the form factors f + and f 0 . Here we adopt the results of the QCD light-cone sum rules (LCSR) [34], where the form factors are constructed as The corresponding parameters are presented in Table II.   TABLE II. Parameters in the form factors of the B → π(K) processes [34]. For the B − → ρ − (K * − )χχ processes, only the pseudoscalar current gives contribution to the transition amplitude, which has the form, where is the polarization vector of the final meson; M and M f are the masses of the initial and final mesons, respectively. A 0 , A 1 , A 2 , and A 3 are form factors. Within the LCSR method [35], they have the forms And the related parameters are given in Table III.   TABLE III. Parameters in the form factors of the B → ρ(K * ) processes [35]. By finishing the three-body phase space integral, we get the branching ratios where λ(x, y, z) = x 2 + y 2 + z 2 − 2xy − 2xz − 2yz is the Källen function; m χ is the mass of the LDM; θ is the angel between the three-dimensional momenta P χ and P f in the momentum center frame of final two LDMs; Γ B − is the total width of B − meson; Ω = 2 originates from the final two LDMs being identical particles.
For the annihilation processes of B 0 , D 0 mesons, that is B 0 (D 0 ) → χχ, only the pseudoscalar current contributes to the decay amplitude, which has the form, where f M is the decay constant of the initial meson, which has the values: GeV and f D 0 = 0.230 GeV [36]. By finishing the two-body phase space integral we get the branching ratio, By comparing the theoretical predictions and the experimental upper bounds (the fourth column of Table 1) of the branching ratios for these decays, we can set the upper bounds for the effective coupling constants g s1 and g s2 with specific mass of the LDM. The results are shown in Fig. 2. One can see that as m χ increasing, the upper limits of the effective coupling constants get more and more larger. The reason is simple: larger m χ means more suppression from the phase space. So from these decay channels we can set more stringent upper limits for the effective coupling constants when 2m χ is not close to the threshold. The B → π(ρ)χχ channel gives smaller bound of |g s1 | 2 (|g s2 | 2 ) compared with the B → K(K * )χχ channel. The B → χχ mode gives the most stringent upper bound of |g s2 | 2 , because the two-body phase space is larger than the three-body case. We also present the result from D 0 decay, which is larger due to its smaller mass.

B. χ is a pseudoscalar
If the LDM is a pseudoscalar, χ and χ † represent different fields. The effective Lagrangian which describes the FCNC processes q → q f χχ † has the form [16], where we have used the definition The last two terms disappear when χ is a scalar.
For the decays of B meson, when the final meson is a pseudoscalar, the second and the fourth terms in Eq. (9) will not contribute to the decay. The FCNC process can be induced by the scalar or vector current, and the transition amplitude has the form where P 1 and P 2 are the four-dimensional momenta of χ and χ † , respectively; h − is π − or K − . The hadronic transition matrix element is parameterized the same as that in Eq. (2) or Eq. (4), and the form factors are expressed in Eq. (3) or Eq. (5).
The transition amplitude receives the contribution from two terms in the effective Lagrangian, and the partial width can be written as Here we have defined , which are independent of the effective coupling constants. The numerical calculation shows that Γ 13 and Γ 31 are much smaller than Γ 1 and Γ 3 . So the contribution of the cross terms can be safely neglected. Comparing the theoretical predictions and the experimental upper bound of these channels, we give the possible relations of the modulus square of the effective coupling constants, which are presented in Fig. 3. In this figure, the area below the colored line is allowed experimentally with a specific mass of the LDM. One can see that as m χ increasing, the allowed region gets larger and larger.  If the final meson is a vector, the situation is a little more complicated, because this time the decay processes can be induced by the second, the third, and the fourth operators in the effective Lagrangian. The transition amplitude is where h * − represents ρ − or K * − . Here we need to consider two kinds of hadronic transition matrix elements. h * − |q f γ µ γ 5 q|B − is parameterized the same as Eq. (4).
where the form factor V (s) has the form [35], with the parameters listed in Table IV.  [35]. The relationship between three effective couplings can be achieved by comparing the theoretical results and the experimental upper limits. Numerical calculation indicates the cross terms can also ben neglected. In Fig. 4, we show that the experimentally allowed region is that under the colored plane which corresponding a specific mass of the LDM.
where P is the momentum of the meson; p 1 and the p 2 are the momenta of the quark and antiquark, respectively; m 1 and m 2 are the masses of the quark and antiquark, respectively; q is the relative momentum between quark and antiquark; χ P (q) is the BS wave function; V is the interaction kernel.
For B c meson, we can safely make an instantaneous approximation for V , that is V (P, k, q) ≈ V (P, k ⊥ , q ⊥ ), where q ⊥ = q − P ·q √ P 2 P , and the same is for k ⊥ . By defining the Salpeter wave function ϕ(q ⊥ ) = i dq 0 2π χ P (q), we reduce Eq. (14) to the three-dimensional form, which can be solved numerically. ϕ(q ⊥ ) is constructed from P , q ⊥ , Dirac gamma matrices, and some scalar function of q 2 ⊥ . We take the 0 − and 1 − states as examples, whose Salpeter wave functions are [38] where ϕ ++ (q ⊥ ) = Λ + where G F is the Fermi coupling constant; α is the fine structure constant; θ W is the Weinberg angle; V q 1 q 2 is the Cabibbo-Kobayashi-Maskawa (CKM) matrix element; X l (x t ) is the Inami-Lim function [40], which has the form with x t = m 2 t /M 2 W . The transition amplitude is The hadronic transition matrix element is calculated by Eq. (16). The branching fraction is achieved by finishing the three-body phase space integral, which is presented in Table V to compare with the results of other models.
There are also the B c → B u νν processes, which is induced by c → u at the quark level.
The Feynman diagrams for such channels are given in Fig. 6. The corresponding effective Lagrangian is In the above equation, we have defined X l (x q ) =D(x q , y l )/2, and the Inami-Lim function D(x q , y l ) is expressed as [40] D(x q , y l ) = 1 8 where we have used x q = m 2 q /M 2 W and y l = m 2 l /M 2 W . The transition amplitude has the form and respectively, where h ( * ) can be D ( * ) , D ( * ) s , or B ( * ) , and the hadronic transition matrix elements are calculated with Eq. (17). By finishing the three-body phase space integral, we get the decay widths expressed as the product of the squared effective coupling constant |g si | 2 and the quantity Γ i . Γ i is independent of the coupling constants, and can be calculated by taking a specific value of m χ . In Fig. 7, we plot them as functions of m χ . One can see that they all decrease when m χ gets larger, because the phase space gets smaller. With the same value of m χ , Γ 1 from the B c → P χχ channels is larger than Γ 2 from the B c → V χχ channels due to the different effective vertex in the amplitude. For Γ 1 , the B c → D s channel gives a larger result than that of the B c → D channel when m χ < 1.9 GeV. When m χ gets even larger, the phase space suppression will be important. For Γ 2 , this turning point is about 1.53 GeV.
We also notice that Γ i of the c → u processes is two orders of magnitude less than that of the b → d(s) processes. This comes from both the smaller phase space and smaller m q for the former case. respectively (see Fig. 2(a)). For the B c → D * s χχ channel, the result of B → K * χχ is used. For the B c → D * χχ case, there are two processes available to set the upper limits, namely B → ρχχ and B → χχ. The later one gives more stringent constraint, which is applied here. As there is no experimental result now for the the D → π + / E, we cannot give the constraints for the B c → Bχχ channel. The |g s2 | 2 extracted from the D 0 → χχ is used set the upper limit for the branching ratio of the B c → B * χχ channel. The upper limits of the branching ratios of B c → P (V )χχ channels are of the order of ( * ) s χχ, a specific feature appears. When m χ is less than about 1.5 GeV, the branching ratios increases slowly with m χ ; after that, the branching ratios decreases rapidly to zero. It is the result of a combination of the increasing |g si | 2 and decreasing Γ i . One notices that in Fig. 7(a) the Γ of the D * (s) case is smaller than that of the D (s) case, however, the branching ratios of the former are several times larger than that of the later, because the experimental upper bound of B and D mesons in Table I  If the LDM is a pseudoscalar, the Lagrangian L 2 is applied. When the final meson is a pseudoscalar, the transition amplitude has the form and for the vector meson case, the transition amplitude can be written as The hadronic transition matrix elements are also calculated with Eq. (17). But this situation is more complicated, because there are two or three operators contribute to the decay.
Similar to the subsection II.B, we will neglect the cross terms, and keep the ones proportional to |g pi | 2 . These terms are named as Γ i , which are independent of the effective coupling constants but depend on the mass of the light dark matter.
To compare the contribution of different terms, we calculate them by giving a specific value of m χ . The results are presented in Fig. 9. We can see all the Γ i s are decreasing with m χ , which is due to the suppression of phase space. When the final meson is D or D s , Γ 1 is large than Γ 3 , while for B, the situation is very different. When the final meson is a vector, Γ 2 and Γ 4 are close to each other, but both larger than Γ 4 . The approach we get the upper limits of the decay width is as follows. When m χ is given, there is an experimental allowed region for the effective coupling constants which are presented in Fig. 3 and Fig. 4.
So we scan the parameter space to get the maximum value of the partial decay width. The results are given in Fig. 10. For the B c → B ( * ) χχ † channels, as there are no constraints for the effective coupling constants available now, so the upper limits of their branching ratios cannot be calculated.
One notices that, the Fig. 10(a) and Fig. 8(a) are almost exactly the same though the calculation processes are very different. For the later, χ is a scalar particle, and only the scalar coupling operator takes effect. For the former, χ is a pseudoscalar, and two effective operators make contribution to the branching ratios. When we scan the parameter space in Fig. 3, which are right triangular regions, we find that if the right endpoint of the hypotenuse is taken, the branching ratio will achieve the maximum. This also means only the scalar operator should be considered. So when we calculate the upper limit of the branching ratios, the operator contributes to the decay modes of Fig. 10(a) is just the same as that contributes to Fig. 8(a), which makes their results are the same. But the condition in Fig. 10(b) is quite different. One can see the upper limits of the branching ratios are larger than those in Fig. 8. When we scan the parameter space in Fig. 4, we find the axial vector operator provides most of the contribution, which is different with the case in Fig. 8(b), where only the pseudoscalar operator takes effect. For the B − c → D * − χχ † channel, there is a kink when m χ is larger than 2 GeV. This is because the corresponding operator which has the most important contribution turns to the vector from the axial vector.
In Fig. 11 we present the differential distribution of the upper limits of the widths as a function of s. As examples, two cases with m χ = 0 GeV and 0.25(M − M f ) are considered both for χ being a scalar or a pseudoscalar particle. In Fig.11(a) and Fig.11   is very small, there is no peak. For comparison, we also plot the differential widths for the s νν channels, which are smaller than those of the LDM channels in most regions of s. For the channels B c → D s νν, B c → D * s νν, and B c → B ( * ) νν, their decay widths are too small to be shown in Fig. 11(c), Fig. 11(d), and Fig. 11(e), respectively.