A description of charmonia decays $J/\psi \rightarrow B\bar B$ within the QCD factorisation framework

Decays $J/\psi\to B\bar B$ into octet baryon-antibaryons pairs are studied within the QCD factorisation framework. The next-to-leading power correction to the leading-order decay amplitude is calculated for the first time. The phenomenological analysis also includes the interference with the electromagnetic decay $J/\psi\to \gamma^*\to B\bar B$. The branching fractions and the angular coefficients $\alpha_B$ are described with an accuracy of about $20\%$. The obtained results allow us to conclude that the calculated amplitudes provide a dominant contribution and therefore QCD factorisation provides a suitable basis for understanding of these decays.


Introduction
A systematic description of exclusive charmonia decays can be performed using QCD factorisation framework, see e.g. Refs. [1,2] and references therein. Within this approach one constructs the expansion of an amplitude in inverse powers of 1/m c combining non-relativistic QCD and hard-collinear effective field theory. As a rule, most of these calculations are performed with leading power accuracy. In many cases this gives a reliable description, but in many cases this approximation is insufficient even for a good qualitative description of the data. Potentially such discrepancies can be associated with various large corrections and with a complicated underlying QCD dynamics. Therefore, a description of exclusive charmonia decays still includes many open questions, see, for example, reviews in Refs. [3,4].
Decays of S-wave charmonia into proton-antiproton pairs have been discussed already long time ago within the QCD factorisation framework in Refs [1,5]. It was found that with the appropriate models of baryon matrix elements one can obtain a reasonable description of the branching ratio for J/ψ → pp. This observation was further used in order to constrain the nucleon DAs in Refs. [6,7] In recent years, the BESII and BESIII collaborations have provided a lot of new and accurate data in Refs. [8][9][10][11][12][13]. The description of these data is a non-trivial task, since requires calculations beyond the leading-power approximation, which are often are very problematic. In addition, the phenomenological consideration suggested in Refs. [14,15] indicates potential problems for the description based on QCD factorisation.
The new data provide not only the values of the cross sections, but also interesting information about the angular behaviour of the e + e − → J/ψ → BB cross section. This behaviour is sensitive to the amplitude, which is power suppressed within the systematic QCD description. The contribution of this amplitude to the decay width is suppressed by additional factor 1/m 2 c and therefore is expected to be a small correction. On the other hand, new data show that the effect of this amplitude on the description of the angular behaviour is very significant.
The various attempts to describe the angular coefficient α B have been considered in Refs. [16][17][18] using different phenomenological models and crude approximations. For the first time the power suppressed amplitude was calculated in the QCD factorisation framework in Refs. [19,20]. The qualitative phenomenological analysis carried out in these works indicates that the corresponding description works sufficiently well. However, in order to better understand the role of the power corrections some theoretical uncertainties need to be clarified.
In this paper we continue the investigation of the decays into baryon-antibaryon pairs within the QCD factorisation framework. Two additional effects are included into phenomenological consideration, which was ignored in Ref. [20]. First, we calculate the next-to-leading power correction ∼ Λ 2 /m 2 c to the dominant leading-power amplitude, where Λ is a typical hadronic scale. This correction is interesting because the mass of charm quark is not very large, and such contribution could provide rather large numerical effect. In addition, this correction is described by the same non-perturbative parameters as the subleading amplitude, which gives them better phenomenological constraints.
This power correction is closely associated with the higher Fock states of the baryon wave function and do not involve relativistic corrections with respect to the small heavy quark velocity in NRQCD. The decay mechanism associated with annihilation into the three hard gluons suggests that such contributions can also be described within the collinear factorisation framework that gives us an opportunity of a more detailed understanding of the systematic expansion in the inverse mass 1/m c . The higher twist baryon light-cone distribution amplitudes (DAs), associated with the non-perturbative collinear matrix elements, have already been investigated in Refs. [21][22][23][24][25][26][27][28]. These results allows one to reduce theoretical uncertainties associated with the non-perturbative QCD input.
In this work we also discuss the effect that occurs due to the contribution of the electromagnetic amplitude J/ψ → γ * → BB. The given amplitude is sensitive to the baryon time-like electromagnetic form factors. For a given energy scale these form factors can not be computed within the QCD framework because of relatively large effect from the soft-overlap rescattering (or Feynman) mechanism. Time-like baryonic FFs are also poorly studied experimentally. Different cross sections for the e + e − → BB channels have been measured by BESIII recent years, and the so-called effective FFs are obtained from data fitting, see Refs. [29][30][31][32][33][34][35][36]. We use these results in order to estimate the numerical influence of the electromagnetic admixture in the angular distribution α B .
Our paper is organised as follows. In Sec. 2 we discuss the higher twist baryon DAs and describe the models, which are used in our calculations. The calculation and the analytical expressions for the hard kernels and for the total amplitudes are discussed in the Sec. 3. The phenomenological analysis is presented in Sec. 4. In Sec. 5 we summarise our findings. In Appendix we provide additional useful details and a more detailed discussion of some technical issues.

.1 The basic set of the twist-5 DAs and their properties
In this section we describe the properties of twist-5 DAs and define the twist-5 auxiliary matrix elements, which are used in our calculation of the power corrections. We consider contributions with only three quark operators, which greatly simplifies the calculations. The total power correction is given by the sum of various contributions with the different combinations of DAs. The leading-power (lp) contribution to amplitude A 1 only includes the integral with the twist-3 DAs ϕ 3 , schematically: A (lp) 1 ∼ ϕ 3 * H * ϕ 3 , where the asterisk is used as a short notation for the convolution integral with respect to quark momentum fractions. The next-to-leading power (nlp) correction is more complicated, it depends on twist-4 ϕ 4 and twist-5 ϕ 5 DAs and can be represented as where the last term can be associated with the kinematical contribution. The various higher twist DAs and their properties have been discussed in Refs. [23,[25][26][27] and we use these results in the present work. Let us define the light-cone matrix elements, which define the light-cone DAs. The analysis of the twist-4 auxiliary DAs are given in Refs. [19,20,25]. Therefore, here we will focus only on the discussion of the twist-5 matrix elements.
We use the same notations as in Refs. [19,20]. The light-like vectors n,n ( n 2 =n 2 = 0, (nn) = 2 ) define the longitudinal and transverse components of a Lorentz vector V µ V + = (kn), V − = (Vn), V µ = V +n In the following we also assume that the baryon momentum k is directed along z-axis and can be expanded as For the baryon spinor N (k, s) we define the large and small components Nn and N n , respectively: The similar relations are also used for the antibaryon spinors. The relevant three-quark light-cone operators are constructed from the QCD quark fields q = {u, d, s} and from the light-cone the Wilson lines The basic operator can be written as where we use the short notation q(x − ) ≡ q(x − n/2). Below, for brevity, we also simplify the notation of the Dirac indices setting {α 1 , α 2 , α 3 } ≡ {1, 2, 3} and we also will not show the colour structure. Following to Ref. [26] we assume the following flavour content of the operators 0|q α 1 q α 2 q α 3 |B ≡ 0|q 1 q 2 q 3 |B : For simplicity, this will not be shown explicitly. Let us write the light-cone matrix element for the baryon state B as where the subscript "tw3" indicate that we only count in this term the contributions of twist-3 and so on. The leading twist-3 contribution read The twist-4 contribution is defined as The twist-5 contribution read where we always assume that σn n = σ αβn α n β . The symbol "FT" denotes the Fourier transformation with respect to the quark momentum fractions with the measure , the quark momenta in (14) are defined as The DAs defined in Eqs.(11)- (13) are symmetric or antisymmetric with respect to interchange with The DAs in Eqs(11)-(13) can be rewritten in terms of the basis light-cone DAs as [26,27] ( below we use short notation ( where The coefficients c ± B are chosen in such a way that these DAs satisfy the simple relations in the SU (3)-flavour symmetry limit, see e.g. Ref. [26].
For the nucleon state, the basic set of DAs can be further simplified by the isospin relationships, namely, it is sufficient to determine The twist-4 and twist-5 DAs can be decomposed into components corresponding to the contributions of the light-cone operators of geometrical twist-3, twist-4 and twist-5. Let us write such decompositions in the following form denote the genuine twist-4/5 contributions. The moments of genuine twist-5 DAs are not known.
In the present work we assume that the corresponding matrix elements are sufficiently small and therefore be neglected that yields The normalisation couplings f B , f B ⊥ , λ B 1,2 and λ Λ ⊥ have dimension of GeV 2 and their values will be discussed later.
The conformal expansion and evolution of the moments for different three-quark DAs was discussed in Refs. [22][23][24]. The explicit analytical expressions for the different terms in rhs in Eqs.(42)-(46) can be found in Refs. [23,26,27]. For convenience, in Appendix A we provide the explicit expressions for the models of basic DAs, which are used in present work.

Auxiliary DAs and their properties
Let us also consider the auxiliary set of the matrix elements, which is more convenient for the actual calculation of the hard kernels. More details about this calculation will be given in the next section. Here we only discuss the formal definition of the auxiliary matrix elements and their relations with the basic set of DAs.
The calculation of the decay amplitudes implies that the QCD fields can be splitted into sum of hard and collinear fields. The letter describe the non-perturbative overlap with the baryon-antibaryon state. Therefore performing a matching onto collinear operators one implies the effective field theory constructed from the hard (p h ∼ m Q , x ∼ 1/m Q ) and collinear fields. Such theory provides a definite power counting, which is used for the ranging of the operators. The power counting of the collinear fields is associated with their momentum, if we assume the collinear sector associated with the momentum k in Eq.(3) then (λ 2 ∼ Λ/m Q ) Performing the multipole expansions of collinear fields one also obtains derivatives of the collinear fields, which counts as Remind, that we do not consider the contributions with the quark-gluon operators, therefore terms like A ⊥ , A − , ∂ ⊥ A + and (n∂)A + will be excluded from our consideration. The auxiliary DAs can be defined in terms of the matrix elements of 3-quark twist-5 operators with the transverse derivatives. These operators are defined using the gauge invariant collinear fields, which have a definite scaling behaviour within the effective field theory framework Obviously The most general 3-quark twist-3 operator is where for simplicity we do not show the colour and flavour structure. The 3-quark twist-4 operators include one transverse derivative where the derivative is applied inside the brackets only. Notice that we do not need to write the covariant derivative because the field χ(x − ) is gauge invariant. The matrix elements of the twist-4 operators in (55) are considered in Ref. [19]. The results read where the color and the flavour structure of the operators are the same as described before. The DAs in the rhs of these equations can be expressed in terms of DAs {V i , A i , T i , S i , P i } B using the Lorentz symmetry and QCD equations of motion. This gives 1 [19,25] 1 If the all function in an equation have the same standard arguments (xi) ≡ (x1, x2, x3) we use the simplification F (xi) ≡ F to simplify the equations.
where the quark masses m i correspond to quarks q 1 q 2 q 3 in the matrix elements in Eqs.(7)-(9), respectively.
The additional relation is derived in Appendix B Combaining Eq.(69) with Eqs.(65), (66) and definitions in Eqs. (22) and (24) one finds for Υ B 4 the following equation Therefore the DA Υ B 4 is not an independent function in the 3-quark approximation. The 3-quark twist-5 operators include two transverse derivatives The longitudinal derivatives (n∂)χ can be excluded with the help of QCD equation of motion Therefore the basis set of the 3-quarks operators can be defined as The operators with [∂ ⊥ χ 3 ] can be excluded because where it is used that the matrix element from the operator with the total transverse derivative vanishes because the baryon momenta have longitudinal components only. The matrix elements of the three operators in Eq.(73) define the auxiliary set of the twist-5 DAs. It is convenient to write their definitions for the operators with projected Dirac indices {12}. We introduce the following set of the matrix elements with where we use the short notaion for the following combination This set of DAs can be splitted into two groups: chiral even and chiral-odd DAs defined in Eqs.(75)-(80) and Eqs.(81)-(83), respectively. These groups can be considered independently.
In order to establish relations of the auxiliary DAs with the basic ones, one needs to relate the operator O 123 defined in Eq.(6) with the operators in Eq.(73). This can be done with the help of the QCD equations of motion. Let us write the quark field as a sum where The small field η(x − ) can be written as, see e.g. Ref. [43] η( where dots denote the terms with insertions of the transverse gluon components W † (x − )i / D ⊥ W (x − ) , which give the quark-gluon operators and therefore can be omitted. Hence It is clear that each field η gives the one transverse drivative, therefore the twist-5 operators in the expansion of O 123 in Eq.(6) have two fields η The main idea is the following: we rewrite the terms in Eq.(90) using Eq.(89) in terms of the basis operators (73) and compute the matrix element using Eqs.(75)- (28). On the other hand, using the definitions (87) and one can express the [O 123 ] tw5 in terms of different projections of the three quark operator defined in Eq.(6) and use Eq.(10) in order to get the matrix element. Comparing the two results one finds the relations, which allows one to establish a connection of the auxiliary DAs with the basis ones. It is convenient to perform such calculations for the chiral-even and chiral-odd operators separately. Let us first consider the chiral-even sector, which includes 8 DAs: Let us consider as example the following operator where the Dirac indices {12} are contracted with the matrix C / n. Using Eq. (89) one finds where we use the short notation The matrix element of this operator can be computed with the help of Eqs.(79) and (80). One finds where the factor 1/(x 1 x 2 ) appears due to the inverse derivatives (in∂) −1 acting on the first and second fields. Remind, that the inverse derivatives (in∂) −1 must be undestood as nonlocal operator On the other hand, using Eqs.(87) one finds Computing the matrix element of this operator with the help of Eq.(13) one gets Comparing Eqs. (97) and (101) one obtains where, for simplicity, the arguments are not shown: F (x 123 ) ≡ F . The similar consideration for the operators ηCγ ⊥ ξη 3 and ξCγ ⊥ ηη 3 2 yields respectively. The analysis for the axial operators ηCnγ 5 ηξ 3 gives Two more equations can be obtained using the QCD EOM, which can be rewritten as and describes the free collinear quark. Consider the following operator with the total derivative Taking the matrix elements from the both sides one finds The same consideration for the axial operator χC / nγ 5 χχ 3 gives Therefore we obtain six linear algebraic relations. The investigation of other possible 3-quark light-cone operators using the same technique does not provide any new information. We did not find any method in order to get two more pure algebraic equations. However one needs at least two more equations in order to express the eight auxiliary DAs (91) in terms of the basis ones. In order to derive these equations we introduce off light-cone correlation function (CF). Performing the expansion of this CF around the light-cone we derive two more equations, which allow us to solve the problem. However these relations are already not algebraic. The derivation of the corresponding relations is more involved and the details are given in Appendix B. The obtained equations read where the factors 2(kx) = k + x − , 2(ky) = k + y − and the coordinates x and y enter in the definition of the FT in Eq.(14) 3 . Remind, that the DAs V i and A i on the rhs denote the auxiliary twist-4 DAs from Eqs.(63)-(64). For instance, for the nucleon, taking the asymptotic contributions of twist-3 and twist-4 DAs one finds Solving two equations (110) and (111) one finds expressions for T xy and G xy . Then one can easily solve the system of the linear equations. The properties of Eqs. (110) and (111) imply that in case of x = 0 or y = 0 the one lhs vanishes and therefore must vanish the corresponding rhs. This gives two non-trivial conditions ( These conditions allows one to perform the integrations by parts on rhs of Eqs. (110) and (111), that reduce these equations to a simple algebraic form. Notice, that the DAs on the rhs of Eqs. (113) and (114) depends on the twist-3 and twist-4 moments only. We have found that Eqs. (113) and (114) are satisfied if one uses for description of the DAs in Eqs. (113) and (114) the first moments (or equivalently the asymptotic contributions) for twist-3 and twist-4 DAs only. We believe that this observation is closely related to the disregard of the quark-gluon contributions and the genuine twist-5 moments Eq.(47). Such simple recipe does not work for the description of the twist-5 contributions associated with higher moments of the DAs and this makes such analysis to be more complicated in comparison with the twist-4 DAs. In order to illustrate this point we discuss a simple example, which is considered in the Appendix C. Therefore in this work we only use the asymptotic expressions for the twist-3 and twist-4 DAs for the description of auxiliary twist-5 DAs. Let us also mention that obtained description of the twist-5 DAs ensures their proper normalisations, which are consistent with the Lorentz structure for the local matrix elements. This is automatically provided by the algebraic equations derived above.
It is convenient to solve the equations for twist-5 auxiliary DAs using the explicit expressions for the asymptotic DAs. Then the solutions are given by the simple analytical expressions. Assuming below B = Λ we write the final result as The similar discussion also holds for the chiral-odd set of the auxiliary DAs. The details are given in Appendix B. The results read For the nucleon set one has to substitute in the formulas f ⊥ B → f N . The coefficients c ± Λ are defined in Eq. (29). It is useful to note that all auxiliary DAs are proportional to the factor x 1 x 2 x 3 that simplifies an investigation of analytic properties of the collinear integrals.
In addition to the identities, which follow from the Lorentz symmetry, one can also derive symmetry relations that follow from the condition that a baryon state has a certain isospin. For the auxiliary matrix elements this can be done in the same way as for the basis one, see e.g. Ref. [21]. For instance, for the nucleon (I N = 1/2) this gives 15 different relations for the auxiliary DAs. We provide these relations in Appendix D. These relations can also be used for the finding of the relation between the auxiliary and basic DAs but different baryons has different isospin, which complicates the analysis. At the same time, the relations that follow from the Lorentz symmetry are universal and valid for all baryons. Nevertheless, the isospin relations provide a powerful check of the obtained result. We have verified that the obtained expressions for the auxiliary DAs in Eqs.

Calculation of hard kernels
The decay amplitude can be written as [19] where ψ is the charmonium polarisation vector,N (k) and V (k ) denote baryon and antibaryon spinors, respectively. The scalar amplitudes A B i can be expanded within the effective field theory framework in small heavy quark velocity v 2 (NRQCD) and in powers λ 2 ∼ Λ/m c (collinear factorisation), where Λ is the typical hadronic scale. In this work we consider the leading-oder contribution with respect to small velocity v and compute the next-to-leading power (NLP) correction to the amplitude A B 1 , which is provided by the diagrams in Fig. 1a). The expansion of the amplitudes A B i in powers of λ is associated with the power counting of the collinear matrix elements, which describe the overlap with outgoing baryons. The structure of the collinear expansion allows one to write the following expansions The leading-power (LP) contribution scales as A lp 1 ∼ O(λ 8 ) because it is described by the two light-cone matrix elements, each of which is of order O(λ 4 ). The NLP correction A nlp 1 includes collinear operators, which are suppressed by the overall factor λ 4 and therefore A nlo 1 ∼ O(λ 12 ). For the second amplitude one finds A lp 2 ∼ O(λ 12 ) 4 . The expansion in powers of 1/m c mixes up power corrections from the expansions of the scalar amplitudes A B i and from the baryon spinors. The calculation of the amplitude in the factorisation framework yields the result in the following form where the coefficients δA i scale as The expressions for δA i are given by the convolution integrals of the hard kernels with the different baryon DAs.
In order to perform the interpretation of δA i in terms of the scalar amplitudes A B 1,2 one needs to perform the expansion of the rhs of Eq.(132). This gives where we use that k + ∼ k − ∼ m c . Comparing Eqs. (134) and (136) one finds The expressions for observables have a compact analytical form if we rewrite them in terms of the linear combinations G M and G E , which read Notice that the total power correction to A B 1 also includes the simple kinematical term δA lo Refs. [1,5]. The contribution δA lo 2 is obtained in Ref. [19]. The correction δA nlo 1 will be considered for the first time. In this calculation we only consider the contributions of light-cone 3-quark operators and do not include the various quark-gluon matrix elements.
The contribution of the sum of diagrams in Fig. 1 a) to δA nlo 1 can be written as where D i (123)(1 2 3 ) denotes the analytical expression of the diagram i in momentum space. It can be written as Here the matrices γ α ⊗γ β ⊗γ γ corresponds to the light quark-antiquark vertices, D (i)αβγ Q denotes the remnant part of the diagram, the trace in Eq.(140) is taken with respect of Dirac indices of the heavy quark fermion line. For simplicity, the colour structure of diagrams is not shown. The coupling f ψ describes the matrix element in NRQCD and can be related with the radial wave function at the origin The definition of the NRQCD matrix element is the same as in Ref. [19]. The expression (1− / ω)ε / ψ in the trace in Eq.(140) represent the Dirac part of the projector on the NRQCD matrix element, ω = (1, 0, 0, 0) and ε µ ψ denote the velocity and the polarisation vector of charmonium, respectively.
The operators P (x i ) andP (y i ) describe the projections on the matrix elements of the lightcone 3-quark operators for baryon and antibaryon, respectively. For the twist-3 case these operators depend on the nucleon spinors, Dirac matrices and twist-3 DAs and can be obtained just from the definition in Eq. (11). For instance, for the matrix element 0| O 123 |B(k) tw3 one easily obtains For the higher twist matrix elements these projections also include the derivatives with respect to quark and antiquark momenta k i and k i , which are assumed to be off-shell k 2 i = 0, k i = 0. Only after the differentiations one can neglect the small momentum components and set k i x i k +n /2 and k i y i k − n/2. The derivation of the twist-4 projections have been considered in Ref. [19] and here we follow the same technique. The analytical results for the different projectors are given in Appendix E. Substituting the various projectors in Eq.(140) one obtains the expression, which allows one to calculate the hard kernels. The analytical calculations were performed using the package FeynCalc [37,38].
The obtained result can be written as where the normalisation coefficient A 0 read The dimensionless integrals J B ij have subscripts "44" and "35" , which are related with the structure of the collinear matrix elements: twist-4×twist-4 and twist-3×twist-5. Each integral in Eq.(145) is given by the sum of simpler integrals J[X, Y ] where X and Y denote the DAs, which are included in the integrand where K XY (x i , y i ) is the universal hard kernel. The calculation of the traces in Eq.(140) gives the sum of the collinear integrals, which can be simplified using the symmetry properties of the DAs and the hard kernels with respect to interchange x 1 ↔ x 2 and x i ↔ y i . Below we provide the analytical expressions for the integrals in Eq.(145), where such simplifications are taken into account. Let us write each term in Eq.(145) as The subscript "even/odd" corresponds to the chiral-even and chiral-odd subsets of the DAs. Then  The analytical expressions for the hard kernels K XY (x i , y i ), as defined in Eq.(147), are given in the Appendix F. It has been analytically verified that these kernels correctly describe the relations between the amplitudes in the SU (3) limit. The corresponding convolution integrals can be relatively easily calculated. Substituting the explicit expressions for DAs with arbitrary moments and summing the singular integrands gives analytical expressions with the well defined four-dimensional integrals. These integrals can be easily computed numerically using the standard mathematica packages for numerical integrations.

Phenomenology
A systematic numerical analysis of existing data is still very challenging due to various theoretical and experimental uncertainties. Therefore, the main goal of the following discussion is to understand how well the existing data can be described using plausible assumptions about various nonperturbative parameters.
In this analysis we take into account the electromagnetic amplitudes, which are shown by the diagram in Fig.1 (b). With these contributions the decay amplitudes read where F B 1,2 are the time-like e.m. form factors. Expression for the width reads The amplitudes G M and G E are defined in Eq.(139), the time-like Sachs form factors G B E,M read The coupling ζ in Eq.(161) is defined as Therefore, the observables depend on the time-like form factors of baryons, which are still not well known. For the given value of the q 2 = M 2 ψ their values are still depend on the large nonperturbative corrections and can not be sufficiently well estimated within the QCD factorisation framework. This introduces additional uncertainty into the phenomenological analysis.
The hadronic amplitudes calculated to this accuracy within the factorisation framework are real. The e.m. FFs have real and imaginary parts, therefore to this accuracy one finds of data in Ref. [31]. In other cases we assume that |G B M | ≈ G B eff . The value of cos ϕ M can be estimated from the data for the width (160), where γ B is fixed by the experimental value. An estimate of G B E is the most challenging point. We assume that G B E ≤1.5 G B eff , where the numerical coefficient 1.5 is our conservative estimate. Numerical value of ζ 2 |G E | 2 / |G E | 2 is quite small and can be neglected, which yields Therefore, assuming the estimate we get an interval, which is considered as an estimate of the possible numerical effect provided by R B em . Let us now discuss the non-perturbative input, which is used for calculation of the hadronic amplitudes. Recall that we consider only three-quark DAs. It is convenient to describe the non-perturbative input in terms of basic DAs as discussed in the Sec.2.1. Analytical expressions for the corresponding models are described in Appendix A. The values of DAs parameters are summarised in Table 1. The parameters of nucleon chiral-even DAs correspond to the model ABO-I from Ref. [25]. This set of parameters have beed obtained from the fit of the space-like electromagnetic FF data using the light-cone sum rules. The moments λ 2 and ξ 10 are evaluated in Ref. [21] using the QCD sum rule approach. Some twist-3 moments for different baryons are also obtained using QCD sum rules in Ref. [5]. Recently, the different twist-3 and twist-4 moments have also been computed on the lattice, see Ref. [28]. We use these results in order to fix parameter values in DA models. In particular, the values for the moments φ B 10 , φ B 11 , λ B 1 ,λ B ⊥ , λ B 2 , π B 10 and π B 11 (B = N ) are fixed according to data in Ref. [28].
For the nucleon twist-3 normalisation coupling we take the value f N (1 GeV) = 5.5 × 10 −3 GeV 2 , which is in agreement with the sum rule estimate f N (1 GeV) = (5.3 ± 0.5) × 10 −3 GeV 2 obtained in Ref. [5]. The other twist-3 normalisation couplings f B are fixed from the fit of data. The values obtained in Ref. [28] were used as the first estimate and then changed if this improves the description of the data. Therefore the resulting value f Λ in Tab.1 is larger than the estimates in Refs. [5,28], which read 0.50 and 0.48, respectively. 6 We also find a relatively small difference for values of f Σ (0.54 [28] and 0.48 [5]) and f Ξ ⊥ (0.64 [28] and 0.50 [5]). The values for f Σ ⊥ and f Ξ are the same as in Ref. [28]. The values of the twist-4 moments η B 10 , η B 11 , ζ B 10 , ζ B 11 and ξ B 10 are fixed from the data, assuming that their values are not significantly different from the nucleon ones because of approximate SU (3) symmetry.
One of the elements of the current analysis is the estimation of the power-law correction to G B M , see Eq.(139). In Table 2 we show the numerical results obtained for these corrections. The total effect of the power corrections is quite moderate due to the partial cancellation in the total sum, see Eq.(139). Therefore the total numerical effect from the next-to-leading power correction varies from 1.7% for Ξ up to 16% for nucleon. The numerical results for branching fractions and ratios γ B are presented in Table 3. The obtained values of γ B include the error bars, which describe the estimate of the uncertainty occurring from to the interference with the electromagnetic FF's as described by the coefficient R em . The given choice of the different parameters for baryon DAs allows one to describe the data for the γ B within the accuracy 20%, which is sufficiently reasonable taking into account various higher order effects. The largest differences are obtained for neutron and for the Σ decay channels. In considered description, the neutron channel differs from the proton one only in the contribution of electromagnetic form factors. However this interference can not provide a sufficient numerical impact. On the other hand the neutron data in Ref. [9] have very large systematic error, which gives γ n = 0.95(6) ± 0.27 syst . Therefore the obtained value of γ n can be accepted as quite reasonable result. For Σ decay channel, it is possible to obtain the larger value of γ Σ but corresponding parameter set gives a worse description of the branching ratio because of relatively large and negative power correction.
Comparing the current analysis with that in Ref. [20] we conclude that the NLP corrections improve the phenomenological description and make it possible to better constrain the unknown parameters of baryon DAs. The negative effects of the power contribution reduces the absolute value of the amplitude G M and increases the value of the ratio γ B that improves the phenomenological description. However, this effect cannot be very large, since this can significantly worsen the description of the branching fractions.

Conclusions
In conclusion, in this work we continue to develop the analysis proposed in Ref. [20]. The decay amplitudes J/ψ into baryon-antibaryon pairs are considered within the NRQCD and collinear factorisation approach, power suppressed corrections and electromagnetic contributions are taken into account. This allows us to better appreciate some of the theoretical uncertainties. The next-to-leading power correction, which is described by twist-4 and twist-5 three-quark baryon DAs is calculated and discussed for the first time. It is shown that this contribution is well defined and this allows one to study their effects in a more systematic way. From the phenomenological point of view the NLP correction allows to improve the description of the data and provide more rigorous constrains on the baryons DAs parameters.
In the given analysis we also take into account the interference with the electromagnetic amplitude, which depends on the time-like baryon form factors. Unfortunately, these values are little known due to the lack of experimental information. Therefore, an evaluation of their numerical impact can only be made under certain assumptions.
The existing experimental data for decays J/ψ → BB provide an interesting information about two observables: branching fractions and angular distribution coefficients α B . The value of α B is determined by the ratio of the two independent decay amplitudes, which is especially interesting because it does not contain many theoretical uncertainties that are significant for the branching fractions. This also provides an attractive opportunity to study the baryon distribution amplitudes.
One of the main goal of the given phenomenological analysis is to verify that the current data can be described within the factorisation framework. In order to constrain the non-perturbative baryon DAs we use results obtained in the framework of QCD sum rules [5,25] and from the lattice calculations [28]. Some of unknown moments have been estimated using flavour SU (3)symmetry and from data fit. It is found that under reasonable assumptions about the nonperturbative parameters, the data can be described with an accuracy of 5% − 20%.
It is shown that the power corrections arising from higher twist baryonic DAs are moderate. This conclusion is, of course, sensitive to the models of baryonic DAs. The maximal numerical effect is about 16% is obtained for nucleon only. It is interesting to note that the model ABOI proposed in Ref. [25] for describing nucleon electromagnetic form factors also works well for describing charmonium decay, if we use the value for the normalisation couplings f N obtained from the QCD sum rules [5]. The values of the branching fractions are very sensitive to the normalisation scale, a reliable description is obtained only for the relatively low scale µ 2 = 1.5 GeV 2 .
Note also that a large value of γΣ ∼ 2 leads to a large numerical effect, which is crucially important for the correct estimate of the decay width. Unlike other baryons, in the case of Σ the contribution of the power suppressed term in the expression for the width, see the equation (160 ), is of the order of one, because it is enhanced by a factor of four due to large γ Σ .
The angular coefficient γ B is quite sensitive to the interference with the electromagnetic decay amplitudes. Under the described assumptions, the corresponding numerical effect is estimated in the range of 10% − 25%. This definitely reduces the accuracy in the constraining of DA parameters. Nevertheless, the obtained results allow us to conclude that the calculated hadronic amplitudes make a dominant contribution and therefore the QCD factorisation provides a suitable basis for understanding dynamics of these decays.

Acknowledgements
This work is supported by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) Project-ID 445769443.

Appendix A Analytical expressions for the models of baryon DAs
In this Appendix we provide the explicit expressions for the models of the basic light-cone baryon DAs, which define the QCD nonperturbative input in the compact form. These DAs are used in order to compute the auxiliary DAs, which are convenient for the calculation of the decay amplitudes. Recall, that in this work we follow the definitions from Ref. [25].
We start from the description of the twist-3 DAs. For the nucleon DA introduced in Eq.(30) we use the following expression For other baryons B = N the DAs are defined in Eqs. (19) and (20). Corresponding models read The orthogonal polynomials P 21 (x i ) read The moments φ B ik are multiplicatively renormalisable and the corresponding anomalous dimensions can be found in Ref. [25].
The set of the models for the twist-4 DAs are defined as following. The nucleon DAs from Eqs. (31) and (32) have the following twist decomposition The kinematical twist-3 contributions read where the orthogonal polynomials R 1i read The genuine twist-4 contributions read The basic twist-5 DAs are used in order to obtain the auxiliary DAs, which are discussed in Sec.
Numerical values of all moments are given in the Tab.1. All these moments are multiplicatively renomalisable and the corresponding anomalous dimensions can be found in Refs. [23,25].

B Relations between the twist-5 DAs
Here we discuss the technical details that clarify the derivation of relationships between the various DAs introduced in Sec.2. Let us start from discussion of the chiral even DAs. There are two ways to get additional equations, which relate the auxiliary and the basic DAs. One way is to consider the light-cone operator with one derivative and to repeat the same analysis as in Sec.2. The second method is to use the off light-cone correlation function (CF). Below we use this one. For simplicity we consider the nucleon case. We also do not write explicitly the colour indices and the Wilson lines.
Consider the following CF (x 2 = 0, y 2 = 0) The twist-5 operators resulting from the expansion can be divided into three groups. The first group includes the operators with two transverse derivatives, which arise from the expansion of the arguments of the fields There second group includes the terms with one transverse derivative and one small field η The required relations can be obtained from the matrix elements of the operators from the second group (B.4). Consider the terms with x ⊥ . There are three twist-5 light-cone operators, which appear in this case where we assume qqq 3 ≡ q(x − )q(y − )q 3 (0). The matrix elements of these operators can be obtaned using QCD EOM (88) and definitions Eqs.
where we assume 2(ky) = k + y − . The similar analysis for the operators with y ⊥ allows one to obtain one more relation Consider now the axial CF The analogous consideration of terms with x ⊥ and y ⊥ yields the two relations 14) It is convenient to combine these two relations with ones in Eqs.(B.9) and (B.11) and simplify the obtained expressions using the relations (102)-(105). This gives the relations given in Eqs. (110) and (111). Consider now the auxiliary chiral-odd DAs. The set of these functions is defined in Eqs.(81)-(83) and reads The algebraic relations can be obtained in the same way as described for the chiral-even DAs in Sec.2.2. Therefore we skip the details and only provide a schematical description. The consideration of the scalar operators ηCξη 3 and ξCηη 3 yields respectively. From the analysis of the pseuvdoscalar operators ηCγ 5 ξη 3 and ξCγ 5 ηη 3 one finds One more relation follows from the consideration of the tensor operator ηCnγ ⊥ ηξ 3 : Using the same consideration as in Eq.(107) we consider the matrix element of the operator (n∂) [ξCnγ ν ⊥ ξξ 3 ]. This yields All other operators give the linear combinations of these six relations.
In order to get one more relation we use the method of CF, which is described above. Consider the following CF where the rhs includes all terms up to twist-5 accuracy. Expanding the operator in the lhs we consider the following linear in x ⊥ and y ⊥ twist-5 terms The consideration of the corresponding matrix elements yields The sum of these two relations gives The consideration of the matrix elements 0| O ασ x |B(k) and 0| O ασ y |B(k) yields respectively. The sum of these two equations gives  gives where we assume x 3 = 1 − x 1 − x 2 . This provides a good check of the obtained expressions.
Finally, let us also mention that the additional relation for the twist-4 chiral-odd DAs given in Eq.(69), can be obtained from the consideration of the matrix element of the operator ξ(x − )Cσ µν ξ(y − )η 3 (0).

C On the determination of quark and quark-gluon contributions for the next-to-leading power correction
In this section, we discuss a simple example that demonstrates a difficulty in the determination of the pure quark contribution for the subleading power corrections. We assume the same kinematical notations as for the baryon with momentum k. Let us consider the light-cone matrix element for the vector meson state (J P C = 1 −− ), which can be parametrised as where e µ ≡ e λ µ (k) denotes the polarisation vector, FT denotes the appropriate Fourier transformation All defined DAs are normalised to one The twist-3 DA g v and the twist-4 DA g 3 in the rhs have been studied in Refs. [41,42]. The obtained results read whereḡ v andḡ 3 denote the contributions of quark-gluon matrix elements. The local limit in Eq.(C.1) gives where it is used that The off light-cone CF can be written as [41] 0|ψ The derivation of the expression for the DA A is given in Ref. [41], the result reads whereĀ again denotes the quark-gluon contributions. Such notation for the contributions of quark-gluon matrix elements is also accepted below. All results in Eqs.(C.4),(C.5) and (C.10) are obtained with the help of differentiation of the CF, see details in Refs. [41,42]. Let us now consider the twist-4 DAs with help of the technique discussed in this work. The twist-4 operators, which appear in the expansion of the off light-cone operatorψ(x)γ σ ψ(0) read where, remind, the collinear fields ξ and η are defined in Eq.(87). The operator O σ 4 [A] denotes the contribution of the quark-gluon operators.
The matrix element of the twist-4 operator with two derivatives can be obtained from the expansion of the CF Consider now the operators with the longitudinal derivative i(n∂) in (C.11). Using QCD EOM where dots denote the quark-gluon contributions. Similarly The sum of these terms gives the operator m.e. of the operator i(n∂) ξ (x − )/ nξ(0) with the total derivative, which can be computed from Eq.(C.9). This gives the following relation A(u) = 4uūφ(u) +Ḡ 1 (u), (C. 16) whereḠ 1 (u) again schematically denotes the contributions of quark-gluon DAs. Consider now the m.e.
where we show the two-particle quark-antiquark operator only. On the other hand (C.18) Therefore this gives whereḠ 3 denotes a contribution from quark-gluon DAs. Using (C.16) one obtains Neglecting the quark-gluon termsḠ 1 ∼Ḡ 3 → 0 one finds the result that agrees with Eq.(C.5) assummingḡ 3 → 0. At the same time the results for A in Eqs.(C.10) and (C. 16) do not agree if one simply neglects the quark-gluon contributions This result can be explained by the fact that the neglected quark-gluon contributions still implicitly depend on the twist-2 moments. For the local matrix elements this can be seen from the QCD equation of motions, see e.g. Ref. [41]. Therefore the simple recipe prescribing just to neglect the quark-gluon contributions, as for the twist-3 operators, is not applicable for the twist-4 calculations. Notice however that Eq.(C.22) is valid for one particular case: for the asymptotic shape φ as (u) = 6uū. The corresponding local operator can not be related with any quark-gluon one. This observation shows that consistent twist-4 calculation, which only account the twoparticle quark operator requires some prescription, which makes it possible to unambiguously and systematically distinguish between the twist-4 quark and the quark-gluon contribution.
In fact, the Lorentz symmetry of the CF (C.9) implies that DA A must satisfy to a certain condition. Let us briefly discuss this point. Consider the following matrix element where we used Eq.(C.1). Now the operator in lhs can be rewritten as where we again skip the quark-gluon operators. The matrix element of this operator gives where Aq q denotes the pure quark part of the DA A A = Aq q +Ā. (C.27) Therefore one finds whereḠ 4 schematically denotes the quark-gluon contributions. Comparing with Eq.(C.23) one finds the following relation In conclusion, the given consideration demonstrates that in general case the definition of the pure quark contribution for the twist-4 DAs is not universal due to the implicit contributions from the matrix element of the quark-gluon operators. In the given example the universal part can only be defined in terms of the asymptotic part of the twist-2 DA φ because the corresponding local operator do not appear from the quark-gluon ones.
The similar arguments also applicable for the three-quark operators. Therefore in our calculations we only consider the asymptotic contributions of the twist-3 and twist-4 3-quark matrix elements in order to define the corresponding universal three-quark components of twist-5 DAs.

D Isospin relations between the auxiliary twist-5 DAs
The standard analysis of the light-cone matrix elements with two transverse derivatives allows one to obtain many relations between the auxiliary DAs defined in the Sec.2.2. Here we consider such relations for a nucleon state. For simplicity, in this section we omit the baryon subscript "B=N" assuming for all DAs F N ≡ F . Let us define the operator basis as following In these formulas we assume that each quark field q = u, d is projected onto the large component ξ as defined in Eq.(87). Recall, that we use the short notation for the Dirac indices q α i ≡ q i . From the technical point of view the derivation of such relations is the same as in Ref. [21] . This method gives the following three relations for the light-cone matrix elements of defined operators For the DAs, which describe the symmetric traceless matrix elements the set of relations read T 2xy −P 23 (T 2xy + T 2xx ) −P 13 (T 2xy + T 2yy ) = 0. (D.8) T 2yy +P 13 T 2yy +P 23 (T 2xx + T 2yy + 2T 2xy ) = 0. The relations for the symmetric matrix elements are described by the matrix elements of the operators The analysis of the operator ∂ 2 ⊥ u ud gives The analysis of the operator u ∂ 2 ⊥ u d gives The analysis of the operator [∂ ⊥ u] [∂ ⊥ u] d gives From these equations it follows that the relations between the chiral-odd and chiral-even auxiliary DAs have the following form

E Expressions for the projections onto higher twist matrix elements
Here we give expressions for the projections P andP that enter Eq.(140). In the following we assume the matrix element with the baryon in the initial state as in the definitions in Sec.2.1. This is convenient for a comparison of definitions of baryon DAs, which exist in the literature. This also implies that one calculates the amplitude of the time reversal process B +B → J/ψ.
Since the corresponding amplitude is real this also gives the correct result. The antibaryon projectors can be obtained from the baryon ones using the charge conjugation. For partonic diagrams we assume the configuration of the external momenta as shown in Fig.2. After differentiations with respect to momenta k i and k i we set The derived projectors depend on the auxiliary DAs, which are defined in Sec.2.2. The twist-4 projectors have already been derived in Ref. [19]. The chiral-even projection reads where (i = 1, 2) The chiral-odd projection reads We assume that derivatives in these formulas are applied to the expressions on the right side, i.e. the derivatives are not applied to the momenta from the left side and therefore corresponding momenta can be simplified neglecting the power suppressed components. For a convenience, we also show the Dirac matrices in the square brackets: [/ k], [/ kγ α ⊥ ] etc. The expressions for the twist-5 projectors are more complicated and include some special contribution, which is completely defined in terms of twist-3 DAs.
The chiral-even projector is defined as The terms in the second line of Eq.(E.10) read . (E.18) The terms from the third line of Eq.(E.10) The terms from the last line of Eq.(E.10) Notice that terms C V 1 and C A1 include the derivatives with respect to initial momentum P , see a more detailed explanation below.
The twist-5 chiral-odd projector is defined as where the coefficients read terms only depend on the twist-3 DAs V B 1 , A B 1 and T B 1 and they are related to the matrix element of the leading-twist operator with the total derivative (z 3 n) 0| (n∂) [ξ 1 ξ 2 ξ 3 ] |B(k) = −i(nk)(z 3 n) 0| ξ 1 ξ 2 ξ 3 |B(k) , (E. 38) where (nk) ≡ k − = 2m 2 B /k + . This contribution does not vanish, in contrast to an operator with the total transverse derivative. It is clear that such contributions generate power corrections associated with the power suppressed component of the total baryon momentum k. Notice, that the calculation of the contributions with derivatives ∂/∂k 1,2 implies that k 3 = k − k 1 − k 2 as it follows from the momentum conservation. The contributions with the matrix elements as in Eq.(E.38) can be rewritten in a special way. This is related with the details of the Fourier transformation and with the derivation of the momentum conservation δ-function.
Let us consider, as example, the contribution with V 1 (x) and V 1 (y i ) and the diagram as in Fig.1(a). The following discussion is also applicable for other cases. Performing the standard steps of the calculation one obtains the following result In order to extract the momentum conservation δ-function from the expression for T V 1V 1 one uses that x 3 = 1 − x 1 − x 2 , y 3 = 1 − y 1 − y 2 and the translation invariance of the diagram D[z 123 , z 123 ]. Performing redefinition z i → z i + z 3 and z i → z i + z 3 for i = 1, 2 one defines the Fourier exponent with the coordinate z 3 e iP ·(z 1 +z 2 )/2 e −i(k 1 z 1 )−i(k 2 z 2 )−i(k 3 z 3 ) e −i(k 1 z 1 )−i(k 2 z 2 )−i(k 3 z 3 ) → e i(P z 3 )−i(k − z 3+ )/2−i(k + z 3− )/2 e iP ·(z 1 +z 2 )/2 e −i(k 1 z 1 )−i(k 2 z 2 ) e −i(k 1 z 1 )−i(k 2 z 2 ) , The integrals over dz i dz j in this expression can be understood as FT of the diagramD dz 1 dz 2 dz 3 e iP ·(z 1 +z 2 )/2 e −i(k 1 z 1 )−i(k 2 z 2 ) dz 1 dz 2 e −i(k 1 z 1 )−i(k 2 z 2 )D [z 12 , z 12 ] = D(k 12 , k 12 , P ). (E.44) The first line in Eq.(E. 43) gives the δ-function if the factor k − (−iz 3 · n) is rewritten as where l is just an auxiliary momentum. This gives T V 1V 1 (x i , y i ) = k − n α ∂ ∂l α δ(P − k − n/2 − k +n /2 − l) D(k 12 , k 12 , P ) l=0 (E.46) = k − −n α ∂ ∂P α δ(P − k − n/2 − k +n /2) D(k 12 , k 12 , P ) P =2mcω (E.47) = δ(P − k − n/2 − k +n /2) 2m 2 B k + n α ∂ ∂P α D(k 12 , k 12 , P ) P =2mcω where one assumes that P 2 is not fixed before the differentiation. This result explains the computational prescription for the terms with derivatives ∂/∂P . Let us also to note that such terms also include the endpoint singularities, which cancel in the sum of the all twist-5 contributions. Therefore, it is important to take them into account properly. The twist-3×twist-5 integrals have the following kernels.