Resummation of non-global logarithms at finite $N_c$

In the context of inter-jet energy flow, we present the first quantitative result of the resummation of non-global logarithms at finite N_c. This is achieved by refining Weigert's approach in which the problem is reduced to the simulation of associated Langevin dynamics in the space of Wilson lines. We find that, in e+e- annihilation, the exact result is rather close to the result previously obtained in the large-N_c mean field approximation. However, we observe enormous event-by-event fluctuations in the Langevin process which may have significant consequences in hadron collisions.


Introduction
In certain search channels of the Higgs boson and new particles at the Large Hadron Collider (LHC), it is often desirable to be able to control QCD radiation from tagged jets in order to suppress large backgrounds. A prime example is the Higgs boson production in association with di-jets. The two competing production mechanisms, the gluon fusion and the vector boson fusion processes, have different patterns of soft gluon radiation due to the difference in their color structure. Discriminating these processes quantitatively using some measure of radiation is therefore a useful strategy to determine the Higgs couplings [1,2].
A related class of observables which are particularly sensitive to soft radiation is the cross section with a veto on unwanted jets in the full or partial region of the phase space. This generally requires the resummation of logarithms in p veto T , the threshold transverse momentum of vetoed jets. Steady progress in this direction has been made for global observables which involve all the particles and jets in the final state including those close to the beam axis. The state-of-the-art is that one can resum the leading logarithms (LL) (α s ln 2 p veto T ) n , the nextto-leading logarithms (NLL) (α s ln p veto T ) n and even the next-to-next-to-leading logarithms (NNLL) (α 2 s ln p veto T ) n [3,4,5].
However, in contrast to such progress, there exists a severe limitation in our ability to resum non-global logarithms [6,7] which arise when measurements are restricted to a part of the phase space excluding the beam and jet regions. In this case, double-logarithms (α s ln 2 p veto T ) n are absent due to the lack of the collinear singularity. The leading contribution is then comprised of single-logarithms (α s ln p veto T ) n which originate from the soft singularity. The problem is that these logarithms do not exponentiate, and because of this difficulty their resummation has been hitherto done only in the large-N c limit [6,7,8,9]. In other words, even the leading logarithms cannot be fully satisfactorily resummed. This could be a potentially serious drawback in actual experiments considering the fact that, strictly speaking, any vetoed cross section at hadron colliders is inevitably non-global due to the finite acceptance of detectors.
In fact, there is a single work by Weigert [10] which did discuss the resummation of nonglobal logarithms at finite N c in a simpler setup of e + e − annihilation where complicacies from the initial state radiation do not arise. His approach is based on an analogy with another, seemingly unrelated resummation in QCD, namely, that of the small-x (or 'BFKL') logarithms in Regge scattering. In this context, a very similar issue arises as to how one can generalize the equation which resums small-x logarithm in the large-N c limit to one at finite N c . It turns out that these equations bear a striking resemblance to the equation which resums non-global logarithms in the large-N c limit. Since technologies to solve the former problem are welldeveloped, they may be suitably adapted to address the latter problem as well. Somewhat surprisingly, however, Weigert's approach has not been pursued for a decade. Part of the reason of this may perhaps be that, as we shall point out, there is actually a flaw in his formulation which deters a straightforward numerical implementation. In this paper we overcome this difficulty and present the first quantitative results of the resummation of non-global logarithms at finite N c .
In Section 2, we quickly review the nonlinear evolution equations which resum the small-x logarithms in high energy QCD. The subject may seem utterly unfamiliar to the readers whose primary interest is jet physics. However, the similarity (or even equivalence) to the resummation of non-global logarithms will soon become apparent in Section 3 where we introduce the relevant evolution equations. Using this similarity, we discuss how to solve the equation for non-global logs at finite N c in Section 4, and present numerical results in Section 5. Finally, we examine the results and conclude in Section 6.

B-JIMWLK equation
Consider a 'dipole' consisting of a quark and an antiquark located at transverse coordinates x and y, respectively. The S-matrix of the dipole moving in the and scattering off some target in the eikonal approximation is where U x is the Wilson line in the fundamental representation and τ is the rapidity of the dipole. The expectation value · · · is taken in the target wavefunction. In deep inelastic scattering (DIS), τ is essentially the logarithm of the Bjorken-x variable The dipole S-matrix (1) (more precisely, 1 − S ) then measures the total cross section of the subprocess γ * p → qqp → X at the corresponding values of x and the photon virtuality The B-JIMWLK 1 equation [11,12,13] resums the small-x logarithms τ n ∼ (α s ln 1/x) n which arise in the high energy evolution of S . It reads where the dipole kernel is given by The equation (4) is not closed because the right-hand-side contains the double dipole S-matrix SS . This non-linearity reflects the gluon saturation effect in the target. The consequence is that (4) is actually the first equation of an infinite hierarchy of coupled equations. However, it becomes a closed equation-the Batlisky-Kovchegov (BK) equation [11,14]-if one truncates the hierarchy and assumes factorization SS → S S While the BK equation (6) can be solved numerically in a straightforward manner, solving the B-JIMWLK equation (4) had been difficult until the ingenious reformulation of the problem as random walk [15]. To explain this, it is essential to write the equation in the operator form The effective HamiltonianĤ takes the form whereŨ ab = (Ũ † ) ba is the Wilson line in the adjoint representation and K is the gluon emission kernel In the last line of (8), we defined The derivative ∇ acts on the Wilson line in the fundamental representation as and also on the adjoint Wilson line with ( The latter operation has been used in obtaining the σ-term in (8). An important observation is that (7) may be viewed as the Fokker-Planck equation treating Wilson lines U as dynamical variables. As is well-known, there is an associated Langevin equation which describes the random walk of these variables. The latter can be simulated numerically on a lattice, and the solution of the B-JIMWLK equation has thus been obtained in [16]. We shall later discuss this approach in detail.
Before leaving this section, it is useful to show another form of the Hamiltonian derived in This Hamiltonian can be used only when it acts on 'gauge invariant' operators of the form and generates the same equations as those obtained from (8). Eq. (13) features the dipole kernel (5) instead of the gluon emission kernel (9). The following relation between the two kernels is worth noting

Non-global logs at finite N c
We now turn to the resummation of non-global logarithms which is our primary interest. Consider, as in the final state of e + e − annihilation, a pair of jets that is overall color singlet and pointing in the direction of the solid angle Ω α,β = (θ α,β , φ α,β ) measured with respect to the positive z-axis (see Fig. 1). Let C in be the region inside a pair of back-to-back cones with opening angle θ in which include the jets, and let C out be its complementary region. We then ask what is the probability P (Ω α , Ω β ) that the total flow of energy into C out is less than E out . In the perturbative calculation of P in the regime Q ≫ E out ≫ Λ QCD where Q is the hard scale, logarithms of the form (α s ln Q/E out ) n appear which have to be resummed. These logarithms are non-global because the measurement is done only in C out . To leading logarithmic accuracy, one may identify E out with the jet veto scale p veto T mentioned in the introduction. It has been shown by Banfi, Marchesini and Smye (BMS) that P satisfies the following evolution equation [8] where we abbreviated P αβ = P (Ω α , Ω β ) and the parameter τ is defined as The equation (16) resums logarithms τ n ∼ (α s ln E out ) n to all orders. The integral kernel is composed of null vectors n µ = (1, sin θ cos φ, sin θ sin φ, cos θ) = (1, n) proportional to the four-vector of hard partons. The 'step function' ensures that real gluons are emitted only in C in . 2 A trivial rewriting of (16) where Θ out (γ) ≡ 1 − Θ in (γ) illuminates the physical meaning of the right-hand-side. The first term represents the familiar Sudakov suppression, whereas non-global logarithms are resummed by the second term describing the emission of an arbitrary number of gluons into C in which then coherently emit the softest gluons into C out . For the purpose of the present work, it is indispensable to note that the equation (16) has been derived in the large-N c limit. Its finite-N c generalization was discussed by Weigert [10] and the result reads where 2Nc . However, (21) itself is highly formal in that the meaning of the averaging · · · can be specified only indirectly, as a proxy of certain complicated functional integrals [10]. Nevertheless, putting this qualification aside, the striking similarity between Eqs. (16), (21) and Eqs. (4), (6) is unmistakable. In fact, it is possible to establish a rigorous mathematical equivalence between the two problems. As shown in [18,19], a conformal transformation known as the stereographic projection 2 In principle, real gluons can be directly emitted from hard partons into C out provided their energy is less than E out . However, to leading logarithmic accuracy such contributions may be omitted [8].
exactly maps the respective kernels onto each other Aside from the kinematical constraint factor Θ in , the map (23) dispels any structural difference between the two sets of equations. 3 This equivalence probably has a deep geometrical origin which goes beyond the perturbative framework. Indeed, such a correspondence persists even in the strong coupling limit of N = 4 supersymmetric Yang-Mills theory [18]. Fortunately, a powerful machinery to solve the JIMWLK problem (4) is available, and this brings hope that the jet problem (21) can be solved in a similar manner. For this purpose, one seeks the operator form of (21) with a formal identification Here the Wilson lines U α,β represent eikonal (jet) lines starting from the space-time origin and extending to infinity in the direction of Ω α,β . 4 Eq. (21) can then be written as a Fokker-Planck equation in group space where [10] In (26), the derivative ∇ is defined by and similarly to (12) in the adjoint case (t a → T a ).

Equivalent Langevin dynamics
The Fokker-Planck equations (7) and (25) can be solved by making use of an equivalent Langevin formulation. To illustrate the idea, consider the following Fokker-Planck equation for some probability distribution where χ ab = χ ba and σ a ≡ 1 2 ∂ b χ ba . We assume that χ is factorized in the form χ ab = E ac E cb . The equivalent Langevin equation is then where ξ is the Gaussian white noise characterized by the correlator The distribution can then be obtained by averaging over an ensemble {x a (τ )} of random walk trajectories To be more precise, the equation (31) makes sense only in a τ -discretized form. There is a well-known ambiguity in how we discretize the equation, the so-called Itô-Stratonovich dilemma. The appropriate choice corresponding to (30) is the Itô scheme where ε is the time step and the argument of E is evaluated at the previous time τ , respecting causality [20]. In (34), we have rescaled the noise as √ εξ → ξ in order to make explicit the fact that the typical variation ∆x is O( √ ε) in a random walk. With this normalization, the noise correlator in discrete time reads The operator form of the B-JIMWLK equation (7) has precisely the structure (30) with a factorized kernel (8). It can thus be described by the Langevin dynamics (34), with the SU(3) matrices U x,y playing the role of {x a } [15]. One can simulate the random walk in group space on a transverse lattice, and this is how the solution to B-JIMWLK equation has been originally obtained [16].
In [10], Weigert suggested to follow the same strategy in solving (25). A comparison of (26) and (30) implies that σ = 0. The kernel of (26) can be written as a sum of three factorized terms. [Note that (Θ in ) 2 = Θ in .] Exploiting the factorized form of M as seen in the last line of (18), Weigert deduced an analog of the stochastic term E ac ξ c in (31) by introducing three independent noises ξ (I) , (I = 1, 2, 3) which generates a random walk U → Ue it a E ac ξc in group space. However, the fact that M is factorized in four-vector space means that the noises must have the correlator which is negative for the spatial components µ, ν = 1, 2, 3, hence they cannot be simulated in practice.
The resolution of this problem again comes from the correspondence with the B-JIMWLK evolution. Firstly, the identity (23) clearly shows the two-dimensional nature of the dipole kernel, so introducing four-vectors is an excess. We then notice the close similarity between (26) and (13), the latter being equivalent to (8). This suggests that we can rewrite the effective Hamiltonian in a form analogous to (8) also for the jet problem. We thus look for a kernel K satisfying (cf. (15)) The solution is found to be (cf. (9)) which has a factorized structure on the unit sphere (|n| = 1) embedded in three spatial dimensions. Reversing the argument in [17] which led (8) to (13), we arrive at an effective Hamiltonian equivalent to (26) where Moreover, differently from (36), we write the expression in the brackets as a sum of two factorized terms This reduces the number of independent noises from three to two ξ (I) (I = 1, 2). They are characterized by the following correlator (in discretized 'time' τ , cf. (35)) where k, l = 1, 2, 3 are the spatial indices. The problem of negative metric has been circumvented. We can now write down the associated Langevin evolution for U α with Ω α ∈ C in (cf. (34)) where All the matrices on the right-hand-side of (45) are evaluated at τ according to the Itô scheme. Expanding the exponential up to O(ε), we get The second line of (48) can be simplified as follows. In the σ-term, we use the identitiesŨ ab t b = U † t a U and t a Ut a = 1 2 (trU) − 1 2Nc U to obtain To the accuracy of O(ε), the terms quadratic in noise ξξ may be replaced by their expectation values using (44). After these manipulations, (48) takes the form Note that there is no singularity at Ω γ = Ω α ∈ C in . By computing the difference to O(ε) and using (44), (39) and the relation K αα (γ) = −1/(1 − n α · n γ ), one can recover (21) after the identification (24). For the sake of numerical simulations, it is more economical to express the evolution (50) in a left-right symmetric form 5 where In this representation only the terms proportional to noise are kept in the exponential. It is easy to check that (52) and (50) are equivalent to O(ε) under the identification ξξ ≈ ξξ . Eqs. (44) and (52) will serve as the starting point of our numerical simulation.

Numerical simulation
We simulate the random walk (52) of Wilson lines living on the unit sphere by discretizing the coordinates 1 ≥ cos θ ≥ −1 and 2π > φ ≥ 0 with lattice spacings a c and a φ , respectively. The SU(3) matrices U α are defined only at the grid points belonging to C in , whereas the noises ξ (I) α are defined at all grid points. 6 The initial condition at τ = 0 is simply given by for all Ω α ∈ C in , or equivalently, P αβ = 1 for all pairs (Ω α , Ω β ) corresponding to no radiation before evolution. 7 We then update {U α } after each time step ε according to the formula (52) with noises ξ (I) (I = 1, 2) randomly generated from the Gaussian distribution In order to ensure that U's remain unitary during the evolution, we need to evaluate the exponential of matrices e iA L/R accurately (although the equation (52) makes sense only to O(ε)). In practice, we use an approximation e iA = (e iA/2 n ) 2 n ≈ 1 + iA 2 n + · · · + 1 m! iA 2 n m 2 n with m, n large enough. On top of this, we perform the 'reunitarization' of U's using polar decomposition method after every 100 steps of evolution. The effect of this latter operation is actually very small due to our accurate evaluation of e iA L/R . The above procedure is repeated a desired number of times N = τ /ε, and at the end of this random walk trajectory we compute the trace We then average it over many such trajectories and identify the result with P αβ τ . In practice, we choose ε = 10 −5 and average over 500 independent trajectories. Fig. 2 shows the evolution of P τ (Ω α , Ω β ) for cos θ α = 1 and cos θ β = −1, corresponding to back-to-back jets in the beam direction. 8 The opening angle of the cones (see Fig. 1) is fixed to cos θ in = 1 2 . The simulation was done on a lattice with 80 grid points in the cos θ direction, and 40 grid points in the φ direction. The exact N c = 3 solution to (21) (solid red line) is compared 6 At cos θ α = ±1, U α and ξ α are independent of φ, as they should. In the case of U α , this is guaranteed by our initial condition (55) and the structure of the evolution (50) which preserves this property. 7 In the JIMWLK case, the initial condition is the value of the S-matrix (1) at small, but not too small value of x. This is model dependent and its initial sampling is non-trivial [16]. 8 We actually plot the real part of the average trU α U † β . trU α U † β is in general complex-valued for each trajectory, but we have checked that the imaginary part of the average is consistent with zero within errors. with the solution of the large-N c BMS equation (16) (dotted green line) previously obtained in [22], 9 and also with the solution of the 'mean field approximation' to (21) (dashed blue line) which differs from the BMS equation only by the coefficient of the Sudakov term N c = 3 ↔ 2C F = 8/3. The latter serves as an indicator of the quality of the mean field approximation P P → P P . For the sake of reference, we also plot the solution obtained by keeping only the Sudakov term (first term on the right-hand-side) in (58) (dash-dotted yellow line). 9 Note that the definition of τ in [22] differs from (17) by a factor of N c .
A comparison of the solid (red) and dashed (blue) lines in Fig. 2 shows that the exact solution is rather close to the mean field solution with the finite-N c corrected Sudakov term for all values of τ explored in this work. This may come as no surprise to those who are acquainted with the solution of the JIMWLK equation which agrees well with the BK solution [16,23]. Nevertheless, we find the present result quite intriguing because we actually observed enormous trajectory-by-trajectory fluctuations. In Fig. 3 we show 50 (out of 500) individual random walk trajectories used in the computation of the average. It turns out that the fluctuations are so large that the standard deviation δP ≡ P P − P P is of the order of P itself for not-sosmall values of τ . Actually, this is the reason why we needed to run O(100) trajectories to obtain a reasonably stable result. Such large fluctuations have not been seen in the previous simulation of the JIMWLK equation where 'already one trajectory gives a good estimate of the final result' [16].
In our opinion, the crucial difference between the two problems which has resulted in such different behaviors of fluctuations is the initial condition. In the jet problem, the initial condition P τ =0 = 1 means that there are no partons besides the qq pair. In the parlance of saturation physics, the system is initially very 'dilute'. The evolution then produces soft gluons which multiply exponentially in τ according to the BFKL formula [9,18]. It has been demonstrated that these gluons have very strong number fluctuations [24] and spatial correlations [25,26,19]. We thus find it natural to attribute the observed large values of δP to such fluctuations and correlations. On the other hand, when solving the JIMWLK equation, one often uses 'dense' or 'classically saturated' initial conditions; S xy τ =0 ≪ 1 for |x−y| larger than some value. Since there are many uncorrelated gluons in the system from the beginning, there is little room left for pure-BFKL evolution, i.e., it is suppressed by the saturation effect. Accordingly, fluctuations and correlations can develop only weakly, and this is consistent with what has been found in the previous simulations.
In conclusion, we have demonstrated for the first time the resummation of non-global logarithms at finite N c to leading logarithmic accuracy. Our study shows that, at least in e + e − annihilation and for phenomenologically interesting values of τ 0.2 ∼ 0.3, event-averaged non-global observables may be reliably computed in the mean field approximation by solving the (modified) BMS equation (58), or equivalently, by Monte Carlo simulations [6]. However, the observed large fluctuations imply that the situation may be drastically different for hadronhadron collisions (cf. [27]). Since four partons are involved in hard scattering, one has to deal with multiple products of Wilson lines such as tr(UU † )tr(UU † ) and tr(UU † UU † ) (cf. [28]). Moreover, if there are gluons in the initial and final states, each of them picks up an adjoint Wilson lineŨ which further increases the number of (fundamental) Wilson lines. [Roughly,Ũ acts like the square of U.] The proper treatment of fluctuations laid out in this paper is potentially very important for such observables. We leave this problem for future work.