NNLO Vertex Corrections in charmless hadronic B decays: Imaginary part

We compute the imaginary part of the 2-loop vertex corrections in the QCD Factorization framework for hadronic two-body decays as B ->pi pi. This completes the NNLO calculation of the imaginary part of the topological tree amplitudes and represents an important step towards a NNLO prediction of direct CP asymmetries in QCD Factorization. Concerning the technical aspects, we find that soft and collinear infrared divergences cancel in the hard-scattering kernels which demonstrates factorization at the 2-loop order. All results are obtained analytically including the dependence on the charm quark mass. The numerical impact of the NNLO corrections is found to be significant, in particular they lead to an enhancement of the strong phase of the colour-suppressed tree amplitude.


Introduction
Charmless hadronic B decays provide important information on the unitarity triangle which may help to reveal the nature of flavour mixing and CP violation. In order to exploit the rich amount of data that is currently being collected at the B factories, a quantitative control of the underlying strong-interaction effects is highly desirable. In the QCD Factorization framework [1] the hadronic matrix elements of the operators in the effective weak Hamiltonian simplify considerably in the heavy-quark limit. Schematically, where the non-perturbative strong-interaction effects are encoded in a form factor F BM 1 + at q 2 = 0, decay constants f M and light-cone distribution amplitudes φ M . The shortdistance kernels T I i = O(1) and T II i = O(α s ) provide the basis for a systematic implementation of radiative corrections; the former contain the short-distance interactions that do not involve the spectator antiquark from the decayingB meson (vertex corrections) and the latter describe the ones with the spectator antiquark (spectator scattering).
The next-to-leading order (NLO) corrections to the kernels T I,II i , which constitute an O(α s ) correction to naive factorization, are already known from [1]. Recently, the nextto-next-to-leading order (NNLO) corrections to T II i have been computed [2,3,4,5,6]. Due to the interaction with the soft spectator antiquark, the spectator scattering term receives contributions from the hard scale ∼ m b and from an intermediate (hard-collinear) scale ∼ (Λ QCD m b ) 1/2 . Both types of contributions are now available at O(α 2 s ) (1-loop), indicating that the NNLO corrections are numerically important.
In this work we compute NNLO corrections to T I i for the so-called topological tree amplitudes (which arise from the insertion of current-current operators). In contrast to the spectator scattering term, the vertex corrections are dominated solely by hard effects and amount to a 2-loop calculation. In particular, we address the imaginary part of the hard-scattering kernels which is the origin of a strong rescattering phase shift that blurs the information on the weak phases. As an imaginary part is first generated at O(α s ), higher order perturbative corrections are expected to significantly influence the pattern of strong phases and hence direct CP asymmetries. Our calculation represents an important step towards a NNLO prediction of direct CP asymmetries in QCD Factorization.
The outline of this paper is as follows: In Section 2 we present our strategy for the calculation of the topological tree amplitudes by introducing two different operator bases. Section 3 contains the technical aspects of the 2-loop calculation. In Section 4 we show how to extract the hard-scattering kernels from the hadronic matrix elements. Our analytical results can be found in Section 5. The numerical impact of the NNLO vertex corrections is investigated in Section 6 and we finally conclude in Section 7. A more detailed presentation of the considered calculation can be found in [7].

Choice of operator basis
In view of the calculation of topological tree amplitudes, we restrict our attention to the current-current operators of the effective weak Hamiltonian for b → u transitions Due to the fact that we work within Dimensional Regularization (DR), we also have to consider evanescent operators [8]. These non-physical operators vanish in four dimensions but contribute at intermediate steps of the calculation in d = 4 − 2ε dimensions. As the imaginary part has effectively 1-loop complexity with respect to renormalization at O(α 2 s ), the considered calculation only requires 1-loop evanescent operators. For our purposes the complete operator basis is thus given bỹ where i, j are colour indices and L = 1 − γ 5 . The operator basis in (3) has been used in all previous calculations within QCD Factorization [1,2,3,4]. We refer to this basis as the traditional basis for convenience and denote the corresponding Wilson coefficients and operators with a tilde. It has been argued by Chetyrkin, Misiak and Münz (CMM) that one should use a different operator basis in order to perform multi-loop calculations [9]. Although the deeper reason is related to the penguin operators which we do not consider here, we prefer to introduce the CMM basis in view of future extensions of our work. This basis allows to consistently use DR with a naive anticommuting γ 5 to all orders in perturbation theory. In the CMM basis the current-current operators and corresponding 1-loop evanescent operators read (indicated by a hat) with colour matrices T A and colour indices i, j, k, l.
Comparing the operator bases in (3) and (4) we observe two differences: First, the two bases use different colour decompositions which is a rather trivial point. More importantly, they contain slightly different definitions of evanescent operators. Whereas the definitions in the CMM basis correspond to the simplest prescription to define evanescent operators, subleading terms of O(ε) appear in the one of the traditional basis. These terms have been properly adjusted such that Fierz symmetry holds to 1-loop order in d dimensions.
We follow the notation of [10] and express the hadronic matrix elements of the effective weak Hamiltonian in terms of topological amplitudes α i (M 1 M 2 ). E.g. the B − → π − π 0 decay amplitude is written as where  [10].
According to this definition, the left (right) diagram in Figure 1 contributes to the tree amplitude α 1 (α 2 ). On the technical level these two insertions of a four-quark operator correspond to two different calculations. Instead of performing both calculations explicitly, we may alternatively compute the amplitude α 2 by inserting Fierz reordered operators into the left diagram of Figure 1. To do so, it is essential to work with an operator basis that respects Fierz symmetry in d dimensions. As we have argued above, Fierz symmetry indeed holds to 1-loop order in the traditional basis from (3) which allows us to derive α 2 directly from α 1 .
We conclude that the CMM basis is the appropriate choice for a 2-loop calculation whereas the traditional basis provides a short-cut for the derivation of the coloursuppressed amplitude. We therefore pursue the following strategy for the calculation of the imaginary part of the NNLO vertex corrections: We perform the explicit 2-loop calculation in the CMM basis using the first type of operator insertion in the left diagram from Figure 1. From this we obtain α 1 (Ĉ i ). We then transform this expression into the traditional basis which yields α 1 (C i ) and finally apply Fierz symmetry arguments to derive α 2 (C i ) from α 1 (C i ) by simply exchangingC 1 ↔C 2 .

2-loop calculation
The core of the considered calculation consists in the computation of the matrix elements Q 1,2 ≡ u(p ′ )d(q 1 )ū(q 2 )|Q 1,2 |b(p) (6) to O(α 2 s ) which represents a 2-loop calculation. As will be described in Section 4.2, only (naively) non-factorizable diagrams with at least one gluon connecting the two currents in the left diagram of Figure 1 have to be considered here. The full NNLO calculation thus involves the 2-loop diagrams depicted in Figure 2, but only about half of the diagrams give rise to an imaginary part. It is an easy task to identify this subset of diagrams since the generation of an imaginary part is always related to final state interactions.
In our calculation we treat the partons on-shell and write q 1 = uq, q 2 =ūq and p ′ = p − q satisfying q 2 = 0 and p 2 = 2p · q = m 2 b (withū ≡ 1 − u). We use DR for the regularization of ultraviolet (UV) and infrared (IR) divergences and an anticommuting γ 5 according to the NDR scheme. We stress that we do not perform any projection onto the bound states in the partonic calculation. We instead treat the two currents in the four-quark operators independently and make use of the equations of motion in order to simplify the Dirac structures of the diagrams. In order to calculate the large number of 2-loop integrals we proceed as follows: Using a general tensor decomposition of the loop integrals, we essentially deal with the calculation of scalar integrals. With the help of an automatized reduction algorithm, we are able to express several thousands of scalar integrals in terms of a small set of so-called Master Integrals (MIs). The most difficult part finally consists in the calculation of these MIs.
In the remainder of this section we present some techniques that we have found useful for the considered calculation. We sketch the basic ideas of the aforementioned reduction algorithm and discuss several techniques for the calculation of the MIs. We refer to the references quoted in the following subsections for more detailed descriptions (see also [7]).

Reduction to Master Integrals
Any scalar 2-loop integral in our calculation can be expressed as where the S i are scalar products of a loop momentum with an external momentum or of two loop momenta. The P i denote the denominators of propagators and the exponents fulfil n i , m i ≥ 0. The scalar integrals themselves depend on the convolution variable u in the factorization formula (1). Very few integrals, arising from diagrams with a charm quark in a closed fermion loop, depend in addition on the ratio z ≡ m c /m b . We have suppressed this dependence in (7) for simplicity. Notice that an integral has different representations in terms of {S, P, n, m} because of the freedom to shift loop momenta in DR. It is the underlying topology, i.e. the interconnection of propagators and external momenta, which uniquely defines the integral. In the following we loosely use the word topology in order to classify the integrals. An integral with t different propagators P i with n i > 0 is called a t-topology. The integrals in the considered calculation have topology t ≤ 6.
The reduction algorithm makes use of various identities which relate integrals with different exponents {n, m}. The most important class of identities are the integrationby-parts identities [11] which follow from the fact that surface terms vanish in DR In order to obtain scalar identities we may contract (8) with any loop or external momentum under the integral before performing the derivative. In our case of two loop and two external momenta we generate in this way eight identities from each integral. A second class of identities, called Lorentz-invariance identities [12], exploits the fact that the integrals in (7) transform as scalars under a Lorentz-transformation of the external momenta. In this way we may generate up to six identities from each integral depending on the number of external momenta. In our example with only two linearly independent external momenta p and q there is only one such identity given by In total we obtain nine identities from a given integral, each of the identities containing the integral itself, simpler integrals with smaller {n, m} and more complicated integrals with larger {n, m}. It is important that the number of identities grows faster than the number of unknown integrals for increasing {n, m}. Hence, for large enough {n, m} the system of equations becomes (apparently) overconstrained and can be used to express more complicated integrals in terms of simpler ones. Not all of the identities being linearly independent, some integrals turn out to be irreducible to which we refer as MIs.
In the considered calculation we typically deal with systems of equations made of several thousands equations. The solution being straight-forward, the runtime of the reduction algorithm depends strongly on the order in which the equations are solved. As a guideline for an efficient implementation we have followed the algorithm from [13].
The reduction algorithm enables us to express the diagrams of Figure 2 as linear combinations of MIs which are multiplied by some Dirac structures. As the coefficients in these linear combinations are real, we may extract the imaginary part of a diagram at the level of the MIs which is a much simpler task than for the full diagrams. As depicted in Figure 3, we find 14 MIs that contribute to the calculation of the imaginary part of the NNLO vertex corrections.

Calculation of Master Integrals
Some MIs in Figure 3 can be solved easily e.g. with the help of Feynman parameters. For the more complicated MIs the method of differential equations [14] in combination with the formalism of Harmonic Polylogarithms (HPLs) [15] turned out to be very useful. In this section we give brief reviews of these techniques and conclude with a comment on the calculation of the boundary conditions to the differential equations.
The analytical results for the MIs from Figure 3 can be found in Appendix A.1 of [7]. As an independent check of our results we evaluated the MIs numerically using the method of sector decomposition [16].

Method of Differential Equations
The MIs are functions of the physical scales of the process which are given by scalar products of the external momenta and masses of the particles. In our calculation the MIs depend on the dimensionless variable u as in (7).
For a given MI we perform the derivative with respect to u and interchange the order of integration and derivation The right-hand side being of the same type as equations (8) and (9), this procedure again leads to a sum of various integrals with different exponents {n, m}. With the help of the reduction algorithm, these integrals can be expressed in terms of MIs which yields a differential equation of the form where we indicated that the coefficients a and b j may depend on the dimension d. The inhomogeneity of the differential equation typically contains MIs of subtopologies which are supposed to be known in a bottom-up approach. In case of the MIs from the third line of Figure 3, one MI in the inhomogeneous part is of the same topology as the MI on the left hand side of (11) and thus unknown. Writing down the differential equation for this MI, we find that we are left with a coupled system of linear, first order differential equations.
We are looking for a solution of the differential equation in terms of an expansion Expanding (11) then gives much simpler differential equations for the coefficients c ij which can be solved order by order in ε. In addition, in the case where we were left with a coupled system of differential equations, the system turns out to decouple in the expansion. The solution of the homogeneous equations is in general straight-forward. The inhomogeneous equations can then be addressed with the method of the variation of the constant. This in turn leads to indefinite integrals over the inhomogeneities which typically contain products of rational functions with logarithms or related functions as dilogarithms. With the help of the formalism of HPLs these integrations simplify substantially.

Harmonic Polylogarithms
The formalism of HPLs [15] allows to rewrite the integrations mentioned at the end of the last section in terms of familiar transcendental functions which are defined by repeated integration over a set of basic functions. We briefly summarize their basic features here, focussing on the properties that are relevant for our calculation. The HPLs, denoted by H( m w ; x), are described by a w-dimensional vector m w of parameters and by its argument x. We restrict our attention to the parameters 0 and 1 in the following. The basic definitions of the HPLs are for weight w = 1 and for weight w > 1 where the basic functions f (a; x) are given by In the case of m w = 0, the definition in (14) does not apply and the HPLs read The HPLs form a closed and linearly independent set under integrations over the basic functions f (a; x) and fulfil an algebra such that a product of two HPLs of weight w 1 and w 2 gives a linear combination of HPLs of weight w = w 1 + w 2 . As described above, the solution of the differential equations leads to integrals over products of some rational functions with some transcendental functions as logarithms or dilogarithms. More precisely, we find e.g. integrals of the type It turns out that all these integrals can be expressed as linear combinations of HPLs of weight w + 1. This is obvious for the first integral as it simply corresponds to the definition of a HPL, cf. (14) with a = 1. For the other integrals in (17), an integrationby-parts leads either to a recurrence relation or again to integrals of the type (14). Not all integrals in our calculation fall into the simple pattern (17), but a large part of this calculation can be performed along these lines.
In the considered calculation we encounter HPLs of weight w ≤ 3. Our results can be expressed in terms of the following minimal set of HPLs .
The situation is more complicated for the last two MIs in the third line of Figure 3 where the internal charm quark introduces a new scale. However, a closer look reveals that these MIs depend on two physical scales only, namely um 2 b and m 2 c = z 2 m 2 b . The MIs can then be solved within the formalism of HPLs in terms of the ratio ξ ≡ z 2 /u if we allow for more complicated arguments of the HPLs as e.g. η ≡ 1 2 1 − √ 1 + 4ξ .

Boundary conditions
A unique solution of a differential equation requires the knowledge of its boundary conditions. In the considered calculation the boundary conditions typically represent singlescale integrals corresponding to u = 0 or 1. It is of crucial importance that the integral has a smooth limit at the chosen point such that setting u = 0 or u = 1 does not modify the divergence structure introduced in (12). In some cases the methods described so far can also be applied for the calculation of the boundary conditions since setting u = 0 or 1 leads to simpler topologies which may turn out to be reducible. If so, the reduction algorithm can be used to express the integral in terms of known MIs. If not, a different strategy is mandatory. In this case we tried to calculate the integral with the help of Feynman parameters and managed in some cases to express the integral in terms of hypergeometric functions which we could expand in ε with the help of the Mathematica package HypExp [17]. Finally, the most difficult single-scale integrals could be calculated with Mellin-Barnes techniques [18].

Renormalization and IR subtractions
The matrix elements Q i which we obtained from computing the 2-loop diagrams in Figure 2 are ultraviolet (UV) and infrared (IR) divergent. In this section we show how to extract the finite hard-scattering kernels T I i from these matrix elements.

Renormalization
The renormalization procedure involves standard QCD counterterms, which amount to the calculation of various 1-loop diagrams, as well as counterterms from the effective Hamiltonian. We write the renormalized matrix elements as where contains the wave-function renormalization factors Z b of the b-quark and Z q of the massless quarks, whereasẐ is the operator renormalization matrix in the effective theory. We introduce the following notation for the perturbative expansions of these quantities ren/bare , and rewrite (19) in perturbation theory up to NNLO which yields The full calculation thus requires the operator renormalization matricesẐ (1,2) . For the calculation of the imaginary part, the terms proportional to the tree level matrix elements do not contribute andẐ (2) drops out in (21) as expected for an effective 1-loop calculation. Mass and wave function renormalization are found to be higher order effects. For the renormalization of the coupling constant we use according to the MS-scheme. The expression for the 1-loop renormalization matrixẐ (1) can be found e.g. in [19] and readŝ 6 0 1 0 where the two lines correspond to the basis of physical operators {Q 1 ,Q 2 } and the four columns to the extended basis {Q 1 ,Q 2 ,Ê 1 ,Ê 2 } including the evanescent operatorsÊ 1,2 defined in (4).

Factorization in NNLO
In this section it will be convenient to introduce the following short-hand notation for the factorization formula (1) where F denotes the B → M 1 form factor, T i the hard-scattering kernels T I i and Φ the product of the decay constant f M 2 and the distribution amplitude φ M 2 . The convolution in (1) has been represented by the symbol ⊗ and the ellipsis contain the terms from spectator scattering which we disregard in the following.
Formally, we may introduce the perturbative expansions Up to NNLO the expansion of (24) then yields In LO the comparison of (21) and (26) gives the trivial relation which states that the LO kernels T (0) i can be computed from the tree level diagram in Figure 4a. In order to address higher order terms we split the matrix elements into its (naively) factorizable (f) and non-factorizable (nf) contributions The corresponding 1-loop diagrams are shown in Figure 4b and 4c respectively. To this order (21) and (26) lead to which splits into for the calculation of the NLO kernels T (1) which shows that the factorizable diagrams and the wave-function renormalization are absorbed by the form factor and wave function corrections F (1) and Φ (1) . This suggests in NNLO the following structure These terms are thus irrelevant for the calculation of the NNLO kernels T (2) i which justifies that we could restrict our attention to the non-factorizable 2-loop diagrams from Figure 2. In NNLO the remaining terms from (21) and (26) contain non-trivial IR subtractions This equation can be simplified further when we make the wave function renormalization factors in the form factor and the distribution amplitude explicit Notice that the resulting amputated form factor F amp and wave function Φ amp contain UVdivergences by construction. Using (30), we see that the wave function renormalization factors cancel and arrive at the final formula for the calculation of the NNLO kernels T As the tree level matrix elements and the factorizable 1-loop diagrams do not give rise to an imaginary part, these terms can be disregarded in the present calculation.

IR subtractions
We now consider the IR subtractions on the right hand side of (35) in some detail. Let us first address the NLO kernels T (1) i which can be determined from equation (30). The renormalization in the evanescent sector implies that the left hand side of (30) is free of contributions from evanescent operators up to the finite order ε 0 . However, as the NLO kernels enter (35) in combination with the form factor correction F (1) amp which contains double (soft and collinear) IR divergences, the NLO kernels are required here up to O(ε 2 ). Concerning the subleading terms of O(ε), the evanescent operators do not drop out on the left hand side of (30) and we therefore have to extend the factorization formula on the right hand side to include these evanescent structures as well. Schematically, with a kernel T i,E = O(1) and an evanescent tree level matrix element F Similarly, the right hand side of (35) has to be modified to include these evanescent structures Notice that the term with the kernel T i,E = O(1) is not required to extract the finite piece of the physical NNLO kernel T (2) i . From the calculation of the 1-loop diagrams in Figure 4c, we find that the NLO kernels vanish in the colour-singlet case, T 2,E = 0, whereas the imaginary part of the colour-octet kernels is given by , † In the notation of [2], the right hand side of (36) corresponds to matrix elements of SCET I operators of the form [(ξW c1 )Γ 1 h v ][(χW c2 )Γ 2 (W † c2 χ)]. In NNLO we match onto two SCET I operators with Diracstructures Γ 1 ⊗ Γ 2 given by O 1 = n / + L ⊗ n /− 2 L and O 2 = n / + γ µ ⊥ γ ν ⊥ L ⊗ n /− 2 γ ⊥µ γ ⊥ν L (in our notation p ′ = 1 2 m b n − and q = 1 2 m b n + ). The matrix element of O 1 defines our structure F (0) Φ (0) and the evanescent combination where L ≡ ln µ 2 /m 2 b and we recall thatū ≡ 1 − u.

Form factor subtractions
We now address the form factor corrections which require the calculation of the diagram in Figure 5 (for on-shell quarks) and its counterterm. According to the definition of F amp in (34), we do not have to consider the wave function renormalization of the quark fields.
We again have to compute the corrections for physical and evanescent operators. Concerning the physical operator, the counterterm is found to vanish which reflects the conservation of the vector current. The evaluation of the 1-loop diagram in Figure 5 gives which contains double IR singularities as mentioned at the beginning of this section. On the other hand, the 1-loop diagram with an insertion of the evanescent operator yields a contribution proportional to the evanescent and the physical operators. We now have to adjust the counterterm such that the renormalized (IR-finite) matrix element of the evanescent operator vanishes (which ensures that the evanescent structures disappear in the final factorization formula). We obtain The form factor subtractions in (37) then follow from combining (38), (39) and (40). We emphasize that the corrections related to the evanescent operator do not induce a contribution to the physical NNLO kernel T amp,E Im T

Wave function subtractions
Concerning the wave function corrections we are left with the calculation of the diagrams in Figure 6 for collinear and on-shell partons with momenta uq andūq. However, as in our set-up q 2 = 0 all these diagrams vanish due to scaleless integrals in DR. We conclude that the wave function corrections are determined entirely by the counterterms. We compute these counterterms by calculating the diagrams from Figure 6 with an (IR-finite) off-shell regularization prescription in order to isolate the UV-divergences (which are independent of the IR regulator). The counterterm for the physical operator is found to be with the familiar Efremov-Radyushkin-Brodsky-Lepage (ERBL) kernel [20] V . For the evanescent operator we obtain where V E (u, w) denotes the spin-dependent part of the ERBL kernel given by Notice that the evanescent operators do induce a finite contribution to the physical NNLO kernel T (2) 1 in this case as the convolution with the corresponding NLO kernel implies We finally quote the result for the convolution with the physical NLO kernel

Vertex Corrections in NNLO
We now have assembled all pieces required for the calculation of the NNLO kernels T (2) i from (37). We indeed observe that all UV and IR divergences cancel in the hard-scattering kernels as predicted by the QCD Factorization framework. Since this is the result of a complicated subtraction procedure, this can also be seen as a very stringent cross-check of our calculation.

α 1 in CMM basis
The procedure outlined so far leads to the colour-allowed tree amplitude in the CMM operator basis defined in (4). We write where the ellipsis denote the terms from spectator scattering which are irrelevant for our purposes. In the CMM basis, the imaginary part of the vertex correctionsV (1,2) takes the form In writing (49) we have made the dependence on the renormalization scale explicit and disentangled contributions that belong to different colour structures. The NLO kernel g 1 (u) is given by and the NNLO kernels h i (u) will be specified below. The kernel h 4 (u; z f ) stems from diagrams with a closed fermion loop and depends on the mass of the internal quark through z f = m f /m b . We keep a non-zero charm quark mass and write z ≡ z c = m c /m b for simplicity.
The NNLO kernels were so far unknown. They are found in this work to be h 0 (u) = 155 4 + 4ζ 3 + 4Li 3 (u) − 4S 1,2 (u) − 12 ln u Li 2 (u) + 4 3 ln 3 u − 6 ln 2 u lnū where the last kernel has been given in terms of In the massless limit h 4 (u; z) becomes

α 1 and α 2 in traditional basis
Following our strategy from Section 2, we compute the colour-suppressed amplitude α 2 by rewriting the colour-allowed amplitude α 1 in the traditional operator basis from (3). Manifest Fierz symmetry in this basis relates the two amplitudes via where the upper (lower) signs apply for i = 1 (i = 2) and the ellipsis correspond to the terms from spectator scattering. In order to deriveṼ (1,2) we have to transform the Wilson coefficients in the CMM basisĈ i into the ones of the traditional basisC i . To NLL approximation this transformation can be found e.g. in [9] and readŝ Combining (48), (54) and (55) we obtain 1 π ImṼ (1) = 1 π ImV (1) Equation (56) represents the central result of our analysis, specifying the imaginary part of the colour-allowed tree amplitude α 1 and the colour-suppressed tree amplitude α 2 according to (54). The expression forṼ (1) is in agreement with [1], whereas the expressions forṼ 1,2 are new. The kernels g 1 (u) and h i (u) can be found in (50) and (51), respectively. The terms proportional to n f have already been considered in the analysis of the large β 0 -limit in [21]. Our results are in agreement with these findings.

Convolutions in Gegenbauer expansion
Our results in (56) have been given in terms of convolutions of hard-scattering kernels with the light-cone distribution amplitude of the emitted meson M 2 . We may explicitly calculate these convolution integrals by expanding the distribution amplitude into the eigenfunctions of the 1-loop evolution kernel where a M 2 n and C (3/2) n are the Gegenbauer moments and polynomials, respectively. We truncate the Gegenbauer expansion at n = 2 and perform the convolution integrals analytically. We find − 2 1 − (2y 1 + 1)(1 + 22z 2 ) ln y 1 y 2 − 4 ln y 2 + 7 + 164z 2 + 180z 4 + 144z 6 − 288z 4 F (z) + 12z 4 (3 + 16z 2 + 12z 4 ) ln 2 y 1 y 2 − 2 1 − (2y 1 + 1)(1 + 22z 2 + 84z 4 + 72z 6 ) ln where we defined F (z) ≡ Li 3 (−y 1 ) − S 1,2 (−y 1 ) − ln y 1 Li 2 (−y 1 ) + 1 2 ln y 1 ln 2 y 2 − 1 12 In the massless limit the function H 4 (z f ) simplifies to The finiteness of the convolution integrals in (58) completes the explicit factorization proof of the imaginary part of the NNLO vertex corrections. We summarize our results for the vertex corrections in the Gegenbauer representation of the light-cone distribution amplitude of the meson We thus find large coefficients for the NNLO vertex corrections and expect only a minor impact of the higher Gegenbauer moments, in particular in the symmetric case with a M 2 1 = 0. Notice that all contributions add constructively in α 1,2 due to the relative signs of the Wilson coefficients,C 1 ∼ 1.1 andC 2 ∼ −0.2. In the case of α 1 the contribution fromṼ (2) 1 is found to exceed the formally leading contributionṼ (1) due to the fact that the latter is multiplied by the small Wilson coefficientC 2 . Concerning α 2 the NNLO vertex corrections are also substantial, roughly saying they amount to a 50% correction. In both cases the impact of the charm quark mass is small, we find a correction of ∼ 3% compared to the massless case. A more detailed numerical analysis including the contributions from spectator scattering will be given in the following section.
We finally remark that the large β 0 -limit considered in [21] fails to reproduce the imaginary part of α 1 as it completely misses the leading contribution fromṼ (2) 1 . In the case of α 2 the approximation turns out to be reasonably good with a deviation of ∼ 10% compared to the full NNLO result. 6 Numerical analysis 6

.1 Implementation of Spectator Scattering
In the numerical analysis we combine our results with the NNLO corrections from 1loop spectator scattering obtained in [2,3,4,5,6]. In contrast to the vertex corrections considered in this work, the spectator term receives contributions from the hard scale µ h ∼ m b and the hard-collinear scale µ hc ∼ (Λ QCD m b ) 1/2 . According to this, the kernels T II i from (1) factorize into hard functions H II i and a (real) hard-collinear jet-function J || . Evaluating both kernels at the same scale µ would imply parametrically large logarithms which may spoil the convergence of the perturbative expansion.
In order to resum these logarithms we follow reference [2] and perform the substitution where U || = e −S U || consists of a universal Sudakov factor S and a non-local evolution kernel U || . As an imaginary part is first generated at O(α 2 s ) in the spectator term, we implement the resummation in the leading-logarithmic (LL) approximation. In the traditional operator basis from (3) the respective imaginary part takes the form where we made the scale dependence of the parameters explicit and introduced the first inverse moment of the B meson light-cone distribution amplitude λ −1 B . We further wrote 1 = a M 0 (µ) in order to simplify the notation. In (64) the resummation is encoded iñ with the Sudakov factor S given in LL approximation in equation (106) of [6] and the kernels r i in equations (38) and (39) of [2]. Following [4] we defined which can be computed by solving numerically the integro-differential equation with initial condition C m (z; µ hc , µ hc ) = 6z C (3/2) m (2z − 1) and γ || from equation (99) of [6].
In order to illustrate the numerical importance of the resummation we compare the values of the imaginary part ofR mn i for m, n ≤ 2 and µ h = µ hc (line I, without resummation) and µ h = 4.8GeV, µ hc = 1.5GeV (line II, with resummation) in Table 1. We observe that the resummation leads to a suppression of the spectator term of ∼ 10% due to the universal Sudakov factor (e −S ≃ 0.89 for our choice of input parameters). The resummation effects induced by U || turn out to be of minor numerical importance.
According to (64) we must evolve the Gegenbauer moments of the mesons M 1 and M 2 to the hard-collinear and the hard scale, respectively. In LL approximation the Gegenbauer moments do not mix and the evolution reads with anomalous dimensions γ 1 = − 64 9 and γ 2 = − 100 9 . We are left with the evolution of the B meson parameters to the hard-collinear scale. We convert the HQET decay constantf B into the physical one f B using the LL relation The evolution of λ B is more complicated. The solution of the integro-differential equation, which governs the LL evolution of the B meson light-cone distribution amplitude, can be found in [22]. Here we adopt a model-description for the B meson distribution amplitude to generate the evolution of λ B . We take the model from [22] which has the correct asymptotic behaviour and is almost form-invariant under the evolution. Finally we implement the BBNS model from [1] in order to estimate the size of power corrections to the factorization formula (1). This results in an additional contribution from spectator scattering related to subleading projections on the light-cone distribution amplitudes of the light mesons. It is given by  where r M χ (µ) = 2m 2 M /m b (µ)/(m q +mq)(µ), ∆ M = 1+ n (−1) n a M n and X H parameterizes an endpoint-divergent convolution integral. The latter is written as which may generate an imaginary part due to soft rescattering of the final state mesons. We take ln m B /Λ h ≃ 2.3, ρ H = 1 and allow for an arbitrary phase ϕ H .

Tree amplitudes in NNLO
The numerical implementation of the vertex corrections from (54) is easier since they can be evaluated at the hard scale µ h ∼ m b and depend only on few parameters. As can be read off from (61), the first Gegenbauer moment of the meson M 2 enters the leading termṼ (1) . We therefore require its next-to-leading logarithmic (NLL) evolution (which can be found in [23]) but as we restrict our attention to B → ππ decays in the following the first moment does not contribute at all. Since the second Gegenbauer moment does not enterṼ (1) , it is only required in LL approximation given by (68). Our input parameters for the B → ππ tree amplitudes are summarized in Table 2. The value for the B meson decay constant is supported by QCD sum rule calculations [24] and recent lattice results [25]. The form factor F Bπ + (at large recoil) has been addressed in the light-cone sum rule (LCSR) approach [26]. As we implement the model from [22] for the B meson distribution amplitude, we take the respective value for λ B . This value is somewhat larger than the one used in previous QCD Factorization analyses [1,2,3,4], but it is supported by a QCD sum rule and a LCSR calculation [27]. The value for the second Gegenbauer moment of the pion can be inferred from a LCSR analysis [28] and lattice results [29].
In order to estimate the size of higher-order perturbative corrections we vary the hard scale in the range µ h = 4.8 +4.8 −2.4 GeV and the hard-collinear scale independently between µ hc = 1.5 +0.9 −0.5 GeV. Throughout we use 2-loop running of α s with n f = 5 (n f = 4) for quantities that are evaluated at the hard scale µ h (hard-collinear scale µ hc ). The quark masses are interpreted as pole masses except for those entering r M χ .
In these expressions we disentangled the contributions from the NLO (1-loop) vertex corrections V (1) , NNLO (2-loop) vertex corrections V (2) and NNLO (1-loop) spectator scattering S (2) . In both cases the NNLO corrections are found to be important, although small in absolute terms. For the imaginary part of α 1 the NNLO corrections exceed the formally leading NLO result which can be explained by the fact that the latter is multiplied by the small Wilson coefficientC 2 , cf. the discussion after (62). We further observe that the individual NNLO corrections come with opposite signs which leads to a partial cancellation in their sum. The phenomenologically most important consequence of our calculation may be the enhancement of the imaginary part of the colour-suppressed tree amplitude α 2 . In our error estimate in (72) we distinguished between uncertainties which originate from the variation of the hard and the hard-collinear scale (scale), from the variation of the input parameters in Table 2 (param) and from the BBNS model which we used to estimate the size of power-corrections (power). By now we expect the power corrections to be the main limiting factor for an accurate determination of the amplitudes. However, although the inclusion of NNLO corrections has reduced the dependence on the renormalization scales, it still remains sizeable (in particular for µ h ). For our final error estimate in the third line of each amplitude, we added all uncertainties in quadrature.
Finally we remark that our numerical values for the NNLO spectator terms are much smaller than the ones quoted in [2]. This is partly related to the fact that we use different hadronic input parameters, in particular a much larger value for λ B . In addition to this, the authors of [2] essentially evaluate the hard functions H II i at the hard-collinear scale in order to partly implement the (unknown) NLL resummation of parametrically large logarithms. The NLL approximation is indeed required for the real part of the amplitudes, but as long as we concentrate on the imaginary part it is consistent to work in the LL approximation as discussed in Section 6.1. As the spectator term is the main source for the uncertainties from the hadronic input parameters, we also obtain smaller error bars than [2].

Conclusion
We computed the imaginary part of the 2-loop vertex corrections to the topological tree amplitudes in charmless hadronic B decays. Together with the 1-loop spectator scattering contributions considered in [2,3,4,5,6], the imaginary part of the tree amplitudes is now completely determined at NNLO in QCD Factorization.
Among the technical issues we showed that soft and collinear infrared divergences cancel in the hard-scattering kernels and that the resulting convolutions are finite, which demonstrates factorization at the 2-loop order. In our numerical analysis we found that the NNLO corrections are significant, in particular they enhance the strong phase of the colour-suppressed tree amplitude α 2 . Further improvements of the calculation still require a better understanding of power corrections to the factorization formula.
Our calculation represents an important step towards a NNLO prediction of direct CP asymmetries in QCD Factorization. As the topological penguin amplitudes, which also affect the direct CP asymmetries, have not yet been computed completely at NNLO (the contribution from spectator scattering can be found in [4]), we refrained from discussing them already in this work. Moreover, the strategy outlined in this work may also be applied for the calculation of the real part of the topological tree amplitudes, which is, however, technically more involved [7,30].