Mueller Navelet jets at LHC - complete NLL BFKL calculation

We calculate cross section and azimuthal decorrellation of Mueller Navelet jets at the LHC in the complete next-lo-leading order BFKL framework, i.e. including next-to-leading corrections to the Green's function as well as next-to-leading corrections to the Mueller Navelet vertices. The obtained results for standard observables proposed for studies of Mueller Navelet jets show that both sources of corrections are of equal, big importance for final magnitude and final behavior of observables. The astonishing conclusion of our analysis is that the observables obtained within the complete next-lo-leading order BFKL framework of the present paper are quite similar to the same observables obtained within next-to-leading logarithm DGLAP type treatment. This fact sheds doubts on general belief that the studies of Mueller Navelet jets at the LHC will lead to clear discrimination between the BFKL and the DGLAP dynamics.


Introduction
The large center of mass energy of hadron colliders like the Tevatron and the Large Hadron Collider (LHC) is not only of interest for the production of possible new heavy particles, but also allows to investigate the high energy regime of Quantum Chromodynamics (QCD). An especially interesting situation appears if two different large scales enter the game. If the two scales are ordered, large logarithms of the ratio of the two scales compensate the smallness of the coupling and therefore have to be resummed to all orders. One famous example is the case of high energy scattering with fixed momentum transfer. If the center of mass energy s is much larger than the momentum transfer |t| -the so-called Regge asymptotics of the process -the gluon exchange in the crossed channel dominates and logarithms of the type [α s ln(s/|t|)] n have to be resummed. This is realized by the leading logarithmic (LL) Balitsky-Fadin-Kuraev-Lipatov (BFKL) [1][2][3][4]   One of the most famous testing ground for BFKL physics are the Mueller Navelet jets [5], illustrated in Fig. 1. The predicted power like rise of the cross section with increasing energy has been observed at the Tevatron pp-collider [6], but the measurements revealed an even stronger rise than predicted by BFKL calculations. Beside the cross section also a more exclusive observable within this process drew the attention, namely the azimuthal correlation between these jets. Considering hadron-hadron scattering in the common parton model to describe two jet production at LO, one deals with a back-to-back reaction and expects the azimuthal angles of the two jets always to be π and hence completely correlated. This corresponds in Fig. 1 to φ J,1 = φ J,2 − π . But when we increase the rapidity difference between these jets, the phase space allows for more and more emissions leading to an angular decorrelation between the jets. In the academical limit of infinite rapidity, the angles should be completely uncorrelated. In the regime of large, but realizable rapidity differences the resummation of large logarithms calls for a description within the BFKL theory. Unfortunately, the leading logarithmic approximation [7,8] overestimates this decorrelation by far. Improvements have been obtained by taking into account some corrections of higher order like the running of the coupling [9,10]. In particular, the effect of energy-momentum conservation, which is beyond BFKL approximation, was shown to have an important impact for large rapidity separation of jets [9]. In our present study this can affect the reliability of the predictions at the borders of the phase space, as discussed in Sec. 3.3.2. Some earlier calculations with the next-to-leading (NLL) BFKL Green's function have been published in Ref. [11,12].
In this paper we present the full NLL BFKL calculation where also the NLL result for the Mueller Navelet vertices [13,14] will be taken into account. In Sec. 2 we recall the LL BFKL calculation deriving also the key formulas which are then used in Sec. 3 where the NLL calculation is presented. In Sec. 4 we present results and give a summary in Sec. 5. Technical details of the numerical implementation are given in an appendix.

Kinematics
The kinematic setup is schematically shown in Fig. 2. The two hadrons collide at a center of mass energy s producing two very forward jets, the transverse momenta of the jets are labeled by Euclidean two dimensional vectors k J,1 and k J,2 , while their azimuthal angles are noted as φ J,1 and φ J,2 . We will denote the rapidities of the jets by y J,1 and y J,2 which are related to the longitudinal momentum fractions of the jets via x J = |k J | √ s e y J . At any real experiment transverse momenta as well as rapidities are measured within certain intervals. A proper theoretical calculation should take this into account and integrate |k J,i | and y J,i over the according interval. However, since at the LHC the binning in rapidity and in transverse momentum will be quite narrow [15], we consider the case of fixed rapidities and transverse momenta.

LL BFKL calculation
Due to the large longitudinal momentum fractions x J,1 and x J,2 of the forward jets, collinear factorization holds and the differential cross section can be written as dσ d|k J,1 | d|k J,2 | dy J,1 dy J,2 = a,bˆ1 0 dσ ab d|k J,1 | d|k J,2 | dy J,1 dy J,2 , (2.1) where f a,b are the standard parton distribution functions (PDFs) of a parton a (b) in the according proton. They depend furthermore on the renormalization scale µ R and the factorization scale µ F .
The partonic cross section at lowest order in the collinear factorization approach would just be described by simple two-to-two scattering processes as they are discussed in standard text books. However, the necessary resummation of logarithmically enhanced contributions calls for a description of the partonic cross section in k T -factorization: dσ ab d|k J,1 | d|k J,2 | dy J,1 dy J, 2 =ˆdφ J,1 dφ J,2ˆd 2 k 1 d 2 k 2 V a (−k 1 , x 1 )G(k 1 , k 2 ,ŝ)V b (k 2 , x 2 ), (2.2) where G is the BFKL Green's function depending onŝ = x 1 x 2 s, and the jet vertex V at lowest order reads [13,14]: In the definition of h (0) a , C A = N c = 3 is to be used for initial gluon and C F = (N 2 c − 1)/(2N c ) = 4/3 for initial quark. Following the notation of Ref. [13,14], the dependence of V on the jet variables is implicit.

NLL calculation
The master formulae of the LL calculation (2.18, 2.19) will also be used for the NLL calculation. Even though the vertices do not simplify as drastically as in the LL case, we gain the possibility to calculate for a limited number of m the coefficients C m,ν as universal grids in ν. In transverse momentum space one would need a two dimensional grid. Moreover, at NLL there are some contributions with an additional transverse momentum integration, such that some contributions would be analytic functions in e.g. k 1 while other would be proportional to distributions like δ (2) (k 1 − k J,1 ).

Strong coupling, renormalization scheme and PDFs at NLL
Based on the MS renormalization scheme, we use the MSTW 2008 PDFs [16] and the two-loop strong coupling in the following form: In the following α s orᾱ s without argument is to be understood as α s (µ 2 R ) orᾱ s (µ 2 R ) respectively. Since in the MSTW 2008 PDFs µ R and µ F are set to be equal, for a consistent calculation we are forced to perform this identification throughout the whole calculation as well.

Jet vertices at NLL
To calculate the coefficients C m,ν (2.19) at next to leading order level, we take for V a (k, x) instead of just the leading order result The matrix elements needed to calculate the Mueller Navelet jet vertex at next to leading order -namely the partonic 2 → 3 process at tree level and the partonic 2 → 2 process at one loop level -are known for a long time. The separation of collinear singularities (to be absorbed by renormalized PDFs) from the BFKL large logarithms in s was performed by Bartels, Vacca and one of us [13,14] in terms of a generic and infrared-safe jet algorithm.
In this paper, we shall apply such procedure to a concrete jet algorithm, namely the cone algorithm, as will be explained in Sec. 3.2.1. We will build on the results obtained in Ref. [13,14] using their notation as well, but we correct an inconsistency in the treatment of the collinear cutoff parameter Λ which later is identified with the factorization scale µ F . In the 'real' C F term the authors rescale the transverse momentum which is integrated over but do not adapt the cutoff parameter Λ. The correction of this point does not change the singular terms, and all the discussion of the arrangement of divergences and subtractions remains unchanged. However the finite part of the subtraction changes such that beside the cutoff functions also the 'virtual' part of the vertex changes, e.g. the term proportional to ln(1−z) The final expressions for the NLL correction to the vertices read: Here N f denotes the number of active quark flavors, b 0 = (11N c − 2N f )/(12π), and N = α s / √ 2. A priori, the factorization scale µ F = Λ and the renormalization scale µ R = µ are independent of each other even though in the end we will set them equal.

Jet definition
For a concrete calculation of Mueller Navelet jet production one also has to choose a concrete jet algorithm obeying the property of infra-red safety, as required by the general procedure of Refs. [13,14]. Two of the most common ones are the cone algorithm as it has been adapted for NLL calculation in Ref. [17], and the k T algorithm where In our study we will use the cone algorithm with a cone size of R cone = 0.5 as it probably will be used in a CMS analysis at the LHC [15].

LL subtraction and s 0
The requirement of a BFKL calculation that the two scattering objects have a similar hard scale is reflected by the fact that in this standard situation of BFKL physics the energy scale s 0 can be written as a product of two energy scales each assigned to one of these scattering objects.
In Ref. [13,14] the energy scale s 0,i (assigned to the Mueller Navelet jet) was chosen as (|k J | + |k J − k|) 2 . While k is integrated over, it is preferable to let s 0 depend only on external scales. Alsoŝ = x 1 x 2 s is in fact not an external scale since the longitudinal momentum fractions x 1 and x 2 are integrated over as well. Therefore, we want to change to a new s ′ 0 : 10) where we introduced the relative rapidity Y = y J,1 − y J,2 .
The energy scale s 0 is a free parameter in the calculation. However, like for the renormalization scale at NLL level a change of it does not go without consequences. In fact, a change of s 0 → s ′ 0 in the Green's function has to be accompanied by an according correction term to the impact factors [18,19]: with K LL being the LL BFKL kernel. Due to the Dirac delta distribution δ(1 − x J,i /x i ) in the jet algorithm inside Φ the ratio of longitudinal momentum fractions in s ′ 0,i reduces to 1 and hence the logarithm in Eq. (3.13) vanishes for k ′ i = k i such that only the real part of the kernel contributes.
To study the role of s 0 , we will investigate the effect when changing it. A subsequent change of s 0,i by just a factor λ can be easily performed at the very end because of the use of BFKL eigenfunctions: (3.14) The LL subtraction, i.e. the terms multiplied by Θ(|k − k ′ | − z(|k − k ′ | + |k ′ |)) in Eqs. (3.5, 3.4), cancels some part in the limit of the additional emission having a big rapidity distance to the jet. In fact, numerically this cancellation works very poorly due to an azimuthal averaging which has been performed for the LL subtraction. A significant improvement can be obtained by omitting this averaging and introducing new LL subtraction terms As a consequence s 0,i changes from It is also possible to use which are slightly inferior concerning the numerical performance but give a s 0 change from s 0,i = (|k J,i | + |k J,i − k i |) 2 to s 0,i = k 2 J,i which already is close to the final s ′ 0 making a correction term (Eq. (3.13)) needless since the ratio of longitudinal momentum fractions in s ′ 0,i effectively reduces to 1 as described above.
We have checked that all three possible subtraction terms after combining them with the according correction term (3.13) lead to the same result. For reasons of numerical performance we have chosen Eqs. (3.15) for the final calculation. As we nevertheless aim for the final s ′ 0 defined in Eq. (3.12) we still have to use the correction term introduced in Eq. (3.13) additionally.

BFKL Green's function at NLL
Last but not least, we also have to take the BFKL Green's function at NLL level. The key to the Green's function is the BFKL kernel at NLL [20,21]. While at LL the BFKL equation is conformally invariant, at NLL it is not such that in fact the LL eigenfunctions E n,ν (2.10) are strictly speaking not eigenfunctions of the NLL kernel. Nevertheless, the action of the NLL BFKL kernel on the eigenfunctions has been calculated in Ref. [22]. The status of the E n,ν being eigenfunctions formally can be saved if one accepts the eigenvalue to become an operator containing a derivative with respect to ν [11,23,24]. In combination with the impact factors the derivate acts on the impact factors and effectively leads to a contribution to the eigenvalue which depends on the impact factors [11,[23][24][25]: where with the constant S = (4 − π 2 + 5β 0 /N c )/12. ζ(n) = ∞ k=1 k −n is the Riemann zeta function while the function φ reads φ(n, γ) = ∞ k=0 (−1) k+1 k + γ + n 2 ψ ′ (k + n + 1) − ψ ′ (k + 1) (3.20) At NLL accuracy, only the leading order vertex coefficients (2.20) enter in the derivative term of (3.17): (3.21)
In general the additional ν-derivative term makes it necessary to recalculate the coefficients d 1,k (defined in Ref. [30]) but in our case the LL vertex does not contain any poles in γ nor in 1 − γ leaving the coefficients d 1,k unchanged 2 . By introducing an ω dependence in the eigenvalue the pole in (2.14) is no longer a simple one such that the residue in fact reads

Approximate energy-momentum conservation in BFKL
We would like to finish this section with the following important observation. Energymomentum conservation is not fulfilled in any truncated BFKL treatment (i.e. LL BFKL, or NLL BFKL, etc.) while it is preserved in any truncated fixed order treatmentà la DGLAP (LO, NLO, etc.). Our approach uses a semi-analytical resummed solution of the BFKL equation at NLL and does not allow, in a direct way, for the implementation of a procedure based on the iteration of the BFKL kernel (in the spirit of Ref. [9]) in which energy-momentum conservation could be imposed step by step. Exact energy-momentum 2 The same is true for the NLL vertex (except for the s0-correction term (3.13)) as can be seen by using a closed contour in γ-plane around 0 or around 1, and numerically integrating integer powers of γ or 1 − γ times the vertex. Based on Cauchy's formula -used here in reverse manner -one can then obtain a numerical evaluation of the residue of arbitrary order, and show that they actually vanish.
conservation is beyond the scope of our pure BFKL treatment. Therefore, our results are expected to be valid only in a limited range of relative rapidity between the two Mueller Navelet jets, away from the kinematical bounds. Nevertheless, this violation within BFKL approach is less dramatic when going further in the truncation (note that BFKL and DGLAP are expected to converge to same result when going higher and higher in the order of perturbation). We thus expect that in the region far from the kinematical limit, such energy-momentum conservation effects would not introduce significant corrections to our NLL BFKL results.

Results
We now present results for the LHC at the design center of mass energy √ s = 14 TeV.
Motivated by a recent CMS study [15] we restrict the rapidities of the Mueller Navelet jets to the region 3 < |y J | < 5. We shall show the differential cross section with respect to the relative rapidity variable Y = y J,1 − y J,2 which therefore takes values between 6 and 10. Note that, since the maximum possible rapidity of a jet with 50 GeV of transverse energy (see below) is y max = 5.6, values of Y ≃ 10 are quite close to the kinematical boundary and the corresponding predictions may suffer some uncertainties because of the fact that momentum conservation is not exactly fulfilled in this BFKL approach. We consider Mueller Navelet jets with |k J | = 35 GeV, and |k J | = 50 GeV respectively. Due to our method of calculation and the factorization between the two Mueller Navelet jets we can can combine the building blocks to two symmetric scenarios (|k J,1 | = |k J,2 | = 35 GeV or |k J,1 | = |k J,2 | = 50 GeV) and one asymmetric scenario of |k J,1 | = 35 GeV, |k J,2 | = 50 GeV (plus the 'mirrored' process |k J,1 | = 50 GeV, |k J,2 | = 35 GeV) even though in doing so one mixes different choices for µ R . But since α s only varies by ∼ 4% between 35 GeV and 50 GeV we give the according result as well.
In Ref. [35] it is argued based on Ref. [36] that imposing |k J,1 | > E, |k J,2 | > E + D for D → 0 large logarithms of non-BFKL origin make a fixed order calculation unstable. Even though in a pure BFKL framework these logarithms do not show up at all, they -and their resummation -might be of significant impact on the result of a BFKL calculation. Therefore, D = 0 is preferred also in the BFKL framework to be safe from these unknown contributions. However, we start the presentation of our result with these symmetric scenarios but in order to be conservative and to avoid the region where initial state radiations might require peculiar treatment involving resummationsà la Sudakov, which are beyond the scope of our work (and not implemented in fixed NLO calculation neither), we will close in Sec. 4.3 with our results in the asymmetric case for which such resummations effects are clearly not required.
In all cases we choose the number of active flavors to be five (N f = 5) with Λ QCD = 221.2 MeV such that α s (M 2 Z ) = 0.1176. The Monte Carlo integration [37] itself is error-prone. This error in Monte Carlo integration can be reduced to less than 1% with a large number of sample points, so as to be practically negligible in comparison with other sources of uncertainties. In practice, due to hardware/ time limitations we will display results for a Monte Carlo integration setup (for details see Sec. A.2) which aims for an accuracy of the order of 1%.
However, as we show in what follows, there are more serious uncertainties due to the renormalization scale µ R which we choose as µ R = |k J,1 | · |k J,2 |. To study the dependence on it we vary µ R by factors 2 and 1/2 respectively. The same we do for the energy scale √ s 0 . We investigate the effect of the uncertainty of PDFs for asymmetric errors as defined in Eqs. (51,52) of [16] with the eigenvector set ensuring all data sets being described within their 90% confidence level limits.
Below we present our results. The logic of their presentation is the following. For each kinematical situation, we discuss the physical observables and their uncertainties: crosssection encoded in C 0 , azimuthal decorrelation encoded in C 1 /C 0 , C 2 /C 1 and C 2 /C 1 . We then give further details on the additional underlying quantities C 1 , C 2 . We use the same color coding for all plots, namely blue shows the pure LL result, brown the pure NLL result, green the combination of LL vertices with the collinear improved NLL Green's function, red the full NLL vertices with the collinear improved NLL Green's function. Whenever we show curves for scales µ R = µ F or s 0 changed by factors 2 or 1/2, the thick curve corresponds to the scale changed by factor 2.
Thus the first thing to look at is the differential cross section as defined in Eq. (2.8). The result of our calculation is shown in Fig. 3 (the according tabled values are shown in Tab. 1 in the Appendix). The decrease of the cross section at large values of Y 7 is mostly an effect of the upper kinematical cut on the rapidity of the single jets. This is true also for the coefficients C 1 and C 2 in Figs. 16 and 19. The purely numerical error due to the Monte Carlo integration of the NLL vertices is mainly below 2% (see Fig. 5), and only for very large Y of the order of 2 − 5%. We varied the renormalization and factorization scale by factors 2 and 1/2 to investigate the µ R dependence. A full scan over this interval is not possible due to the CPU time consumption of the evaluation for a single choice of µ R . The results are displayed in Fig. 4. As one would expect, the full NLL result depends less on µ R than the LL result or the combination of LL vertices with resummed NLL Green's function which was so far state-of-the-art [11,12]. However, the results obtained for the three choices of µ R studied here seem to suggest that, as expected, after inclusion of the NLL vertices the µ R dependence is no longer monotone and flattens out, resulting in higher stability of NLL results. Another important scale is the energy s 0 introduced by the Mellin transformation from energy to ω space which is necessary to formulate the BFKL equation. Like µ R it is an artificial scale which in an all order calculation would not affect the result. Indeed, the dependence is reduced when the NLL corrections to the vertices are taken into account (see Fig. 4).
The dependence on PDF uncertainties is shown in Fig. 5. This significant sensitivity is almost identical for pure LL, pure NLL, combined LL vertices with collinear improved NLL Green's function and for full NLL vertices combined with the collinear improved NLL Green's function.
The azimuthal decorrelation, described by coefficients defined in Eq. (2.9), has often be predicted to be a striking feature of BFKL physics. Our results displayed in Fig. 6 for cos ϕ and in Fig. 10 for cos 2ϕ explicitly show that our inclusion of NLL vertices leads to an enormous correlation in the azimuthal angle (for completeness our results for C 1 and C 2 coefficients alone are displayed respectively in Fig. 16 and Fig. 19).
In particular, cos ϕ shown in Fig. 6 is rather close to the typical values predicted by LO-DGLAP Monte Carlos Pythia [38] and Herwig [39], used for CMS studies [15]. Note that Herwig has the tendency to predict more decorrelation, presumably because since it implements more radiations than Pythia, it has the phenomenological effect to involve some kind of NLO-DGLAP corrections, which enhance the decorrelation. The various sources of uncertainty of our results are shown for cos ϕ in Fig. 7 (variation of µ R = µ F ), Fig. 8 (variation of s 0 ), and for cos 2ϕ in Fig. 11, and Fig. 12 accordingly.
Not only is the Y dependence much flatter for the NLL vertices. But the mean value for cos ϕ itself is very close to 1, especially for the NLL calculation including the collinear improved Green's function. Actually, in the collinear improved approach and for low values of µ R , cos ϕ = C 1 /C 0 can exceed unity. This feature has to be ascribed to the fact that the collinear improvement is justified and applied only to the n = 0 conformal spin of the kernel, with the effect of lowering C 0 in a large region of the phase space, while keeping C 1 unchanged. For this reason, the predictions of angular dependent quantities are more trustable in the pure NLL approach (without collinear resummation).
The µ R dependence of C 1 /C 0 (see Fig. 7) is puzzling, but has to be considered as a consequence of the C 0 and C 1 dependencies in Figs. 4 and 17 respectively. A similar behavior can be observed for the s 0 dependence (see Fig. 8 for C 1 /C 0 and Fig. 4 and Fig. 17 for C 0 and C 1 respectively). While the µ R and s 0 dependences for C 0 and C 1 are significantly reduced by the inclusion of the NLL vertices, in the ratio this is only true if compared to the pure LL calculation. The combination of LL vertices and NLL collinear improved NLL Green's function is less sensitive on changes of s 0 or µ R . The reason for this surprising 'weakness' of the NLL result is that the changes of the LL vertices when changing s 0 or µ R are not very sensitive on n such that ratios of LL vertices are very stable against scale changes. In contrast, the NLL correction -especially the LL subtraction -is very large (and negative) for the n = 0 component while of minor significance for n > 0 such that effects of a scale change for NLL vertices do not vanish by considering ratios C n /C 0 . The special role of n = 0 becomes apparent if one compares the situation to C n =0 /C m =0 (see Figs. 15,14) where the expected advantage of the full NLL calculation is clearly visible.
For both scales, the curves exceeding 1 belong to smaller scales (in the calculation with the collinear improved Green's function) which seem to be very disfavored in full NLL BFKL calculations as already discussed in [23,40,41]. The dependence on the PDFs completely drops out for LL vertices, and also for NLL vertices it is negligible as can be seen in Fig. 9.
A similar rather large dependency on µ R = µ F and s 0 is obtained for C 2 /C 0 , although it does not lead to values of cos 2ϕ close to 1, as can be seen from Figs. 11, 12, based on detailed studies of coefficients C 0 and C 2 displayed respectively in Figs. 4 and 20.
In Ref. [11] it has been proposed to also study other observables C m /C n with m = 0 = n. This is motivated by the fact that the main source of uncertainty of the Green's function is associated with the n = 0 component. This observation is not altered by the inclusion of NLL vertices, as we display in Figs. 13-15. Moreover, Fig. 14 shows that the effect of changing µ R = µ F leads to modifications of similar size for pure NLL and combined LL vertices with NLL Green's function predictions, while the changing of s 0 (see Fig. 15) leads to a reduced dependency with respect to s 0 of the pure NLL prediction. The PDF dependence for all ratios C m /C n cancels in the same manner such that there is no use plotting it more than once (see Fig. 9 for C 1 /C 0 ). A priori, one would expect the numeric uncertainties for the C n>0 calculations including the NLL vertices to be larger because the coefficients C n,ν contain an azimuthal integration which in the case of n = 0 becomes trivial while for n > 0 has to be carried out. This is indeed true for the NLL corrections alone, but since these corrections are more significant for n = 0 than for n > 0, in the sum together with the error-free LL vertices the opposite turns out to be true (the Monte Carlo errors are smaller than 1% for C 1,2 as shown in Figs. 18, and 21 respectively).
We would like to draw the reader's attention to the analogous effects for the kernel. For both the BFKL kernel and the Mueller Navelet vertices the NLL corrections to the n = 0 component around ν = 0 are very large and negative while the relative corrections for n > 0 are positive, much smaller than for n = 0 and slowly increasing with n.
Finally, for other kinematical configurations to be discussed below, the same kind of PDF uncertainties appear, and Monte Carlo errors are of similar order. Because of that, we will not display the corresponding curves.   Figure 4: Relative effect of changing µ R = µ F by factors 2 (thick) and 1/2 (thin) respectively (left), and √ s 0 (right) by factors 2 (thick) and 1/2 (thin) resp. on the differential cross section in dependence on Y for |k J,1 | = |k J,2 | = 35 GeV. The tabled values are shown in Tabs. 2 and 3.         0.8  0.8   0.8 Figure 14: Effect of changing µ R = µ F by factors 2 and 1/2 respectively on cos 2ϕ / cos ϕ in dependence on Y for |k J,1 | = |k J,2 | = 35 GeV. The tabled values are shown in Tab. 13. 0.8        Going to larger jet scales, we meet more or less the same advantages and problems as for 35 GeV. Again we start with the differential cross section (2.8). The result is shown in Fig. 22 (the according tabled values are shown in Tab. 23 in the Appendix). The dependences with respect to µ R , and s 0 are displayed in Fig. 23.
The azimuthal decorrelation is displayed in Fig. 24 for cos ϕ and in Fig. 27 for cos 2ϕ , again explicitly showing that inclusion of NLL vertices leads to an enormous correlation in the azimuthal angle (for completeness our results for C 1 and C 2 coefficients alone are displayed respectively in Fig. 33 and Fig. 35). Here, the angular correlation even has the tendency to increase with growing rapidity Y . This might be interpreted as the effect of stronger limited phase space for additional emissions at large energies and large transverse momenta of the produced jets (Note, that the cross section is a factor ∼ 10 smaller at Y = 6 compared to the previous configuration, and a factor ∼ 100 smaller at Y = 10).
The various sources of uncertainty of our results are shown for cos ϕ in The scale dependences of C 1 /C 0 (see Figs. 25, and 26) as well as of C 0 (see Fig. 23) and C 1 (see Fig. 34) alone reveal the same basic features as before, namely a non-monotone scale dependence of the NLL corrections and (more serious) unphysical results for cos ϕ in case of the resummed NLL prediction for small s 0 and/ or µ R = µ F scales.
A similar rather large dependency on µ R = µ F and s 0 is obtained for C 2 /C 0 , as can be seen from Figs. 28, 29, based on detailed studies of coefficients C 0 and C 2 displayed respectively in Fig. 23 and Fig. 36.
The problematic behavior for smaller scales of s 0 and/ or µ R is more dramatic for |k J,1 | = |k J,2 | = 50 GeV (see e.g Figs. 25,26,28,29). Especially the µ R dependence (see Fig. 25 seems to indicate that already the a priori natural scale µ R = |k J | is too small.         0.8  0.8   0.8 Figure 31: Effect of changing µ R = µ F by factors 2 and 1/2 respectively on cos 2ϕ / cos ϕ in dependence on Y for |k J,1 | = |k J,2 | = 50 GeV. The tabled values are shown in Tab. 33. 0.8      We end up with the consideration of the asymmetric case, which we investigate in order to provide a comparison with NLO-DGLAP predictions [42] obtained through the NLO-DGLAP partonic generator Dijet [43]. These prediction are very sensitive to the precise compensation between the real and the virtual contribution, and a symmetric cut leads to some kind of Sudakov resummation effects which are not completely under control at the moment [44], even leading to a negative cross-section for |k J,1 | = |k J,2 | = 35 GeV. These prediction are much more stable in the asymmetric configuration. Our own predictions for the cross-section, for cos ϕ , cos 2ϕ , and cos 2ϕ / cos ϕ are given respectively in Figs. 37, 39, 42 and 45. Due to the factorization, the sensitivity of our prediction with respect to s 0 , µ R is similar to the two previous symmetrical configurations, as shown in Figs One sees from Fig. 37 that our pure NLL prediction, as well as our resummed NLL prediction, are a bit below the NLO-DGLAP prediction, while the LL prediction is much higher than the NLO-DGLAP prediction. The combined LL vertices plus resummed NLL Green's function is rather close to the NLO-DGLAP prediction. One may expect that including higher order corrections in both DGLAP and BFKL approaches would make them converging. We note however that comparing both kinds of treatment should be done with some cautious. Indeed, the NLO-DGLAP involves scales which are smaller than the scale which we consider: we take µ R = |k J,1 | · |k J,2 | which is similar to (|k J,1 | + |k J,2 |)/2, while the NLO-DGLAP calculation uses the scale (|k J,1 | + |k J,2 |)/4 . Changing this scale from (|k J,1 | + |k J,2 |)/4 to (|k J,1 | + |k J,2 |)/8 leads to a variation of the order of 5% in the NLO-DGLAP prediction. Our treatment, especially when considering the azimuthal decorrelation, favors higher scales, like |k J,1 | · |k J,2 | ∼ (|k J,1 | + |k J,2 |)/2 or even 2 |k J,1 | · |k J,2 | ∼ |k J,1 | + |k J,2 | .
The azimuthal decorrelation, which is expected to be the best signal, is predicted to be similar in magnitude and shape both from our pure NLL prediction and our resummed NLL prediction and from the NLO-DGLAP approach, as can be seen from Figs. 39 and 42. Note however that the uncertainties of our predictions are rather high. Anyway, the general trend is clear: the azimuthal decorrelation is much lower than expected from a LL BFKL treatment or from a mixed treatment with LL vertices combined with NLL Green's function. It is also rather flat with Y . The only observable which still remain different when comparing pure NLL approaches (the resummed NLL approach makes no difference here since it only affects C 0 ) with NLO-DGLAP is the ratio cos 2ϕ / cos ϕ for which the NLO-DGLAP is still significantly higher than the NLL prediction, as can be seen from           0.8 0.8

Conclusions
We have implemented at full NLL order the Mueller Navelet jets cross-section as well as their relative azimuthal angle dependency. In contrast to the general belief, the effect of NLL corrections to the vertex function is very important, of the same order as the one obtained when passing from LL to NLL Green's function. The importance of NLL corrections to the impact factor observed in the present paper is analogous to recent results obtained at NLL in diffractive double ρ-electroproduction [23,40]. Interestingly, the full NLL calculations for cos ϕ and cos 2ϕ are quite close to a calculation [42] using Dijet [43] which is based on DGLAP dynamics and to a dedicated study [15] using Pythia [38] and Herwig [39]. The uncertainty due to changes in µ R (and s 0 ) is drastically reduced for all C n when one takes into account the NLL Mueller Navelet vertices. The uncertainty due to PDFs are also moderate. As a consequence, our results for the cross-section are very stable.
However, for azimuthal decorrelation the dependence on µ R (and s 0 ) is still sizeable. In the case of the NLL Green's function with collinear improvement one observes that cos ϕ can exceed 1 for certain choices of the parameters, in particular for low values of µ R = µ F , taken to be smaller than the "natural" value |k J,1 | · |k J,2 | . One might also think of a collinear improvement of the vertices [33] but the Mueller Navelet vertex for fixed |k J | does not have poles in γ nor 1 − γ, so there is no room for such a treatment. The resummation of soft initial radiation might be of relevance for the azimuthal correlation as well. This is left for further investigations, and in this work we rather consider the full NLL calculation without additional collinear resummation to be our solid prediction, while the 'collinear improvement' as it stands is not appropriate to study azimuthal dependences.
At present, there is little experience with the effect of NLL impact factors. To the best of our knowledge, up to now, the only full NLL BFKL calculation existing in the literature is the vector meson production in virtual photon collisions [23,40,41], which is very sensitive to NLL corrections to the impact factor and for which very large values for s 0 and µ R are preferred. In [41] it has been shown that a collinear improved treatment combined with the application of the principle of minimal sensitivity [45,46] reduces this large values to more "natural" values. Still, µ R larger than the "natural" values are favored [41]. In the present case, with the scales µ R and s 0 set by the jet scale, we get azimuthal correlations which are rather similar to DGLAP dynamics predictions (although, as we already mentioned, the DGLAP prediction are based on smaller scales). To conclude, contrarily to the expectation, it thus seems that the azimuthal decorrelation is almost not enhanced by an increasing rapidity. This suggests that the study of Mueller Navelet jets is probably not the best place to exhibit differences between BFKL and DGLAP dynamics. own predictions. This work is supported in part by the Polish Grant N202 249235, the French-Polish scientific agreement Polonium, by the grant ANR-06-JCJC-0084 and by the ECO-NET program, contract 12584QK, and by a PRIN grant (MIUR, Italy).

A.1 Programs used
We implemented all numerical calculations in Mathematica. To this purpose we used the according interfaces for the MSTW 2008 PDFs [16] and for version 1.5 of Cuba [37] which we used for numerical integration.

A.2 Choice of Parameters
Cuba provides different integration routines which we also used to cross-check the results of the Monte Carlo integration. However, for the final results we used the Vegas routine of Cuba with an aimed precision of 10 −2 and a maximal number of 500 000 points per integration. To use a Monte Carlo integrator, all integration intervals have to be mapped on finite intervals. For the transverse momentum integrations we used the mapping |k| = |k J | tan(ξπ/2).
The cancellations which analytically have been shown in Refs. [13,14] numerically can corrupt the integration due to the limited precision of a computer. In all these cases, were the integration interval have been mapped to the compact interval [0, 1], we used a cut off of 10 −5 where the cut off dependence becomes negligible.

A.2.1 The ν-grid
Due to the complicate matrix element, the PDF evaluation, and the implementation in Mathematica instead of a dedicated stand-alone code the Monte Carlo integration is very time consuming. Therefor, the choice of the ν-values at which the coefficients C n,ν (2.19) are evaluated is crucial.
We are guided by the shape of the BFKL Green's function which is peaked around ν = 0 and then monotonically falls. The smaller Y the slower the decrease. Even though the minimal Y in this study is 6, we want our coefficients to be prepared also for smaller Y 's. We choose a maximal ν max such that an integration up to ν max of just the NLL BFKL Green's function at Y = 4 for n = 0 reproduces 96% of the integration over the full ν-range. For the case of Y = 6 it reproduces 99.97% of the full integral.
The coefficients C n,ν are oscillating like exp(iν ln k 2 J,i ) but more important is the product of the two which has an oscillating part with a frequency ν oscillation = π/ ln |k J,1 | |k J,2 | . This frequency is zero for |k J,1 | = |k J,2 | but for the example of |k J,1 | = 35 GeV and |k J,2 | = 50 GeV we chose a step width for ν of ν oscillation /4.
For large Y it is really the small region close to ν = 0 which matters. Therefore we sample this region in more detail according to the shape of the NLL BFKL Green's function For the final integration over ν the product C n,ν (|k J,1 |, x J,1 )C * n,ν (|k J,2 |, x J,2 ) is interpolated by cubic splines.

A.3 Grouping the integrand
In this section we describe how the NLL contribution to the coefficients C n,ν , as defined in Eq. (2.19), is arranged.
The jet defining function S It is useful to replace the integration variable k by k → k J − k in the integrands V q [5] and V q [7] (V g [23] and V g [25]). Moreover we split up V q [3] (V g [21]) and create V q [18] (V g [26]) where in the new elementary blocks the integrand k ′ is replaced by k J − k ′ . Then we make the following replacements The elementary blocks are now grouped to 14 minimal basic blocks B[i] which also contain the integrations from Eq. (2.19) where we made use of the short hand notationˇ≡´dφ J d 2 k dx f (x)E n,ν (k) cos(mφ J ). We would like to point out that the inclusion of V q [18] (V g [26]) in B [4] (B [14]) is essential. Even though it is correctly stated after Eq. (88) in Ref. [13] (and repeated in Ref. [14] after Eq. (53)), that in the composite jet configuration the domain of integration shrinks like z 2 for z → 0, the conclusion that this prevents a divergence is wrong. In fact, in the limit z → 0 the integrand scales like z −3 and only the sum of V q [18] (V g [26]) and V q [5] (V g [23]) cancels properly against V q [7] (V g [25]) in the dangerous region.
Note that in case of fix x J for B [1] and B [7] no numerical integration is needed, while for B [2] and B [8] only one integration over z has to be done. All four are proportional to the LL Mueller Navelet vertex regarding the transverse momentum dependences. The Dirac-δ in transverse momenta we always use for the k integration. Only for contributions with δ (2) (k ′ − k J ) it is used for the k ′ integration. The x-integration is trivially performed by evaluating the according Dirac δ-distribution.

B. Tabled values of diagrams
To allow for later accurate comparisons, we give the values for all plots in this work. We mark the pure LL calculation by 'LL', and the pure NLL one by 'NLL'. The combination of LL vertices with the NLL collinear improved Green's function is denoted as 'LL+', and the combination of NLL vertices with the NLL collinear improved Green's function as 'NLL+'. Whenever in a figure the effect of the variation of one parameter is presented, in the according table the first column shows the central value, while the second and third show the change of this central value due to the varied parameter. For brevity we suppress the energy unit GeV in the headings of the tables.                                                             Fig. 51.