Higher twist corrections to $B$-meson decays into a proton and dark antibaryon from QCD light-cone sum rules

The $B$-Mesogenesis framework anticipates decays of $B$ mesons into a dark antibaryon $\Psi$ and various Standard Model baryons. Here, we focus on the exclusive decay process $B\to p \Psi$ observed as a proton and missing energy in the final state and determine the decay width by employing the QCD light-cone sum rule framework. We include all contributions up to twist six to the nucleon distribution amplitudes in order to parameterize the non-perturbative effects in the operator product expansion. We obtain the decay width and branching fraction with respect to the mass $m_{\Psi}$ of the dark antibaryon $\Psi$, normalized to the model-dependent effective four-fermion coupling.


Introduction
Although the Standard Model of particle physics (SM) is up-to-date the best established theoretical framework to describe particle interactions, it lacks explanations for many observed phenomena like the baryon asymmetry or the dark matter abundance of the universe.Hints on the matter-antimatter asymmetry of the universe can for instance be deduced from measurements of the Cosmic Microwave Background (CMB) [1,2] or the Big Bang Nucleosynthesis (BBN) [3,4].In general, there exist many different theoretical approaches in the literature which try to address these problems from a theoretical point of view.However, the typical scales involved in these scenarios lie around the Planck scale and are therefore hard to verify experimentally.The B-Mesogenesis model proposed in [5][6][7][8] has emerged as an elegant solution to these puzzles as its features become apparent at measurable energy scales.Previous studies suggest that the decay B → hadrons + Ψ is expected to possess appreciable branching fractions with an inclusive width of the order of 10 −4 [5,6,8].To explore the feasibility of detecting these modes, a deeper investigation of separate exclusive decay channels becomes important.The original work in [6] roughly estimated the ratios of exclusive to inclusive widths utilizing phase-space counting of quark states.But as it has been shown in [9] for the two-particle decay B → pΨ and later for different decay channels in [10], the QCD light-cone sum rule (LCSR) approach is well suited to provide estimates for these exclusive decays.LCSRs were initially introduced in [11][12][13] and subsequently applied to various hadronic matrix elements.In [9], the hadronic B → p transition posed the central challenge in the computation.In order to address this problem, the authors have rather investigated the p → B transition as it only differs by a global phase to the desired B → p decay.The advantage of inverting this transition lies in the fact that they need to employ nucleon distribution amplitudes to parameterize the non-perturbative contributions in the LCSR approach, which are studied in greater detail in [14][15][16][17][18] than B-meson distribution amplitudes.In addition to that similar computations have already been carried out in the calculation of form factors for the Λ b → p transition [19], which are also applicable for the desired B → p transition.However, only the leading-twist contributions have been considered in order to obtain a first estimate for the corresponding branching fractions.Following the approach from [9], we also focus on the specific two-body decay B + → p + Ψ, but include all contributions to the nucleon distribution amplitudes up to twist six.This allows us to perform a dedicated study on the reliability of the leading-twist contributions and the impact of higher twist corrections on the decay width and branching fractions.Moreover, we obtain an estimate on the convergence of the OPE itself.This works is organized as follows: In section 2, we introduce the formalism including the basic features of the B-Mesogenesis model relevant for the B → p transition.Therein, we state the input parameters of the model as well as the effective four-fermion operators with the new dark matter particle Ψ.Furthermore, section 3 is devoted to the derivation of the LCRSs, while section 4 illustrates the computation of the various OPE contributions.This leads to the expressions for the various form factors, which we extrapolate to the physical timelike region in section 5. Section 6 introduces the parameters of the nucleon distribution amplitude and their dependence on the renormalization scale µ and subsequently shows the numerical evaluation of these form factor expressions in subsection 6.1 and of the branching fractions in subsection 6.2.Finally, we conclude in section 7. Appendix A and B provide further supplementary information.

Effective operators
As we have already pointed out, the decay B → pΨ is one of the simplest decay channels which lead to the introduction of a dark sector restoring baryon number conservation by considering the combination of the SM and the dark sector.In order for the dark matter particle Ψ to be observable and not to decay into SM particles immediately, it is only allowed to interact gravitationally with the SM or via a heavy colour-triplet scalar field Y with a mass of the order of a few TeV.However, the hypercharge of the scalar mediator particle Y is not unique, there exist in general two different possibilities with Q Y = −1/3 and Q Y = 2/3.Throughout this work, we only focus the model with Q Y = −1/3 for simplicity, because these considerations can be similarly applied to the second case.The part of the Lagrangian governing the additional interactions of the Y -field with SM quarks and the dark matter field Ψ is given by ) with c(R) denoting charge conjugated (right-handed) fields, q R = 1 2 (1 + γ 5 )q, while i, j, k indicate the colour indices of the quarks in fundamental representation.Notice that we expect a completely antisymmetric combination of colour charged fields in the interaction of Y with SM quarks in order to preserve gauge invariance [5].Additionally, the quantities y ud , y Ψd and y ub , y Ψb represent the (antisymmetric) Yukawa couplings between the dark sector and the SM sector.Similar to [9], we exploit that the mass of the interaction particle Y , M Y , is much greater compared to the typical momentum transfers of this decay k ∼ m B and therefore integrate out the heavy mediator Y .Effectively, the propagator turns into such that we obtain the following effective Lagrangian including four-fermion interactions We depict the corresponding Feynman diagram for the effective interaction in figure 1.
From the expression in Eq. (3), we can immediately read off the effective Hamiltonian Figure 1: Diagram for the p → B transition taken from [9], which differs from the B → p transition describing the B → pΨ decay by an unobservable global phase.The mediator particle Y has been integrated out such that we obtain an effective four-fermion interaction containing the new dark matter particle Ψ.
For the extraction of the effective three-quark operators from the four-fermion interaction in Eq. ( 4), it is useful to employ the Fierz identity [20] Ψd After that we factorize the field Ψ from the effective four-fermion interaction and obtain with the effective four-fermion coupling G (d) = (y ub y Ψd )/M 2 Y and define the local threequark operator and its conjugate It is also possible for the b-quark to couple to the dark matter particle Ψ, which leads to the following operators: In the remaining analysis, we consider the operators in Eqs. ( 7) and ( 8) as two individual versions of the B-Mesogenesis model and call them (d)-and (b)-model, respectively.This is in analogy to [6], where these operators are referred to as "type-II" and "type-I" operators.These operators in Eq. ( 7), ( 8) constitute the central elements of our framework, since the correlation function includes these operators for the determination of the B → pΨ decay.
The decay amplitude for this particular decay in the (d)-version of the B-Mesogenesis model is given by and similarly for the (b)-model if we replace the operator in Eq. ( 9) by the corresponding operator from Eq. ( 8).We choose the momentum assignment according to figure 1 such that the on-shell conditions read (P + q) 2 = m 2 B , p 2 = m 2 p and q 2 = m 2 Ψ .Furthermore, we decompose the B → p transition into the four different form factors Note that we introduce the factor 1/m p for the last two Dirac structures to render the form factors dimensionless.After replacing (d) → (b), we obtain the corresponding form factors for the (b)-model.These set of form factors will be determined via the light-cone sum rule approach in the following sections.

Derivation of the light-cone sum rules
The correlation function plays the key role in the derivation of the light-cone sum rules, since it directly connects the physical timelike region with the perturbatively calculable spacelike region via the quark-hadron duality (QHD).Therefore, we begin this discussion by stating the form of the correlation function which corresponds to the effective framework introduced in the last section 2 where j B (x) = im b b(x)γ 5 u(x) is the B-meson current.For the second version of the Bmesogenesis model, the operator O (d) needs to be replaced by O (b) from Eq. ( 8).Contrary to the discussion above, we investigate the p → B transition rather than the necessary B → p transition for the decay B → pΨ.These transitions differ at most by a global phase which do not alter physical observables like decay widths or branching fractions.One advantage is that we can parameterize the long-distance contributions in terms of nucleon light-cone distribution amplitudes, which are better known than the B-meson DAs [14][15][16][17][18]21]. Especially the parameters of the nucleon DAs have been determined to better accuracy by advanced lattice computations and sum rules analyses such that we can perform our analysis to twist six accuracy.In this context, the second advantage is that we can closely follow the analysis for Λ b → p form factors from [19].Although the currents inside the correlation function differ for this problem, the computation of the form factors requires to use the same distribution amplitudes.We state the decomposition of the nucleon matrix element, its transformation into distribution amplitudes of definite twist and the relevant shape parameters in appendix A. The next step is to make use of the unitarity condition and the Schwartz reflection principle.While the former corresponds to the insertion of a complete set of states inside Eq. ( 11), the latter expresses this correlation function in terms of a dispersion relation in (P +q) 2 .As it is usually the case for transitions involving B-mesons, we can easily separate the ground state contribution in form of a B-meson pole.Beyond the threshold cutoff s h = (m B +2m π ) 2 , we observe hadronic contributions like excited states and continuum contributions, which we incorporate into the hadronic spectral density ρ h (d) .Employing in the hadronic dispersion relation in (P + q) 2 , we can access the form factors via Notice that we ignore possible subtraction terms in this context as they vanish during a Borel transformation, which additionally suppresses the continuum contributions and leads to a better convergence of the sum rules.The correlation function in Eq. ( 11) consists of different kinematical contributions, which are given by all possible Lorentz-invariant amplitudes: Here, we use that the Dirac spinors u p,{R,L} (P ) are the right-handed and left-handed components of the Dirac spinor u p (P ), i.e. u p,{R,L} (P ) = P {R,L} u p (P ) = 1 2 (1 ± γ 5 )u p (P ).Furthermore, contributions involving structures like / P u p,{R,L} can be included in the first line of Eq. ( 13) due to the Dirac equation.Based on the decomposition into different kinematical contributions, we can derive individual dispersion relations for the various form factors in Eq. (10).For this, we insert (13) into Eq.( 12) and group contributions according to their different Dirac structures such that we can directly access the individual form factors.For example, the dispersion relation for the form factor where we exploit that In addition to that we decompose the hadronic spectral density ρ h(d) in terms of the decomposition in Eq. (13).The other form factors can be obtains similarly by the following replacements: Note that we introduce a factor of 1/m p in front of the form factors B→p R,L (q 2 ) to obtain the correct mass dimensions.Before heading to the calculation of the correlation function, which will be covered in the next section, we need to find a proper expression for the hadronic spectral density.All results from the next section can be expressed via a dispersion relation of the OPE contributions, which we compute in the deep spacelike region by employing perturbative methods.For the amplitude Π where we write ρ . In order to remove the hadronic spectral density, which is hard to describe from a theoretical point of view due to its complicated structure, we replace the integral over the hadronic spectral density ρ h(d) by an integral over the spectral density obtained from the OPE 1 π ImΠ (d),OPE using the (semi-local) quark-hadron duality This step introduces the effective threshold s B 0 , which needs to be determined in the numerical analysis as it is an input parameter in the sum rule framework.Finally, we substitute Eq. ( 16) into Eq.( 14) and perform a Borel transformation in the variable (P + q) 2 to arrive at the desired form of the sum rules: As we have previously discussed, we obtain the other form factors with the replacements from Eq. ( 15).

Correlation Function
Now that we derived the form factors in terms of the sum rules corresponding to Eq. ( 18), we can evaluate the OPE to the desired twist accuracy and to leading order in α s .
In order to apply light-cone sum rules, we need to work in the phase space region where (P + q) 2 ≪ m 2 b and q 2 ≪ m 2 b is valid, which means that the momenta of the involved particles are far off-shell.According to figure 1, the b-quarks inside the correlation function in Eq. ( 19) are connected to form a b-quark propagator.This leaves us with a proton-to-vacuum matrix element containing two uncontracted u-quark fields as well as a d-quark, which encodes the non-perturbative information in the p → B transition.Following the usual procedure of the light-cone sum rule approach, we replace this matrix element by nucleon DAs by decomposing the matrix element in terms of different Lorentz structures based on Lorentz invariance and parity invariance.We obtain in total 24 structures if we consider all contributions up to twist six accuracy, which can be subsequently related to distribution amplitudes of definite twist.We provide all details regarding this procedure and the further parameterization of the shape of the distribution amplitudes in the conformal expansion in appendix A. Contrary to previous investigations [9], see also [10], we consider all contributions up to twist six accuracy including the O(x 2 ) corrections to the leading-twist contributions.O(x 2 ) corrections to twist four contributions, which correspond to a twist six effect, are numerically negligible [15] and not considered in this work.By including these set of contributions, we are able to estimate the reliability of the leading-twist analyses from [9] since we can study the convergence of the OPE and observe the impact of higher twist corrections on the branching fraction of the B → pΨ decay itself.We start the computation by reproducing the known leading-twist results from [9]: In general, the nucleon DAs introduce three different variables α 1,2,3 representing the momentum fractions of the individual quarks inside the proton.In Eqs. ( 20) and ( 21), we have integrated over α 2 and α 3 and renamed α = α 1 .Moreover, we observe that for the (d)-model only one form factor contributes, while there is an additional contribution for the (b)-model from the Lorentz structure / qu p,L (P ).We state the leading-twist distribution amplitudes V 1 , A 1 , T 1 in appendix A, where we also elaborate on the other higher twist distribution amplitudes in more detail.Scalar products of the form P • q are not suitable for Borel transformations, hence we replace them by This allows us to cancel the (P + q) 2 -dependence in the denominators of Eqs. ( 20) and ( 21) after rewriting twist 3 d d x e ikx e i(P +q)x e −iP αx = (2π) d δ (d) (k + q + P ᾱ) twist 4, 5, 6 Table 1: Integration over the position variable x with additional factors of x µ starting at twist four accuracy.For brevity, we introduce the notation ᾱ = 1 − α.
Constant terms independent of (P +q) 2 vanish under the subsequent Borel transformation.Higher twist corrections become more involved, because they explicitly show x-dependencies in the expressions.First of all, they occur as explicit factors x µ in our calculation, which we rewrite according to table 1 as derivatives acting on the momentum-conserving δdistribution.Besides that we deal with scalar products of the form P • x, which need to be introduced in order to relate the 24 invariant functions S i , P i , A i , V i , T i in the decomposition in Eq. (42) to distribution amplitudes of definite twist, see appendix A, and for instance [14,15] for more details.An additional partial integration with respect to the variable α removes these factors and we additionally notice that the occurring surface terms vanish.With these steps in mind, we can perform a similar derivation of the form factors as in [9].We intend to bring the contributions from Eqs. (20) and (21), including now all contributions up to twist six, into the form of an dispersion integral.For this, we perform the following substitution after using Eq. ( 23) Finally, we need to perform the Borel transformation to obtain the final form of the LCSRs for the form factors where B M 2 denotes the Borel transformation.In our case we identify Q 2 = −(P + q) 2 and we encounter the cases k = 1, 2, 3.This is related to the fact that the denominators in Eqs. ( 20) and ( 21) occur to higher powers if we consider the higher twist corrections.Ultimately, this leads to a suppression of these factors in powers of 1/M 2 .Performing all these steps for instance for the F (d) B→p R (q 2 ) form factor results into We state the relevant notation and various functions in appendix A. Notice that we introduce the notation to denote the distribution amplitudes which we integrate in α once or twice in order to remove the scalar product P •x.With the replacements from Eq. ( 15), we can derive similar expressions for the other form factors.It turns out that the (d)-model receives additional contributions for the B→p L (q 2 ) form factor such that we end up with two different form factors in each model.However, we note that the T -structures in the proton matrix element in the decomposition (42) result into sizeable effects for the (b)-model, while they vanish in the (d)-model.We provide the remaining form factor in the (d)-model and the form factors for the (b)-model in appendix B.

Extrapolation to the large m Ψ -region
In the last section, we have obtained the form factors for both the (d)-and (b)-model, which ultimately enter the decay width for the process B → pΨ.These form factors are valid in the limit q 2 ≪ m 2 b since we have used the light-cone sum rule approach, meaning that we are still working on the light-cone with small distances.However, the kinematics of this two-particle decay shows that the upper bound of the particle mass m Ψ is given by m B − m p ≈ 4.34 GeV, while the lower bound lies around m p in order to prevent proton decays [5,6,8].Hence, the mass of the dark matter particle Ψ is theoretically allowed to be in the range m Ψ ∼ m b .
In order to extract reliable results from the form factors expressions derived in the last section, we perform an extrapolation of Eqs. ( 26) and ( 55) to (57) using the BCL-version [22] of the z-expansion [23].For this, we perform a conformal mapping of the variable q 2 onto the complex variable z: with While t 0 is a default parameter, usually chosen as above, in order to minimize the truncation error of the z-expansion, t ± are set by the physics of the decay.The parameter t − = m B − m p is precisely the upper bound on m Ψ dictated by the two-particle decay of the B meson at rest and t + = m B + m p constitutes the threshold for multiparticle states and higher resonances.Starting from this threshold, the timelike form factors become imaginary, which is also represented by the variable z developing an imaginary part.However, the choice of t + introduces another subtlety, namely an isolated pole due to the Λ b -baryon at q 2 = m 2 Λ b .The treatment of these particular issues is further described in [22] such that we finally end up with Here, it is sufficient to truncate the z expansion to O(z 2 ), because the allowed range for the mass m Ψ results in small values of z in the interval 0.077 > z > −0.083.This leaves us with two free parameters which we need to determine further.These two parameters are given by the form factor evaluated at q 2 = 0, which we can determine directly from Eqs. ( 26) and ( 55) to (57), and the slope parameter b B→p R , which we get from the fitting procedure in the next section.For the other form factors, one just has to perform the replacement rule in Eq. ( 15) and additionally (d) → (b).

Numerical analysis
The input parameters for the LCSRs are given in table 2. We perform our analysis at a renormalization scale of µ = 3 GeV.This particular choice is in accordance with recent studies on the B → π transition or B * Bπ strong couplings [25,26], because these analyses indicate that this scale and its uncertainty is optimally suited for B-meson interpolating currents.Additionally, we adopt the b-quark mass in the MS-scheme and use the B-meson decay constant obtained from recent lattice QCD computations with n f = 2 + 1 + 1 [27].−0.312 [16] Table 2: Input parameters in the LCSRs from the references in this Table.
As table 2 shows, the input parameters for the nucleon distribution amplitude are extracted from different sources.For example, a recent lattice calculation [21] determines some input parameters at the scale µ 0 = 2 GeV, whereas a LCSR computation [18] works at µ 0 = √ 2 GeV.Therefore, we make use of the following RGE to run these parameters to required scale µ = 3 GeV with γ φ being the non-cusp anomalous dimension for the DA parameter φ.Its solution to one-loop order is given by where γ 0 φ denotes the one-loop anomalous dimension.The value for Λ QCD = 0.288 GeV is taken from [28] for n f = 4. Analogously, Eq. (32) can be used in order to run the remaining nucleon DA parameters to the desired scale µ = 3 GeV.This requires the one-loop non-cusp anomalous dimensions γ 0 φ for all parameters, which we give in table 3  These parameters can be related to the parameters of the conformal expansion via [18]

Form factors
Next we determine the uncertainties around the central values of the form factors.This requires to consider the form factors before their z-expansion and we individually introduce variations to each input parameter within its predefined range of uncertainty.Subsequently, we perform a z-expansion for each parameter variation yielding two sets of parameters for the slope parameter b B→p R and the normalization F B→p R (0), namely the upper and lower bound on these parameters with respect to the input parameter variation.The values in each set are then added in quadrature to obtain the possible parameter range for both fitting parameters.Apart from the usual correlation between the threshold parameter s B 0 and the Borel parameter M 2 and the impact of the µ variation on the renormalization scale dependent parameters, we assume that the remaining parameters are completely uncorrelated.We state the slope parameters and the form factor at q 2 = 0 within their uncertainties in table 4 for the (d)-model and in table 5 for the (b)-model:   The right panel depicts the individual twist contributions to the different form factors such that a direct comparison between the higher twist corrections and the leading twist-three contribution becomes possible.These leading contributions have been previously examined in [9] and can also be found in [10].In this context, it is notable that the leading contribution to the form factor B→p L (q 2 ) starts at twist four accuracy, contrary to the other three form factors. Nevertheless, the OPE shows in general good convergence for all form factors, which is in accordance with other sum rule analyses including B-meson interpolating currents.
Although the OPE shows good convergence for all form factors, we observe at the benchmark value m Ψ = 2 GeV that the (b)-model exhibits a substantial twist-four correction, violating the typical hierarchy of the OPE since the twist-four contribution is larger compared to the leading-twist contribution.This has its origin in the significant T 2,4 contributions in the (b)-model, which vanish in the (d)-model computation.Nevertheless, the convergence of the OPE is still well established beyond twist four in the (b)-model.
In the left panel, we investigate the complete form factor expressions including all contributions up to twist six at the benchmark value m Ψ = 2 GeV.At this stage, it is interesting to compare these expressions for different choices of the threshold parameter s B 0 based on the upper and lower bounds specified in table 2. Given the small fluctuations of the form factors with respect to different Borel parameters across varying thresholds, we see that the sum rules for the four form factors remain stable at the reference value m Ψ = 2 GeV, thus ensuring their reliability.The choice of the Borel parameter range M 2 aligns with the specifications from table 2 and is in agreement with the findings in [9].Furthermore, an alternative validation of this Borel window can be performed by considering the form factors over a broader range of the Borel parameter M 2 and determining the Borel window based on the stability of the sum rule with respect to M 2 .As illustrated in figures 2 to 5, it is evident that the sum rules remain stable within the specific Borel window defined in table 2. Additionally, this choice of the Borel parameter window ensures that continuum and excited states are reasonably suppressed and contribute approximately around 20% to 30%.Consequently, our findings are not excessively influenced by the Quark-Hadron Duality (QHD) approximation.Moreover, the determination of the threshold parameter s B 0 is accomplished by calculating the derivative of the sum rules presented in Eq. ( 18) with respect to −1/M 2 .By taking the ratio between the derivative outcome with Eq. ( 18), we obtain an approximation for the B-meson mass m B , which is subsequently adjusted by choosing the value of s B 0 such B→p R (q 2 ) in the (b)-model for different Borel parameters M 2 .The left plot shows all twist contributions combined for various choices of the threshold parameter s B 0 , while the right plot illustrates each twist contribution individually.The mass of the dark matter particle m Ψ is set to the benchmark value m Ψ = 2 GeV [6].
that it fits existing literature data [4].The stability of the sum rules in figures 2 to 5 shows again the validity of this approach.B→p L (q 2 ).B→p L ;i defined in Eq. ( 33) and ( 34) with respect to the dark fermion mass m Ψ .The right panel contains the ratios for the form factor B→p L , while The left panel shows the contributions from the form factor So far, our discussion has been centered around the chosen value m Ψ = 2 GeV.Since the LCSR approach works in the limit q 2 = m 2 Ψ ≪ m 2 b , we naturally assume that the OPE also converges when m Ψ remains below 2 GeV.It is important to understand the applicability of the OPE across higher values of m Ψ , especially at which point the OPE starts to break down.This analysis is shown for the (d)-model in figure 6, focusing on the form factor B→p R (q 2 ) in the left panel and on the form factor F (d) B→p L (q 2 ) in the right panel.For these investigation, it is useful to introduce the ratios: where B→p L ;i belong to the twist i contribution of the form factors B→p L .Using (33) and (34), we derive an approximation for m Ψ where the OPE diverges.The underlying notion here is that when the higher twist contributions grow, which is indicated by increased ratios R B→p L ;i , these contributions become dominant and spoil the convergence of the OPE in general.Our observations reveal that for the form factor B→p R (q 2 ) in the left panel of figure 6, the convergence is spoiled around m Ψ ≈ 3 GeV.In comparison to that the form factor F (d) B→p L (q 2 ) turns out to be insensitive to this investigation.Nonetheless, for the earlier considered benchmark value of 2 GeV, the convergence  of the OPE remains robust.It is important to note that probing m Ψ beyond 6.2 GeV is unfeasible since we cross the multihadron threshold at t + , which results into complex z-parameter values.A similar examination can be carried out for the (b)-operator and this is shown in figure 7.This underlines that the hierarchy in the twist expansion remains intact for both form factors confirming that the expansion retains its convergence throughout the entire kinematical range of m Ψ .

Branching fractions
Before calculating the branching fraction, we need to obtain the decay amplitude of the considered B + → pΨ decay.Following [9], we start with the amplitude A (d) (B + → pΨ) from Eq. ( 9) and insert the decomposition into form factors for the B → p transition from Eq. ( 10) with We have expressed the decay amplitude through the four different form factors and the ratios of the proton mass m p as well as the dark matter particle mass m Ψ .In order to obtain the two-body decay width, one has to square the amplitude from Eq. ( 35) and multiply this expression with normalisation factors The phase space integration for a two-body decay can be carried out analytically yielding the Källen function The final expression for the decay width in terms of form factors is given by In comparison to the corresponding expression in [9,10], the form factor B→p L (q 2 ) contributes here due to higher twist contributions.Since this form factor already appears in the (b)-model at the leading-twist approximation, a simple exchange of (d) with (b) allows us to derive corresponding relations for the second model.The observable of interest is the branching fraction for this decay.This parameter can be derived by dividing Eq. ( 40) by the total width of the B meson.Alternatively, we can also multiply Eq. ( 40) by the B-meson lifetime τ B ± = 1.638 ± 0.004 ps from [24]: As shown in figure 8, we present the branching fraction for the (d)-model incorporating contributions to the nucleon distribution amplitudes up to twist six and compare them to the leading-twist contribution from [9,10].Within m Ψ values up to 3 GeV, we identify that both computations agree very well within their uncertainties.This observation aligns with our earlier observation in figure 6 that at m Ψ = 3 GeV the higher twist corrections become dominant and ultimately affect the convergence of the OPE.Generally, the uncertainties on our twist six calculation turn out to be larger compared to the previous leading-twist evaluation.This disparity arises due to the larger error estimates on input parameters of the distribution amplitudes in the conformal and nextto-conformal expansion.In particular, the parameter ξ 10 introduces large uncertainties as we assume a conservative error of 50% based on the value from [16].Nonetheless, we affirm in general that the leading-twist outcome constitutes a reliable approximation for the branching fractions within the m Ψ -range up to 3 GeV.However, a significant discrepancy between the twist-three calculation from [9,10] and our computation arises for the (b)-model.In this case, the branching fraction increases by roughly a factor of 20 and the two computations do not agree within their uncertainties.This deviation can be attributed to the substantial T 2,4 -contributions at the twist four level which we observe for both form factors, as it is evident from figures 4 and 5. Similar to the (d)-model, the uncertainties on the twist six computation are notable and share the same origin as for the (d)-model.Particularly, the upper bound uncertainty becomes enhanced once m Ψ > 3 GeV.This behavior illustrates that the branching fraction loses its reliability as m 2 Ψ ∼ m 2 b , thereby violating a crucial requirement for the light-cone expansion.The blue line with the blue dashed error band illustrates the original twist-three computation from [9], while the black curve shows the computation including contributions up to twist six.The dashed red curves represents the uncertainty on this calculation.Consequently, we conclude that the leading-twist analysis from [9] for this specific model falls short as higher twist corrections shift the value of the branching fractions considerably and spoil the hierarchy of the OPE.However, the branching fraction now lies with these additional contributions in the sensitivity range of Belle-II, estimated to be around 3 • 10 −6 [6].

Conclusion
The SM plays a crucial role in modern particle physics, although there are effects like dark matter or the baryon asymmetry of the universe which are not incorporated into this theory.As there are strong experimental hints for the existence of dark matter, several models are proposed of which we have studied the recently proposed B-Mesogenesis model [5,6,8].It is especially interesting for experimental facilities like Belle-II since it predicts new dark matter particles at energy scales which are in principle within its sensitivity range.
We have focused here on one allowed decay mode within this model, namely the decay B → pΨ.For this, we have used the light-cone sum rule approach to determine the branching fraction of the decay B → pΨ.While the leading-twist contributions have been obtained in previous works [9,10], we include higher twist corrections up to twist six into our analysis in order to check the reliability of the leading-twist analysis and the convergence of the operator product expansion in general.Thus, we have followed the procedure from [9] and computed the branching fractions of the two considered versions (the (d)and the (b)-model) in the B-Mesogenesis scenario, normalized to the effective four-fermion couplings G (d) and G (b) .We observe for the two form factors of the (d)-model, which contribute to the respective branching fraction, that higher twist contributions become increasingly smaller in the parameter range 0.94 GeV ≤ m Ψ ≤ 3 GeV, indicating that the OPE converges and that the leading-twist contribution is dominant.This shows that the higher twist corrections have a minimal impact on the branching ratio in this parameter range and that the behavior mostly follows [9].Beyond 3 GeV, we see that the that higher twist corrections start to dominate and therefore that the OPE breaks down.This behavior is expected since the light-cone approach requires that m 2 Ψ ≪ m 2 b , which starts to get violated beyond 3 GeV.Hence, we conclude that the branching fraction estimate from [9] is reliable in the regime m Ψ ≤ 3 GeV and in agreement with the results from this work.In contrast to that the twist four contributions tend to be the dominant contribution in the (b)-model, albeit the OPE converges beyond this twist correction.This has the consequence that the branching fraction for this model is significantly increased, namely by factor of 20.This observation has a direct impact on experimental searches, since both versions of the B-Mesogenesis model are now in the sensitivity range of Belle-II.From an experimental point of view, the decays B → ΛΨ → pπΨ or B → ΛΨ → pπΨ are however easier accessible as they see two SM particles in the final state combined with the new-physics particle Ψ in form of missing energy in the detector.But the necessary distribution amplitudes for the ∆-and Λ-baryons are less known and only a leading-twist analysis is available in the literature.Nevertheless, ratios of these decays with the computed branching fraction B → pΨ in this work would be independent of the couplings G (d) and G (b) and hence reduce the number of input parameters.

A Nucleon distribution amplitude
Generally, the hadronic matrix element from the previous discussion can be decomposed into different Lorentz structures based on symmetry considerations based on Lorentz covariance, spin and parity of the proton [14]: On the left side of Eq. (42), we investigate a proton to vacuum matrix element with an onshell proton of P 2 = m 2 p in the initial state.The quark fields u(x), u(0) and d(0) correspond to the valence quarks inside the proton.The Greek letters α, β, γ denote spinor indices, while Latin letters i, j, k are colour indices.Gauge link factors between each valence quark are suppressed rendering the expression in Eq. (42) gauge invariant.According to the previous discussion in [9], we can set a 1 = 1 and a 2 = a 3 = 0 in our case and note that the matrix C is the charge conjugation matrix defined as C = γ 2 γ 0 and u p (P ) is the proton spinor.The tensor σ µν is defined in terms of σ µν = i 2 [γ µ , γ ν ].As Eq. (42) indicates, there are in total 24 invariant functions S i , P i , A i , V i , T i , which we can not assign a definite twist.These calligraphic quantities can be related to the twist amplitudes in the following way: F(a 1 , a 2 , a 3 , (P • x)) = dα 1 dα 2 dα 3 δ(1−α 1 −α 2 −α 3 )e −i(P •x) i α i a i F (α i ) (43) The variables α 1,2,3 denote the momentum fractions of the different quarks inside the proton.We start with twist three contributions and relate the calligraphic quantities appearing in Eq. (42) with Eq. (43) to the definite twist amplitudes: F integrand on r.h.s. of (43) In the above table, the definite twist distribution amplitudes are given by At leading twist, there are in total three coefficients ϕ (0,±) 3 , which can be parameterized through the parameters f N , the normalization factor in leading conformal spin, and A u 1 as well as V d 1 , which belong both to the next to leading conformal spin: The remaining contributions can be classified according to their twist and specify their tensor structure [14]: twist 3 twist 4 twist 5 twist 6 Based on these considerations, we can order the calligraphic quantities in Eq. (42) into the different twist contributions.Thus, we the calligraphic quantities contributing to twist four are: F Integrand on r.h.s. of (43) Abbreviation S 1 S 1 For brevity, the renormalization scale dependence is dropped in the following discussion.Moreover, we follow the notation from [14][15][16]19].The twist four DAs in the conformal expansion are given by: Finally, for twist six we obtain the following contributions: F integrand on r.h.s. of (43) Abbreviations

B Form factors
In this section, we state the remaining expressions for the form factors before the zexpansion.

−0. 15 Table 5 :
Parameters of the z-expansion for the (b)-model form factors (in GeV 2 units) including all contributions to the nucleon DA up to twist six.Now one can use these input parameters to extrapolate the four different form factors and obtain these for the two models.Figures 2 to 5 show the form factors with respect to the Borel parameter M 2 .

Figure 2 :
Figure 2: Analysis for the form factor F (d)B→p R (q 2 ) in the (d)-model for different Borel parameter M 2 .The left plot shows all twist contributions combined for various choices of the threshold parameter s B 0 , while the right plot illustrates each twist individually.The mass of the dark matter particle m Ψ is set to the benchmark value m Ψ = 2 GeV[6].

Figure 3 :
Figure 3: Similar analysis as in figure 2 for the form factor F (d)

Figure 4 :
Figure 4: Analysis for the form factor F (b)

Figure 5 :
Figure 5: Similar analysis as in figure 2 for the form factor F (b)

2 R 6 Figure 6 :
Figure 6: Plot for the ratios R (d) B→p R ;i and R (d)

Figure 7 :
Figure 7: Plot for the ratios R (b) B→p R ;i and R (b) B→p L ;i defined in Eq. (33) and (34) after replacing (d) → (b) with respect to the dark fermion mass m Ψ .The right panel contains the ratios for the form factor F (b) B→p L , while The left panel shows the contributions from the form factor F (b) B→p R .

Figure 8 :
Figure 8: Branching fraction for the decay B → pΨ in the (d)-model with respect to the dark matter particle mass m Ψ , normalized to the effective four-fermion coupling |G (d) | 2 .The blue line with the blue dashed error band illustrates the original twist-three computation from[9], while the black curve shows the computation including contributions up to twist six.The dashed red curves represents the uncertainty on this calculation.

Figure 9 :
Figure 9: We show the same as in figure 8, but for the (b)-model.

Table 4 :
Parameters of the z-expansion in the case the (d)-model form factors (in GeV 2 units) including all contributions to the nucleon DA up to twist six.