Higher-twist B-meson Distribution Amplitudes in HQET

We present a systematic study of higher-twist distribution amplitudes (DAs) of the B-meson which give rise to power-suppressed $1/m_B$ contributions to B-decays in final states with energetic light particles in the framework of QCD factorization. As the main result, we find that the renormalization group equations for the three-particle distributions are completely integrable in the large $N_c$ limit and can be solved exactly. General properties of the solutions are studied. We propose simple models for higher-twist DAs which satify all existing constraints and can be used in phenomenological studies.


Contents
1 Introduction B-meson light-cone distribution amplitudes (DAs) are the main nonperturbative input to the QCD description of weak decays involving light hadrons in the framework of QCD factorization [1][2][3]. The simplest light antiquark-heavy quark DA gives the dominant contribution in the expansion in powers of the heavy quark mass and received much attention already [4][5][6][7][8][9][10]. It has become increasingly clear, however, that the leading power accuracy is not sufficient and the power of the QCD factorization approach depends crucially on the possibility to control, or at least estimate, the corrections suppressed by powers of the b-quark mass. This task is very challenging due to infrared divergences which appear in power-suppressed contributions in the purely perturbative framework. In recent years there had been some progress in this direction based on the combination of light-cone sum rule approach with the expansion in terms of B-meson DAs [11][12][13][14][15][16]. In this technique the so-called soft or end-point nonfactorizable contributions to B-decays can be calculated in terms of the DAs of increasing twist. One of the problems on this way is that higher-twist B-meson DAs involve contributions of multiparton states and are practically unknown.
In the recent paper [17] it was pointed out that the structure of subleading twist-three B-meson DAs is simpler than expected. In particular the twist-three DA φ − (ω) evolves autonomously with the scale and does not mix with "genuine" three-particle contributions. In this work we extend this analysis to twist-four DAs. Our main result is that the corresponding renormalization group equations (RGE) are completely integrable and can be solved exactly. Combined with the relations that follow from QCD equations of motion, this structure provides one with a set of robust constraints on the DAs and allows one to build phenomenologically acceptable models with a minimum number of free parameters.
The presentation is organized as follows. In Sect. 2 we present the general classification of the existing two-and three-particle B-meson DAs and discuss the corresponding twist and conformal spin assignment. We also use this Section to introduce our notation. Sect. 3 contains our main results. We discuss here the scale dependence of the three-particle twist-four DAs, find explicit analytic solution and work out the representation of the DAs in terms of the eigenfunctions of the evolution kernels. The derivation uses the formalism of the Quantum Inverse Scattering Method (QISM) [18][19][20] and is sketched in Appendix A; details will be published elsewhere. The two-particle higher-twist DAs that appear in the light-cone expansion of the heavy-light correlation functions are discussed in Sect. 4, based on the earlier work by Kawamura et al. (KKQT) [21]. The relations between different DAs due to QCD equations of motion (EOM) are discussed in detail. Two new relations between the three-particle DAs are derived in App. C using a different technique based on the two-component spinor formalism. Several simple models for the higher-twist DAs that satisfy EOM constraints are introduced and discussed shortly in Sect. 5. The final Sect. 6 contains a summary of our results and a short discussion of the remaining problems.

General classification
Following [4] we define the B-meson DAs as matrix elements of the renormalized nonlocal operators built of an effective heavy quark field h v (0) and a light (anti)quark at a light-like separation: is the Wilson line factor that ensures gauge invariance. Such factors are always implied but will often be omitted for brevity.
In practical calculations the gluon field strength tensor is often contracted with the light-like vector. The definition in (2.7) leads to the following pair of equations: cf. [21]. It is easy to show that but the X, Y -and X, Y -functions are not related to each other. Note that the DAs W and Z do not appear in these expressions.
In what follows we will use shorthand notations for the field coordinates on the light cone:q (z 1 ) ≡q(z 1 n) , G µν (z 2 ) ≡ G µν (z 2 n) .

Collinear twist and conformal spin assignment
The basis of the DAs in (2.7) is convenient because of the simple Lorentz structures. However, it is not suitable for discussion of QCD factorization as the DAs in this basis do not have definite (collinear) twist. Hence terms with different power suppression in the heavy quark expansion get mixed.
Twist t and conformal spin j of the light quark and gluon fields are given by the usual expressions where d is the canonical dimension and s is the spin projection on the light cone. Twist of a nonlocal heavy-light operator can then be defined as the sum of twists of the light constituents plus one unit of twist for the effective heavy quark field h v . Adding one unit of twist for h v is entirely a convention which we adopt in order to match the usual twist hierarchy for light quark-gluon operators; in this way in both cases the leading-twist contributions are defined as twist-two. The DAs of definite twist and conformal spins of the constituent fields are most easily defined by the corresponding projections of the general expression in (2.7). One finds one DA of twist three three independent twist-four DAs where where and one twist-six DA The conformal spin assignment for all DAs is summarized in Table. 1. The twist-five and twist-six DAs are not expected to contribute to the leading power corrections O(1/m B ) in B-decays and will not be considered further in this work.
Note that we also do not consider twist-four four-particle B-meson DAs (with two gluon fields and/or with an extra quark-antiquark pair). By analogy to the DAs of light mesons twist-four four-particle DAs are expected to be small and, most importantly, have autonomous scale dependence i.e. they do not mix with the three-particle DAs. Thus they can be consistently put to zero at all scales and do not reappear via evolution.

Asymptotic behavior at small momenta
Asymptotic behavior of all DAs at small quark and gluon momenta is determined by conformal spins of the fields [24] f (ω 1 , ω 2 ) ∼ ω 2j 1 −1 In particular These expressions can be verified considering correlation functions of the corresponding light-ray operators and suitable local currents, e.g. [11], and are stable against evolution provided the renormalization group equations respect conformal symmetry which is true to the leading logarithmic accuracy.

Spinor representation
Discussion of scale dependence of the DAs is considerably simplified using spinor formalism.
In this work we follow the conventions adopted in [25,26].
Any light-like vector can be represented by a product of two spinors. We write n αα = n µ σ µ αα = λ αλα ,n αα =n µ σ µ αα = µ αμα (2.22) whereλ = λ † ,μ = µ † . We choose which is consistent with our normalization (nn) = 2. 1 The "+" and "-" fields are defined as the projections onto the auxiliary λ and µ spinors, is written in this notation as The equation of motion (EOM) for the effective heavy quark field The gluon strength tensor F µν can be decomposed as Here f αβ andfαβ are chiral and antichiral symmetric tensors, f * =f , which belong to (1, 0) and (0, 1) representations of the Lorenz group, respectively. Rewriting the relevant operators in spinor notation one obtains (f → gf ) Since the contributions of left-and right-handed (chiral and anti-chiral) quarks have to be equal, one can drop half of the terms, e.g., the ones involving ψ-spinor, for most purposes. Note that Φ 4 and Ψ 4 + Ψ 4 contain light quark and gluon fields of the same chirality, whereas in Φ 3 and Ψ 4 − Ψ 4 chirality of the light degrees of freedom is the opposite. Since chirality (and twist) is conserved in perturbation theory, we expect that Φ 4 can get mixed under evolution with Ψ 4 + Ψ 4 , but the scale dependence of the "genuine" twist-four contribution to Ψ 4 − Ψ 4 is autonomous.
For completeness we write also twist-five and twist-six DAs in the spinor representation: and (2.31)

Scale Dependence
DAs are scale-dependent and satisfy renormalization group equations with the evolution kernels that can be found in Refs. [26,27] (in position space). They are collected in Appendix D. The corresponding expressions in momentum space can be found in Ref. [28]. The evolution equation for the leading-twist B-meson DA Φ + was derived by Lange and Neubert [5] and solved in Refs. [8,9]. The evolution equation for the twist-three DAs Φ − and Φ 3 was constructed and solved in the large-N c limit in Ref. [17] using a "hidden" symmetry of this equation called complete integrability. It turns out that evolution equations for the twist-four DAs are completely integrable as well and can be solved in the same manner. The corresponding expressions present the main result of our study. In this section we present the final expressions. The derivation uses the formalism of the Quantum Inverse Scattering Method (QISM) and is sketched in Appendix A; details will be published elsewhere.
For the two-particle twist-two and twist-three DAs one obtains [3,8,9,17] in position space, and in momentum space, respectively. The coefficient functions η + (s, µ) and η (0) 3 (s, µ) contain all relevant nonperturbative information and have to be fixed at a certain (low) reference scale µ = µ 0 . 2 The most important parameter for the QCD description of B-decays is the value of the first negative moment The scale dependence of η + (s, µ) and η 3 (s, µ) is given by where L = α s (µ)/α s (µ 0 ) and Here s 0 = e 5/4−γ E , Γ cusp (α s ) = αs π C F + . . . is the cusp anomalous dimension and we have factored out the scale dependence of the B-meson decay constant (3.6) The three-particle twist-three DA Φ 3 (z 1 , z 2 , µ) satisfies a more complicated renormalization group (RG) equation where H 3 is a certain integral operator that can be written as a sum of two-particle kernels that can be found in Appendix D. This equation was solved in Ref. [17] in the large-N c limit, i.e. neglecting corrections to H 3 that are suppressed by a factor 1/N 2 c . In [17] the eigenfunctions of the evolution equation and the corresponding anomalous dimensions have been found.
The DA Φ 3 (z 1 , z 2 , µ) can be expanded in terms the eigenfunctions of the large-N c evolution kernel as follows [17]: where du uū e is(u/z 1 +ū/z 2 ) . Note that the eigenfunctions Y 3 (s, x | z) are even under reflection x → −x, so that the coefficient functions in this expansion can be chosen even as well, η 3 (s, x, µ) = η 3 (s, −x, µ). They are characterized by two real numbers s > 0 and −∞ < x < ∞. It turns out that the corresponding anomalous dimensions can be written as a sum of terms depending on s and x separately. The s-dependent part can be absorbed in the same universal factor R(s; µ, µ 0 ) as for the leading twist so that one obtains [17] with the anomalous dimension [17] 3 Note that in addition to the integral over all real values of x the DA Φ 3 (z, µ) contains an extra contribution, the first term in Eq. (3.8), corresponding to a particular imaginary value x = i/2. This special term has a lower anomalous dimension separated by a finite number from the rest, continuum spectrum, and can be interpreted as the asymptotic DA. This interpretation fails, however, for large quark and/or gluon momenta ω 1 , ω 2 µ (alias small coordinates z 1 , z 2 1/µ) in which case contributions with all anomalous dimensions have to be included, see Ref. [17] for a detailed discussion. Note also that the twist-three contribution to the two-particle DA Φ − in Eq. (3.1) is determined entirely by this special term, η 3 (s, µ). The "genuine" three-particle twist-three contributions to Φ 3 encoded in η 3 (s, x, µ) decouple from Φ − to the stated 1/N 2 c accuracy. It is useful to have in mind that the evolution kernel H 3 is a hermitian operator so that its eigenfunctions Y 3 (s, x | z) are mutually orthogonal and form a complete set of functions with respect to a suitable scalar product. Explicit construction is given in Appendix B. In this way the coefficient functions η 3 (s, x, µ), η (0) 3 (s, µ) can be calculated as scalar products of the model DA with the corresponding eigenfunctions.
The renormalization group equations for twist-four DAs can be treated along the similar lines, the main difference being that one obtains a 2 × 2 matrix RG equation. The matrix structure is clear for the chiral-even case since the two existing chirality-even DAs Φ 4 and Ψ 4 + Ψ 4 are mixed by the evolution, and also for the chiral-odd case it is necessary to take into account additional operators containing transverse derivatives (A.3), see App. A. The complete solution of the evolution equations is presented in Eq. (A.29).
It turns out that this general result can be simplified under the assumption that twistfour four-particle B-meson DAs (with two gluon fields and/or with an extra quark-antiquark pair) are put to zero. Such DAs can be expected to be small and have autonomous scale dependence i.e. they do not mix with the three-particle DAs [29]. Thus they can be consistently put to zero at all scales and do not reappear via evolution.
With this assumption we are able to derive an exact relation between the two chiraleven DAs, and a similar relation for the chiral-odd DAs that allows one to eliminate contributions with transverse derivatives (A.32). The derivation is presented in App. C. In this way (to this accuracy) the general expressions for the scale dependence of the DAs in (A.29) reduce to and The first term in the expression for Ψ 4 − Ψ 4 (3.13c) can be interpreted as the Wandzura-Wilczek-type contribution of twist-three (related to the DA Φ 3 ) in the same manner as the twist-three two-particle DA Φ − contains a term related to the leading-twist DA Φ + , cf. Eq. (3.1). We can write The remaining "genuine" twist-four contributions are expressed in terms of two nonperturbative functions η Since the scale dependence is the same, one may wonder whether a relation between η exists on a nonperturbative level. Such relation is derived and discussed in the next section; its theoretical status is, however, less solid as compared to the rest of our results.
Going over to momentum space corresponds to a Fourier transform of the eigenfunctions so that, e.g., and similar for all twist-four DAs. 4 Explicit expressions for the eigenfunctions of the evolution equations in momentum space can easily be derived using that One obtains The RG kernels (for all twists) are hermitian operators with respect to the SL(2) scalar product, see App. A, B. As the result, eigenfunctions of the evolution equation corresponding to different anomalous dimensions are orthogonal and form a complete set of functions.
The resulting orthogonality relations are collected in App. B. They can be used to invert the representations in Eq. (3.13) and express the coefficients η 4 , κ 4 at a low reference scale in terms of the models for the DAs in momentum (or coordinate) space.

Two-particle higher-twist DAs and equations of motion
In the parton model language, higher twist effects are due to contributions of higher Fock states but also to nonvanishing parton transverse momenta (or virtuality). Due to QCD equations of motion (EOM) the latter can be expressed in terms of the multiparton configurations as well and can be thought of as contributions of gluon (or quark-antiquark pair) emission from the external lines of the partonic hard-scattering amplitude. In applications to hard exclusive reactions it has become customary to take into account these diagrams through the contributions of two-particle higher-twist DAs that arise as terms ∼ O(x 2 ) in the expansion of the relevant nonlocal quark-antiquark operator close to the light cone.
For the present case we can write, assuming introducing two new DAs, g + (ω) and g − (ω) that are of twist four and five, respectively. Terms O(x 4 ) are neglected. Eq. (4.1) has to be understood as a light-cone expansion to the tree-level accuracy which should be sufficient for the calculation of higher-twist corrections to the leading order in the strong coupling.
Note that the l.h.s. of Eq. (4.1) cannot have a power singularity 1/(vx) at (vx) → 0. This implies the constraints The DAs g + (ω) and g − (ω) can be expressed in terms of the three-particle DAs considered in previous sections. The corresponding expressions were derived by Kawamura et al. (KKQT) [21] starting from the operator identities taking the appropriate matrix elements and comparing the resulting expressions with the definition of the DAs in the limit x 2 → 0. In this way one obtains [21] The first KKQT relation, Eq. (4.4a), only involves twist-two and twist-three contributions. It allows to calculate the twist-three DA Φ − (z) in terms of Φ + (z) (Wandzura-Wilczek contribution [3]) and the "genuine" twist-three three-particle DA Φ 3 (z 1 , z 2 ). This relation can be derived in many ways (see e.g. App. C) and was used to arrive at the representation for Φ − (z) in Eqs. (3.1), (3.2) [17]. The second and the third relation, Eqs. (4.4b) and (4.4c), provide one with the expressions for the two-particle higher-twist DAs G ± (z) in terms of the three-particle DAs of the same twist and lower-twist Wandzura-Wilczek-type terms.
The last KKQT relation, Eq. (4.4d), is a nontrivial constraint relating the higher-twist matrix elements with the leading twist. Using the representation in Eqs. (3.1), (3.13) one obtains from Eq. (4.4d) the following relation for the coefficient functions: This equation presents a nonlocal generalization of the Grozin-Neubert relations [4] ∞ where λ 2 E and λ 2 H are matrix elements of certain local quark-gluon operators (5.1). It is easy to show that Eq. (4.8) correspond to the expansion of Eq. (4.7) at s → 0 and collecting terms O(s) and O(s 2 ), respectively.
Since γ 4 (x = 0) = 0 (3.18), the scale dependence of the higher-twist contributions on the r.h.s. of Eq. (4.7) matches the scale dependence of the leading-twist DA on the l.h.s. of this relation, however, only if the derivatives ∂ s are not applied to the s-dependent R-factor (3.5). This difficulty is due to the fact the light-cone expansion x 2 → 0 in (4.1) beyond tree level requires a careful treatment of the x 2 -dependent cusp anomalous dimension. A detailed investigation of this problem goes beyond the tasks of this work.
In order to tame potentially large corrections ∼ α s ln(sµ) to Eq. (4.7) one can try to enforce this relation for the integrated quantities, ∞ 0 ds, in which case it transforms into a constraint on the low-momentum behavior of the DAs that is most relevant for applications. In this way one obtains after a short calculation where φ ′ (ω, µ) = ∂ ω φ(ω, µ) and (ψ 4 ) tw−4 is the "genuine" twist-four contribution to the DA ψ 4 , In any case, it is important to have in mind that the expressions for the two-particle higher-twist DAs G ± (z) in (4.4b),(4.4c) are obtained in the same approximation as the constraints in (4.4d), (4.7), or, in the minimal version, (4.9). These constraints have to be fulfilled, for consistency, for any model of G ± (z) at a low scale.
In Ref. [21] a model for the leading-twist DA was formulated, called there "Wandzura-Wilczek approximation", by putting all quark-gluon contributions on the r.h.s. of Eqs. (4.7) to zero. We think that this approximation is not viable and the interpretation referring to Wandzura and Wilczek is misleading. 6 In this calculation one has to start with the regularized version of the integral 1 0 du Ψ4(z, uz) → 1 0 du u ǫ Ψ4(z, uz) and take the limit ǫ → 0 at the end.
Similar EOM constraints are familiar and widely used for the light-quark systems, see e.g. [30,31]. The simplest of them is the following operator identity for the divergence of the quark energy momentum tensor (for massless quarks): where Such identities exist for all leading twist operators. The general statement is that the divergence (in mathematical sense) of a multiplicatively renormalizable leading twist operator can be expressed as a sum of quark-gluon operators [32]. If all quark-gluon contributions are put to zero, one obtains an infinite series of conserved currents. This is in fact the symmetry of a naive parton model (free quarks), and indeed sending g → 0 is the only way to get rid of quark-gluon operators in a theoretically consistent way. A free theory is, however, not a viable approximation for modeling of the bound states. Note that the scale dependence must be neglected in this case as well. This approximation has been, therefore, never followed, to the best of our knowledge. Instead, the EOM relations of this type have commonly been used as constraints that allow one to reduce the number of independent twist-four matrix elements, see e.g. [30][31][32]. The situation with the heavy-light operators is analogous; Eq. (4.7) provides one with powerful constraints on the higher-twist DAs by reducing the number of nonperturbative parameters but does not imply any restrictions on the leading-twist DA itself.

Models
Modeling higher-twist DAs requires certain nonperturbative input which is currently very limited. Aim of this section is to present a few phenomenologically acceptable models that satisfy all known constraints. Matrix elements of local operators can be estimated from QCD sum rules. One obtains where in the second calculation some NLO corrections have been taken into account. Note that the ratio is almost the same in both cases and is generally more reliable than the values of the matrix elements themselves as many uncertainties cancel. Assuming that the integrals over quark and gluon momenta at a low scale converge, one obtains normalization conditions for the DAs 5) or, equivalently, Under the same assumptions, at a low scale

Model I: Exponential
The simplest model can be obtained combining the known low-momentum behavior (2.21) with an exponential falloff at large momenta, and using the above normalization conditions (cf. [11]): This construction is similar in spirit to the simple leading-twist DA proposed by Grozin and Neubert [4] φ + (ω, and has the advantage that all relevant coefficient functions in the expansion over the eigenfunctions of the evolution equations can be calculated explicitly in analytic form: The EOM relation (4.7) becomes For the simplest leading-twist DA in (5.9) this equation is satisfied if in agreement with (4.8). If these relations are enforced, there remains to be one free parameter, e.g., the ratio R = λ 2 E /λ 2 H (5.4). Using the above explicit expressions for the three-particle DAs we can compute the two-particle twist-four DA g + (ω) from Eq. (4.4b). Since this relation is derived in the same approximation as the EOM (4.4d), we have to require that Eq. (5.11) is satisfied identically. This means, e.g., that for the Grozin-Neubert leading-twist DA (5.9) the relations in (5.12) have to be enforced. One obtains where Ei(x) is the exponential integral. For small momenta 14) The function g + (ω) (5.13) is plotted in Fig. 2 where we have taken R = 1/2 (5.4). It is interesting that despite a rather elaborate analytic expression in (5.13) this function can be approximated with a very good accuracy by the simple expression Note that all EOM relations between the DAs are linear so that a more general model can be constructed as an arbitrary linear combination of the above expressions with different values of ω 0 . For the leading-twist DA containing a large-momentum "tail" φ + (ω) ∼ ln ω ω the definition of g + (ω) in Eq. (4.4b) and the EOM relation (5.11) both have to be modified. The problem is seen, e.g., using the expansions at µs ≪ 1 Using this expansion it is easy to check that the expression on the l.h.s. of Eq. (5.11) acquires terms ∼ α s sΛ ln(µs) that do not match the assumed O(s 2 ) behavior of the twistfour contributions. Before a more satisfying solution is available, we suggest to use the low-momentum part of φ + (ω) only for the construction of the higher-twist DA models at a low reference scale. Alternatively, the sensitivity to large-momentum contributions can be removed by imposing EOM for the integrated coefficient functions (in s-space), cf. (4.9). In this way one obtains the relation which can be used as a constraint on the twist-four parameters for a more general leadingtwist DA model. The scale dependence of the twist-four DAs ψ 4 +ψ 4 and φ 4 omitting the overall prefactor (λ 2 E + λ 2 E )/ω 2 0 is shown in the upper and the lower figure in Fig. 3, respectively. The DAs at the scale µ = 2.5 GeV, shown in yellow, are overlaid with the input expressions (5.8) where we assumed, for definiteness, µ 0 = 1 GeV. In this plot we use the variables so that ω is the total momentum of the light degrees of freedom and u is the fraction of the total momentum carried by the antiquark. Evolution has two effects. One of them is to suppress the higher-twist DAs at small momenta at larger scales and create the largemomentum "tails" similar to the leading-twist DA. Another effect is that the three-particle DAs at higher scales are tilted towards a larger momentum fraction carried by the gluon so that the region u → 0 is enhanced and the opposite region u → 1 depleted. This shift becomes more pronounced for larger values of the total momentum. This general pattern is seen in Fig. 3: the DAs at scale µ = 2.5 GeV are smaller than the input ones at µ = 1 GeV in the whole u region for ω 4ω 0 whereas for ω 4ω 0 the DAs at µ = 2.5 GeV are larger than at µ = 1 GeV for an increasingly broad interval in u. Effects of the scale dependence on the twist-three DA φ 3 are qualitatively similar and are discussed in detail in Ref. [17].

Model II: Local duality
Another class of models can be constructed using the local duality assumption, which is that the contribution of the B-meson state to the correlation functions of the type is equal to the perturbative spectral density integrated over a certain energy interval (interval of duality). Evaluating the simple quark loop and expanding in powers of x 2 one obtains in this way with ω 0 = (4/3)Λ. This "naive" model for g + (ω) can be used for a rough estimate of the size of the higher-twist effects. Selfconsistent models of this type can be constructed by including three-particle contributions. It turns out that the constraints due to the EOM (4.4d) require using the same threshold and power behavior ∼ (2ω 0 − ω) p for all DAs except for φ 3 which has to be one power lower. Two examples are given below.
Model IIA: Choosing p = 1 we obtain where in the expressions for three-particle DAsω = ω 1 + ω 2 . This set of DAs is consistent with EOM (4.4d), (4.7) provided the parameters satisfy the Grozin-Neubert constraints (4.8): The two-particle twist-4 DA g + (ω) in Eq. (5.23) can to a high accuracy be approximated by a simpler expression where R is defined in (5.4).
Model IIB: Choosing p = 3 we obtain instead and  where, as above, in the expressions for three-particle DAsω = ω 1 + ω 2 . This set of DAs is consistent with EOM (4.4d), (4.7) provided the parameters satisfy the constraints (4.8): 28) The two-particle twist-4 DA g + (ω) in Eq. (5.27) can to a high accuracy be approximated by a simpler expression where R is defined in (5.4). The four models for the two-particle twist-4 DA g + (ω) described in the text are compared to each other in Fig. 4. Note that the model parameter ω 0 has in each case to be adjusted to the same physical scaleΛ. By construction the four considered models have different high-energy behavior but their difference at small momenta proves to be rather moderate. This is encouraging and allows one to hope that the model uncertainties in the power-suppressed contributions to B-meson decay form factors can be kept under control.

Summary
Motivated by the challenge to control power-suppressed 1/m B contributions to B-decays in final states with energetic light particles in the framework of QCD factorization, in this work we present a systematic study of three-particle higher-twist DAs of the B-meson which are main nonperturbative input in such analysis.
We find eight independent three-particle DAs and classify them according to collinear twist and chirality which is related to their properties under conformal transformations. The twist-three three-particle DA φ 3 and the related two-particle DA φ − (ω) has been studied in detail in Ref. [17]. In this work we concentrate on the three existing twistfour DAs. Our principal result, Eq. (3.13), is the expansion of the DAs in terms of the eigenfunctions of the evolution equation and the calculation of the corresponding anomalous dimension, Eq. (3.18). We also derive a new relation between the two chiral-even DAs, Eq. (3.12), that is valid up to four-particle contributions of the typeqGGh v orqqqh v , and a similar relation for the chiral-odd DAs that allows one to eliminate contributions with extra transverse derivatives (A.32). We introduce two-particle higher-twist DAs that appear in the light-cone expansion of the heavy-light correlation functions following the approach of Kawamura et al. [21], and discuss the corresponding EOM relations. We introduce several simple models for the higher-twist DAs with different large-energy behavior that satisfy all tree-level EOM constraints. We find that if the constraints due to EOM are enforced, model dependence of the low-momentum part of the two-particle twist-four DA g + (ω) is rather mild, see Fig. 4. This is encouraging and allows one to hope that the twist-four uncertainties in the predictions for physical observables in B-decays can be kept under control.
Several issues are not covered in this work and require further studies. The first of them concerns the definition of two-particle twist-four (and higher) DAs that are auxiliary objects that arise via the light-cone expansion. The definition due to Kawamura et al. that is employed in our work is sufficient at the tree-level, but it has to be made more precise depending on the assumed power counting and particular factorization technique (e.g. SCET) beyond this accuracy. A related problem is that the EOM relation (4.4d) between the DAs must be modified by O(α s ) corrections as in the present form it is not consistent with the RG equations. It would be very important to find an alternative derivation of this relation that does not require considering off-light cone operators at the intermediate step.
Another problem is of course that the theory of 1/m B corrections to B-decays is at its infancy. It is well known that the contributions of higher-twist DAs are plagued by end-point divergences so that the twist hierarchy is lost, in general. In recent years there had been some progress in this direction based on using dispersion relations [11][12][13][14][15][16]. In this technique (light-cone sum rules) the hierarchy of higher-twist contributions is restored as an expansion in powers of the duality interval in the light hadron (photon) channel. Apart from (many) open theory issues, it remains to be seen, however, whether such techniques can provide one with phenomenologically relevant precision. Acknowledgments Y. J. is grateful to Henry Lamm for useful comments. The work of Y.J. and A.M. was supported by the DFG, grants BR 2021/7-1 and MO 1801/1-1, respectively.

Appendices A Renormalization group equations
The structure of the renormalization group equations for higher-twist operators is much simpler in the spinor representation. To this end we consider the nonlocal (light-ray) twist-three operator and two doublets of twist-four operators corresponding to the quark and gluon field of the same and opposite chirality Matrix elements of these operators define the B-meson DAs of interest, see Sect. 2.3. Note that we introduce an extra operator with a transverse derivative acting on the quark field. The corresponding matrix element defines a new DA, We will show that this DA can be expressed in terms of the other ones using equations of motion (EOM) so that the operator containing the transverse derivative can be eliminated. In this way, however, SL(2) symmetry of the evolution equation becomes obscured. It proves to be advantageous [25,26] to treat both operators as independent at the intermediate step and impose the EOM relations on the solutions. The light-ray operators satisfy renormalization group equations The evolution kernels H are integral operators acting on the light-cone coordinates of the fields [34] and, to one-loop accuracy, can be written as sums of two-particle kernels which describe the interaction between the partons, The and similar for H. Explicit expressions for the two-particle kernels are collected in Appendix D.
It is well known that evolution equations for the light quark and gluon operators enjoy SL(2, R) symmetry at one loop level. Each light field transforms according to a certain representation of the SL(2, R) group which is defined by one number, the conformal spin j. The conformal spins of "plus" and "minus" quark fields χ + , χ − are j = 1 and j = 1/2. The gluon fields f (f ) ++ and f (f ) +− transform according to the representation with j = 3/2 and j = 1, respectively, and the (holomorphic) transverse derivative of the quark field, D −+ χ + , has spin j = 3/2. The SL(2, R) (conformal) symmetry implies that the evolution kernels commute with the generators of symmetry transformations so that the two-particle light kernels can be written in terms of the corresponding quadratic Casimir operators [29]. This SL(2)-invariant representation (cf. App. D) is very convenient for the further analysis.
For the heavy-light operators the SL(2, R) symmetry breaks down. However, the evolution equations remain to be invariant with respect to the special conformal transformations [27]. Hence, e.g., the kernel H 3 commutes with 2 ) is a two-particle generator The corresponding generators for the twist-4 kernels are 2 × 2 matrices, , (A.10) Note that the conformal spins of up and down components in the twist-4 doublets (A.2) are different. This residual symmetry is sufficient to solve the evolution equation for the two-particle leading-twist DA [9], but is not enough for three-particle DAs which are functions of two variables. Fortunately it turns out that the leading large-N c contributions to the highertwist equations possess a hidden symmetry, so that we are able to construct another operator that commutes with the evolution kernel. For the twist-three case, this operator was found in Ref. [17]: The corresponding operators ("conserved charges") for the twist-4 kernels take the form (this is a new result): where , (A.15) . (A.16) The operators Q i , Q i , i = 1, 2 satisfy the commutation relations: The above expressions for Q 2 and Q 2 can be derived and the commutation relations verified most easily using the Quantum Inverse Scattering Method (QISM). This derivation will be presented elsewhere. With this construction, the problem in question becomes mathematically equivalent to a quantum-mechanical system with two degrees of freedom, with H playing the role of the Hamiltonian and Q 1 , Q 2 the conserved charges. Thanks to commutativity all three operators share the same set of the eigenfunctions. Thus, instead of trying to find the eigenfunctions of H directly, one can construct eigenfunctions of the charges Q 1 , Q 2 which are much simpler.
The complete set of twist-three eigenfunctions was obtained in [17], They are labeled by two quantum numbers s > 0 and x ∈ R related to the eigenvalues of the conserved charges, The charges Q i are self-adjoint operators w.r.t the SL(2)-invariant scalar product so that the eigenfunctions Y 3 (s, x | z) are mutually orthogonal, see App. B.
For the twist-four case, in addition to the two conserved charges Q 1 and Q 2 , one can construct an extra, "supplementary" charge Q 3 , such that [Q 1,2 , Q 3 ] = 0. The "supplementary" charges have a rather simple form .20) and are quite helpful for constructing of the eigenfunctions.
To this end we start with the following ansatz where j 1 and j 2 are the conformal spins of the light fields in the corresponding operator. Such functions are, by construction, eigenfunctions of the first charge Q 1 (Q 1 ): Requiring that the supplementary charge Q 3 (Q 3 ) is diagonalized yields a simple relation between the "upper" and "lower" components, ϕ 1 (u) and ϕ 2 (u), after which the eigenvalue problem for Q 2 (Q 2 ) reduces to a second order differential equation for, e.g., ϕ 1 (u). Selecting the solutions with the required analytic properties one obtains the following set of eigenfunctions: • Chiral operators (eigenfunctions of Q i -charges): Note that Z (+) 4 -functions are related to the twist-three eigenfunctions Y 3 as follows: 25) and the same relation is valid between Z (0) 4 and Y 3 . The last step is to calculate the eigenvalues of the Hamiltonians (anomalous dimensions): This can most easily be done by comparing the large-z 1 , z 2 asymptotic behavior on the both sides. In this way one obtains where γ 3 (x) and γ 4 (x) are defined in Eqs. (3.11) and (3.18), respectively. Expansion of the B-meson DAs over the eigenfunctions of the evolution equations reads where L = α s (µ)/α s (µ 0 ) and R(s; µ, µ 0 ) is defined in Eq. (3.5).
The expressions in (A.29) are valid for the most general case. We will show in App. C that neglecting contributions of four-particle (quasipartonic) [29] operators of the typē qGGh v andqqqh v , the following relations hold: The resulting simplified expressions for the DAs are given in Eq. (3.13) in the main text. We remind that quasipartonic four-particle operators do not mix with three-particle operators to the one-loop accuracy [29]. For this reason neglecting such contributions is consistent with the scale dependence.

B Conformal scalar product
Finding a suitable scalar product on the space of the B-meson DAs is an auxiliary, although very useful, tool for solving the evolution equations. The requirements for a scalar product are the following: i) the DAs under consideration have to belong to the corresponding Hilbert space ii) the evolution kernels have to be self-adjoint operators. It is not determined uniquely by the physics of the problem so that there is a certain freedom. The starting observation is that, first, support properties of the B-meson DAs in momentum space, ω > 0, imply that they are analytic functions of the light parton coordinates z in the lower half plane, Im z < 0 and, second, the evolution kernels at one loop order are functions of the SL(2, R) generators. These properties invite for using a formalism where z is treated as a complex number and conformal symmetry of the equations is implemented explicitly. Such a formalism is well known in mathematical literature.
One defines the SL(2, R) invariant scalar product for functions holomorphic in the lower complex half-plane [35] where the integration goes over the lower half-plane C − of the complex plane, Im z < 0, and the integration measure for spin j is defined as The scalar product (B.1) is invariant w.r.t. to the SL(2, R) transformations [35] The generator of special conformal transformation iS (j) + is self-adjoint w.r.t. this scalar product and its eigenfunctions iS (j) are orthogonal and form a complete set of functions in the Hilbert space 7 defined by Eq. (B.1): The expression on the r.h.s. of Eq. (B.5) is the reproducing kernel (unit operator) [36], i.e. for arbitrary function (holomorphic in the lower half-plane) Exponential functions e −iωz , ω > 0 form another complete orthogonal set w.r.t. the same scalar product, The momentum space DAs defined by the Fourier transform can be found making use of (B.7) and the following relation: In this way one obtains In particular for j = 1/2 corresponding to the B-meson DA φ − (ω, µ) the conformal expansion goes over Bessel functions J 0 (2 √ ωs) as compared to J 1 (2 √ ωs) for the leading twist, in which case j = 1.
The coefficient functions of DAs in the s-representation can be written as the scalar products as well, e.g. for the leading twist The SL(2) invariant scalar product for the functions of two variables reads, in position and momentum space representations, respectively. The scalar product for the twist-4 doublets, Φ(z) = (Φ 1 (z), Φ 2 (z)), can be written in the form where the two terms on the r.h.s. are defined by the scalar product (B.12) with the conformal spins matching those of the corresponding -up or down -component of the doublet. The coefficient c is fixed by the requirement that the Hamiltonian and conserved charges are self-adjoint operators. It is easiest to impose this condition on the "supplementary" charges, Q 3 (Q 3 ). For example, for the same-chirality doublet one has to require where from it follows that c = 2. In this way we obtain (for all possible combinations of superscripts ±, 0 in the bra-and ket-states) where the ellipses stand for contributions of quasipartonic twist-four four-particle DAs of the typeqGGh v orqqqh v . We use a technique that has been developed, for light quarks, in Refs. [25,37] and is based on the two-component spinor formalism. To this end it is convenient to relax the normalization condition Eq. (2.23) the auxiliary λ and µ spinors, so that they can be treated as independent.
As a simpler example, let us first re-derive in this approach the familiar relation between the twist-three DAs, Eq. (4.4a). To start with, note that the matrix element involving plus components of the heavy quark and the antichiral light quark vanishes identically: Physics reason is that the quark helicities do not combine to zero. Formally, this matrix element vanishes because is not possible to construct a tensor with two "undotted" spinor indices T αβ = 0|χ α (zn)[zn, 0]h β (0)|B such that T ++ = λ α λ β T αβ = 0 from the two vectors v αα and n αα = λ αλα at our disposal. Alternatively, this can be seen directly from the definition in Eq. (2.1).
The trick is to consider the derivative which, of course, must vanish as well. Since χ + = λ α χ α , h + = λ α h α and n αα = λ αλα , applying the derivative µ∂ λ we obtain three separate contributions that have to sum to zero: where taking the matrix element 0| . . . |B is implied. We stress that n is a light-like vector. Nevertheless, taking the derivative ∂/∂n ρ one can ignore the constraint n 2 = 0 and treat all four components of n µ as independent ones. Indeed, let F (n) be an arbitrary function of n µ defined on the surface n 2 = 0 and extend it formally to n 2 = 0, F (n) → F (n)| n 2 =0 + n 2 F 1 (n). 8 Then (µσ ρλ ) ∂ ∂n ρ n 2 F 1 (n) so that the contribution of the added term vanishes: the answer does not depend on F 1 (n) which means that the chosen extension beyond the light-cone surface does not matter. Note that we have to differentiate both the (anti)quark field and the Wilson line: The difference of the present approach to the one used in [21] is that here we start with a matrix element of the renormalized light-ray operator which is finite by definition, whereas the light-cone expansion x 2 → 0 in [21] produces singularities that have to be isolated in coefficient functions.
D Two-particle evolution kernels is the projector on the states with nonzero conformal spin.

D.2 SL(2)-invariant representation
It is also possible to present above integral forms of the evolution kernels in terms of SL(2) transformation generators. [