$B$-meson decay into a proton and dark antibaryon from QCD light-cone sum rules

The recently developed $B$-Mesogenesis scenario predicts decays of $B$ mesons into a baryon and hypothetical dark antibaryon $\Psi$. We suggest a method to calculate the amplitude of the simplest exclusive decay mode $B^+\to p \Psi$. Considering two models of $B$-Mesogenesis, we obtain the $B\to p$ hadronic matrix elements by applying QCD light-cone sum rules with the proton light-cone distribution amplitudes. We estimate the $B^+\to p \Psi$ decay width as a function of the mass and effective coupling of the dark antibaryon.

1. The B-Mesogenesis scenario suggested in [1,2] (see also [3,4]) simultaneously solves the problems of the baryon asymmetry and relic dark matter abundance in the Universe. The effective four-fermion interactions emerging in the B-Mesogenesis models generate decays of the b -quark into a light diquark and a dark antibaryon Ψ such that the baryon charge is conserved. These processes trigger new exotic decay modes of B mesons and Λ b baryons into final states containing light hadrons and Ψ, the latter revealing itself as an "invisible" i.e. a missing energy-momentum in the particle detector.
According to [1,2], the decays B → hadrons + Ψ should have appreciable branching fractions for a realization of the B-Mesogenesis scenario, with an inclusive width in the ballpark of 10 −4 . It is therefore very important to estimate the branching fractions of separate exclusive decay modes in order to investigate if these modes are within the reach of the ongoing searches for B meson decays to invisibles. In [2], only the ratios of exclusive and inclusive widths were estimated using a phase-space counting of quark states. Here, we consider the simplest two-body decay B + → p + Ψ with a single proton and missing energy in the final state. Since a hadronic B → p transition is involved, albeit an unusual one, we deal with the problem of calculating the hadronic matrix element of this transition within a QCD-based framework.
In this letter, we demonstrate how to solve this problem using the QCD light-cone sum rules. This method introduced in [5,6,7] and applied to various hadronic matrix elements can easily be extended to the meson-baryon matrix elements. To this end, we introduce a proton-to-vacuum correlator of the B-meson interpolating current with the effective three-quark operator coupled to the Ψ particle. We calculate this correlator near the light-cone in terms of nucleon distribution amplitudes (DAs). After that we derive the LCSRs for the B → p form factors, employing dispersion relations and quarkhadron duality in the B-meson current channel. Earlier, the LCSRs with nucleon DAs have been used to calculate the nucleon electromagnetic form factors [8,9,10] and the semileptonic form factors of heavy baryons (see e.g. [11]). On the other hand, the Bmeson interpolating current is widely used in the LCSRs for various B → h form factors, where h is a light meson. Hence, there are practically no new or unknown inputs in the new sum rules obtained below.
2. According to [1,2], the decay of a B + meson into a proton and dark antibaryon Ψ is generated by an exchange of a heavy colour-triplet scalar field Y . This model has two versions which differ by the electric charge Q Y of the heavy mediator. In this paper, for definiteness, we consider only the version with Q Y = −1/3, where the interactions of the Y -field with the quarks and Ψ are encoded in the Lagrangian: where the index c (R) indicates charge conjugated (right-handed) fields, q R = 1 2 (1+γ 5 )q, i, j, k are the colour indices and only the terms relevant for the B → pΨ decay are shown.
The effective four-fermion interactions generating this decay are given in [1,2] in a generic form. Here we derive more detailed expressions with a definite Dirac structure. The effective interactions are obtained by forming the Y exchange diagrams with the two vertices taken from the Lagrangian L (−1/3) and neglecting the momentum transfer m B with respect to a very large mass m Y in the propagator. For the model (1), we obtain the local effective Hamiltonian: We find it convenient to rewrite the above formula with the charge conjugation matrix attributed to the dark fermion field Ψ. To this end, we use the equations 1 which can be proven by employing the explicit form of the charge conjugation transformation for a Dirac field. We then factorize out the "external" Ψ field from the effective Hamiltonian: where G (d) = (y ub y Ψd )/M 2 Y is the effective four-fermion coupling and are the local three-quark operator and its conjugate. The index (d) indicates the terms in the effective Hamiltonian (4) originating from the d-quark coupling to ΨY . Replacing d ↔ b yields terms in (4) with the effective coupling G (b) = (y ud y Ψb )/M 2 Y and, respectively, with the three-quark operators We treat the two parts of H (−1/3) with the operators (5) and (6) as two sample models of the B-Mesogenesis, denoting them as model (d) and model (b), respectively. In [2] they correspond to the "type 2" and "type 1" operators.

3.
We consider first the model (d) in which the decay B → pΨ is generated by the effective interaction involving the three-quark operator O (d) . It is then straightforward to write down the decay amplitude where u c Ψ is the (charge conjugated) bispinor of the Ψ field. The four-momentum assignment is shown explicitly in (7) and the on-shell conditions are P 2 = m 2 p , (P + q) 2 = m 2 B , q 2 = m 2 Ψ . The hadronic matrix element of the three-quark operator in (7) can be parameterized in terms of four independent B → p form factors: where u p R,L = 1 2 (1 ± γ 5 )u p is the bispinor of the right-handed (left-handed) proton. Indeed, there are four independent kinematical structures constructed using the bispinor and two four-momenta P and q. Note that the structuresū p L,R (P ) / P are reduced tō u p L,R (P ) due to the Dirac equation.
Substituting (8) into (7) at fixed q 2 = m 2 ψ and using the Dirac equation for a chargeconjugated bispinor, / q u c Ψ (q) = m Ψ u c Ψ (q) , we obtain the decay amplitude in a compact form: Here, the coefficients in the scalar and pseudoscalar structures are linear combinations of the B → p form factors: In order to square the amplitude (9) and perform the sum over the proton and Ψ polarizations, we use the relation Multiplying the squared amplitude by the kinematical factors, we obtain the two-body decay width in terms of the form factor combinations defined in (10): where λ represents the Källen function. Turning to the model (b) with the operator O (b) , one has to replace d → b in the above equations (7)-(11).

4.
Here we derive the LCSRs for the B → p form factors starting from the model with the three-quark operator O (d) . The key element of our method is the correlation function It contains the time-ordered product of O (d) and the B-meson interpolating current j B = m bb iγ 5 u with the four-momentum P + q. The bilocal operator product in (12) is sandwiched between the on-shell proton state and the vacuum. The correlation function (12) is decomposed in kinematical structures: where Π R,L , Π R,L are Lorentz-invariant amplitudes.
If the external momenta P + q and q are far off-shell, so that (P + q) 2 m 2 b and q 2 m 2 b , the integration region in (12) effectively shrinks to a domain near the light-cone x 2 ∼ 0, where the propagating b-quark is highly virtual. In this region, we calculate the correlation function in terms of perturbative amplitudes convoluted with the light-cone distribution amplitudes (DAs) of the proton. Similar correlation functions were used to obtain the LCSRs for Λ b → p form factors [11].
Below we will employ the definitions of the nucleon three-quark DAs in which the proton is in the initial state, hence our choice for the initial and final state in (12). As a result, we will actually obtain LCSRs for the inverse hadronic transition p → B. Accordingly, (12) contains the three-quark operator O (d) from (5) instead of its Dirac conjugate in the decay amplitude (7). This choice is in fact inessential, because the p → B form factors that will be obtained from LCSRs coincide with the B → p form factors in the decay amplitude up to an irrelevant global phase.
To access the form factors from the correlation function (12), we employ the hadronic dispersion relation in the variable (P + q) 2 : where we isolate the B-meson pole. The contributions of the excited and continuum states located above the lowest threshold s h = (m B + 2m π ) 2 are collected in the integral over the spectral density ρ h(d) . Note that subtractions are neglected in (14) in anticipation of the Borel transformation. Equation (14) yields dispersion relations for the invariant amplitudes defined in (13). To separate them, we rewrite the proton-to-B hadronic matrix element in terms of form factors and use, analogously to (13), the decomposition into invariant functions for the hadronic spectral density ρ h(d) . We obtain for the first invariant amplitude where we use the definition of the B-meson decay constant: The remaining three dispersion relations are obtained from the above by replacing, respectively, Note that in (15) we tacitly replace the p → B form factors by the B → p ones, ignoring an irrelevant global phase. Below, we will calculate the correlation function (12), employing the light-cone OPE in terms of the proton DAs. The results for invariant amplitudes will be cast in a convenient form of dispersion relations, e.g. for the amplitude Π This representation allows us to employ the usual assumption of (semi-global) quarkhadron duality: introducing the effective threshold s B 0 . Substituting (17) into (15), we then subtract from both parts of the resulting equation the integrals that are equal due to the duality approximation (18). Performing the standard Borel transformation (P + q) 2 → M 2 , we finally arrive at the LCSR for the first form factor: The sum rules for the other three form factors are obtained from the above one by the replacements (16).
The main advantage of our method is its apparent universality. Turning to the model (b), we derive the LCSRs for the corresponding form factors, repeating the same steps as above. We only need to replace the three-quark operator O (d) by O (b) in the correlation function (12).

5.
To calculate the correlation function (12), we write down the operators in terms of quark fields: so that the Dirac indices inside the square brackets are contracted. The index T denotes a transposed field, and C is the charge conjugation matrix. We hereafter neglect the u, d-quark masses and assume isospin symmetry.
(P + q) Figure 1: Diagrammatic representation of the correlation function (12). The wavy (dashed) line represents the interpolating current of the B meson (the dark antibaryon).
In the deep spacelike region, (P + q) 2 , q 2 m 2 b , the expansion near the light-cone, x 2 = 0, is applied to the operator product in (20). To leading O(α 0 s ), the correlation function is described by the diagram shown in Fig. 1, in which the perturbative part is reduced to the free b-quark propagator and the long-distance part represents a protonvacuum matrix element of the three quark fields u(x), u(0), d(0). We replace this matrix element by a set of the proton light-cone DAs, equivalent to the nucleon DAs in the isospin symmetry limit.
Since our main goal is to present the LCSR method and to provide an approximate estimate of the B → p form factors, we retain only the lowest twist-3 DAs. The highertwist and less singular terms of the light-cone expansion in the proton-vacuum matrix element are neglected 2 .
We use the definition of the leading twist-3 DAs from [8]: where the variables α 1,2,3 are the fractions of the proton longitudinal momentum carried by the three quarks. The DAs V 1 , A 1 and T 1 have a polynomial form following from the collinear conformal symmetry. Adopting the usual next-to-leading approximation in the conformal spin expansion (see e.g. [8]), we use the following expressions: where the normalization parameter f N is the nucleon decay constant. The two additional parameters V d 1 and A u 1 characterize the deviation from the asymptotic form of DAs. The latter corresponds to V d 1 = 1/3, A u 1 = 0. Below, we will also need the once-integrated DAs: Substituting in (20) the b-quark propagator and rewriting the proton-to-vacuum matrix element via DAs according to (21), we integrate over x and over the virtual b-quark momentum, obtaining: The expression in square brackets is the invariant amplitude Π R ((P + q) 2 , q 2 ). We note that all other amplitudes from the decomposition (13) vanish in the twist-3 approximation.
To proceed, we transform the invariant amplitude to a form without (P + q) 2 in the numerator: The term independent of (P + q) 2 is omitted, because it vanishes after Borel transformation. The transition to a dispersion integral form is performed introducing a new integration variable such that the integral ds. Further steps include the quark-hadron duality approximation, which replaces the upper limit of the s-integration by an effective threshold s B 0 and the Borel transformation 1/(s − (P + q) 2 ) → e −s/M 2 . After that, we obtain the r.h.s. of (19), where it is convenient to return to the original integration variable α 1 ≡ α. Our main result, the LCSR for the form factor F (d) B→p R (q 2 ), reads: where s(α) is given by (27) and the effective threshold transforms to Following the same procedure as above, we derive the LCSRs for the B → p form factors in the model (b). We start, correspondingly, from the correlation function: and the OPE result in the twist-3 approximation reads: Note that in this case we encounter two different invariant amplitudes, hence two form factors, for which we obtain the following LCSRs:  Table 1: The input parameters in the LCSRs.
6. The input parameters for LCSRs are collected in Table 1. The b-quark mass is taken in M S-scheme. The adopted renormalization scale and its variation are optimal for a correlation function with the B-meson interpolating current, following e.g. the recent analyses of LCSRs for the B → π form factors and B * Bπ strong couplings [14,15]. Accordingly, we use the same range for the Borel parameter. We expect this range to be large enough for the validity of the leading twist approximation in LCSRs, since the yet unaccounted higher twist DAs lead to contributions that are typically suppressed by inverse powers of M 2 . We also checked that, with the same choice of the Borel parameter range, the contributions of excited and continuum states to LCSRs do not exceed 10% (20%) for the model (d) (model (b)). Hence, our results are not too sensitive to the quark-hadron duality approximation. The effective threshold s B 0 is then fixed for each value of M 2 , fitting to the B meson mass the LCSR differentiated with respect to −1/M 2 . Furthermore, for the B-meson decay constant we use the average [16] of the N f = 2 + 1 + 1 lattice QCD determinations.
For the input parameters determining the nucleon DAs, we take the most recent lattice QCD results [17] at the scale µ 0 = 2.0 GeV. The scale dependence is calculated using the known one-loop anomalous dimensions, e.g. for the nucleon decay constant: where β 0 = 11 − 2/3n f and n f = 4 for the scales of our choice. The parameters A u 1 and V d 1 are equal to linear combinations of the multiplicatively renormalizable nonasymptotic coefficients ϕ 10 and ϕ 11 (for details see e.g., [10]). We notice that the contributions of ϕ 10 cancel in the linear combination of A u 1 and V d 1 entering both integrated DAs in (23) and (24). As a result, the LCSRs in the leading twist depend on a single nonasymptotic coefficient: equal to 0.148 ± 0.03 at µ = 3 GeV, if we use ϕ 11 (µ 0 = 2 GeV) from Table 1.
With the input specified above, the B → p form factors are then calculated from the LCSRs (28), (31) and (32) at q 2 m 2 b , including also the spacelike region q 2 < 0. Hence, we can directly obtain the B → pΨ amplitude at a dark-antibaryon mass, m Ψ = q 2 , in the lower part of the expected [1,2] 4.34 GeV. To extrapolate the form factors to the upper part of that range, we employ the zexpansion, mapping the complex q 2 -plane onto the unit disk on the plane of the variable More specifically, we use the BCL version [18] of z-expansion, slightly modified according to [14]. For instance, we use for the model-(d) form factor: Note that the pole factor in (35) reflects the pole of Λ b -baryon located below theB + p threshold. Similar expansions parameterized by the value at q 2 = 0 and by the slope parameter are used for the form factors in the model (b). For the fit of the form factors from LCSRs to their z-expansion, we have chosen the interval −5.0 < q 2 < +1.0 GeV 2 (0.113 ≥ z ≥ 0.077). The adopted order of this expansion is quite sufficient, since the extrapolation interval 1.0 GeV 2 < q 2 < (m B −m p ) 2 also maps onto small z-values, 0.077 > z > −0.083. The fitted parameters of (35) and of the analogous expressions for the two form factors in the model (b) are presented in Table 2.  Table 2: Parameters of the z-expansion for the B → p form factors (in GeV 2 units).
The quoted uncertainties are obtained by varying the input parameters within the adopted intervals and adding them in quadrature. We observe a significant suppression of the form factors in the model (b) with respect to the model (d). This effect can be traced to the O(m 2 p /m 2 b ) suppression of the LCSRs (31) and (32) with respect to (28). Without going into further details, we conclude that there is a nontrivial sensitivity of the hadronic B → p form factors to the configuration of quark fields in the effective operator.
Finally, substituting the form factors in the decay amplitude (9) and its analog for the model (b), we calculate the B → pΨ decay widths in both models: where in the last expression we take into account that the form factors are real valued in the region below (m B + m p ) 2 . Multiplying the above expressions by the B ± lifetime, τ B ± = 1.638 ± 0.004 ps [13], we arrive at our main results, the branching fractions of the B + → pΨ decay. They are plotted in Fig. 2. For uniformity, we assume equal effective four-fermion couplings |G (d) | 2 = |G (b) | 2 = 10 −13 GeV −4 . This choice is in the ballpark of the upper limits extracted [2] from the LHC constraints on the heavy coloured scalars Y predicted in the B-Mesogenesis models. Our predictions for the branching fractions can be easily rescaled for any other value of these couplings.

7.
In this letter, we presented the first estimate of the B-meson decay width into a proton and dark antibaryon Ψ emerging in the B-Mesogenesis framework. Our main results are the B → p form factors determining the hadronic amplitude of this decay. They were calculated as a function of Ψ mass from the QCD-based LCSRs with the proton distribution amplitudes. Our results reveal a strong dependence of these form factors on the form of the effective operator originating from the three-quark coupling with a dark antibaryon.
The accuracy of our calculation is limited by the lowest twist of the nucleon DAs and by the leading O(α 0 s ) of the underlying OPE. Both higher-twist and O(α s ) contributions can be systematically calculated in the future. However, we do not expect a significant impact of these contributions on the numerical results obtained here.
It is straightforward to apply the LCSR method suggested here to other versions of the B-Mesogenesis scenario, e.g. to a superposition of the models (d) and (b), or to the models [2] with the hypercharge Q Y = +2/3 of the heavy mediator. Other exclusive B-decays, for example, into a strange baryon and Ψ, are also accessible with a certain modification of the LCSRs. A complete assessment of the B → hadrons + Ψ decays is indispensable without reliable estimates of their inclusive widths. Here, one has to employ quite different methods, e.g. the well established heavy-quark expansion.
Finally, our results can be useful for ongoing and future experimental searches for the B-decays to invisibles. According to [2], the Belle II experiment is sensitive to branching fractions in the ballpark of BR(B → pΨ) = 3×10 −6 . Comparing with our results plotted in Fig. 2, we conclude that, at least for the model (d), such a sensitivity is well within the parameter space formed by the limits on the effective coupling G (d) and the mass m Ψ , provided this mass is fixed by measuring the proton energy in the B meson rest frame. This enables a realistic test of the B-Mesogenesis scenario in the future.