The structure of IR divergences in celestial gluon amplitudes

The all-loop resummation of SU$(N)$ gauge theory amplitudes is known to factorize into an IR-divergent (soft and collinear) factor and a finite (hard) piece. The divergent factor is universal, whereas the hard function is a process-dependent quantity. We prove that this factorization persists for the corresponding celestial amplitudes. Moreover, the soft/collinear factor becomes a scalar correlator of the product of renormalized Wilson lines defined in terms of celestial data. Their effect on the hard amplitude is a shift in the scaling dimensions by an infinite amount, proportional to the cusp anomalous dimension. This leads us to conclude that the celestial-IR-safe gluon amplitude corresponds to a expectation value of operators dressed with Wilson line primaries. These results hold for finite $N$. In the large $N$ limit, we show that the soft/collinear correlator can be described in terms of vertex operators in a Coulomb gas of colored scalar primaries with nearest neighbor interactions. In the particular cases of four and five gluons in planar $\mathcal{N}=4$ SYM theory, where the hard factor is known to exponentiate, we establish that the Mellin transform converges in the UV thanks to the fact that the cusp anomalous dimension is a positive quantity. In other words, the very existence of the full celestial amplitude is owed to the positivity of the cusp anomalous dimension.


Introduction
Infrared (IR) divergences are ubiquitous in gauge theory scattering amplitudes. They arise when charged external states interact with an infinite amount of virtual excitations with arbitrarily low energies. The existence of these infinities suggests that the usual Fock space, constructed out of a unique vacuum, is not properly defined. Consequently, the notions of asymptotic one-particle states, and the vacuum itself, need to be revised.
Despite these issues, the most common strategy to deal with IR singularities is to focus on physical quantities such as the total cross section. Indeed, it was proven that IR-safe observables can be extracted from QED [1] and QCD [2,3] cross sections, by a careful cancellation between virtual and real emissions, order-by-order in perturbation theory.
Nonetheless, there is another way of dealing with infinities arising from these long-range interactions, namely, by extracting the IR-safe information from the amplitude itself. Due to the universal structure of the infrared radiation, it is possible to separate the different energy scales involved in the collision. Using the eikonal approximation, one chooses a scale Q where the finite contribution to the amplitude is consistently separated from the soft and collinear singularities. This type of organization in QCD has been the focus of extensive work [4][5][6][7][8][9][10].
The key idea behind the IR factorization of gauge theory amplitudes is that, since the energy scale is arbitrary, soft factors must be independent of Q. Then, one can use renormalization group techniques to treat the singular piece in terms RG flow equations. This approach imposes strong constraints on the IR structure of the amplitudes, resulting in exponentiated expressions valid at all orders in perturbation theory [11][12][13][14][15].
Similarly, a deep connection has been shown to exist between IR divergences in amplitudes and UV divergences in Wilson line correlators [16][17][18][19], permitting a gauge invariant characterization of these singularities. In particular, it is worth highlighting the work of Korchemsky and Radyushkin [20,21], who showed that the renormalization properties of cusped light-like Wilson loops universally describe the pole structure of the regularized gauge theory S-matrix. The quantity controlling these poles is known as the cusp anomalous dimension, γ(α), a perturbative function of the gauge theory coupling α.
Over the last years, there has been a renewed interest in understanding the role of IR divergences in gauge theory and gravity from a perspective based on symmetries [22][23][24][25][26][27][28][29][30][31][32]. The central point in this discussion is that the degenerated vacuum, observed in theories with massless excitations, is greatly explained by the presence of an infinite set of conserved charges.
These quantities correspond to surface integrals located at the null boundary of the spacetime that were first discovered in [33,34] and further developed in a holographic context in [23,35]. Their importance in this discussion is twofold: on the one hand they do not change the energy of the process, thus inducing an infinite degeneracy of the vacuum, but introduce degrees of freedom helping to account for the universal organization of IR radiation; on the other hand, their two-dimensional nature shows an intriguing holographic realization of generic 4D gauge theories.
Along these lines, celestial amplitudes provide an specific realization of a holographic description of flat spacetime physics. They are a map between 4D scattering amplitudes for massless particles and conformal correlators on the sphere at null infinity [36,37], known as the celestial sphere. The isomorphism between the Lorentz group and SL(2, C) identifies the usual plane-wave asymptotic states with conformal primaries on the celestial sphere whose CFT correlators are built upon [38]. Progress along these lines encompassing different perspectives in these topics can be found in .
Two-dimensional correlators are obtained through a Mellin transform of the bulk amplitude. More concretely, the Mellin transform involves an integration over all the energies of the external particles participating in the scattering process, from the deep infrared to the ultraviolet. Hence, their study may reveal yet another aspect of the UV/IR connection, which could help to further constrain the S-matrix. In other words, they may impose a new kind of high energy consistency relations onto scattering amplitudes.
In this article, we connect the singular behavior exhibited by all-loop IR regulated n-point gluon amplitudes with 2D conformal field theory structures. We follow and extend the approach recently presented in [69], showing that the factorization of divergences in non-abelian gauge theories [12,14,[16][17][18][19][20][21]70] persists in a universal fashion when expressed in the celestial basis. More concretely, we first identify the hard contribution as a celestial correlator with shifted conformal dimensions. This shift corresponds to an integral of the cusp anomalous dimension over the energy scales involved. This quantity becomes divergent once the IR-regulator is removed. Furthermore, we notice that the divergent factor, consisting of soft and collinear singularities, can be expressed as a correlator of colored primary fields on the celestial sphere, given by Wilson lines anchored to each external particle. By carefully analyzing the SL(2, C) covariant properties of these Wilson line operators, we find that their conformal weights have the exact same shift observed in the hard correlator. This leads us to conclude that the celestial-IR-safe gluon amplitude corresponds to a correlator of operators dressed with Wilson line primaries.
The organization of this paper is as follows. In section 2, using the BDS formula, we examine the celestial factorization of gluon amplitudes in the large N limit of maximally supersymmetric Yang-Mills [71]. In 2.4, the divergent piece is shown to be effectively described by a Coulomb gas of scalar colored primaries with nearest neighbour interactions. Moreover, in subsection 2.5, it is shown that for four and five gluons, given that the hard factor of these amplitudes also exponentiates, the convergence of the Mellin transform in the UV is ensured by the positivity of the cusp anomalous dimension. Section 3 generalizes the results of the previous section to SU(N) non-abelian gauge theory. We start by describing, using the color-space formalism [72], generic properties of all-loop IR regulated gluon amplitudes in momentum basis, based on [12,14,70]. Subsection 3.1 is devoted to describe IR divergences in terms of correlators of Wilson lines. In 3.2, we show that the hard celestial correlator is built out of operators with infinitely shifted conformal dimensions as shown in (3.15). At the end of this subsection, it is shown that this exact same shift arises as the anomalous dimension of renormalized Wilson lines. It is then argued that the anomalous dimension corresponds to the conformal weight of these operators. We conclude this section in 3.3, where the large N limit of SU(N) celestial gluon amplitudes is discussed, and we demonstrate that the structure of the divergences matches (2.17). We close this article with section 4, analyzing possible corrections to the results presented here and describing future work. Appendix A is devoted to show that, after a simple observation, the full expression for the celestial BDS formula of [56] it is directly mapped to our results of section (2.3).
Note added: While this paper was being finalized, overlapping results appeared in [73].

Review of the BDS Ansatz
We review known results of all-loop n-gluon amplitudes in maximally supersymmetric SU(N) Yang-Mills theory (SYM), for maximal helicity violation (MHV). Perturbatively, they can be expressed as with the Dirac delta enforcing momentum conservation and α = 2e −ǫγ E g 2 /(4π) 2−ǫ being the Yang-Mills coupling constant. The term A (L) is the L-loop contribution that is regularized by computing momentum integrals in D = 4 − 2ǫ dimensions. The amplitude A (L) is decomposed in single and multi-traces of generators in the fundamental representation of SU(N). In the large N limit, the leading contribution is controlled by single traces of n generators and reads where π runs over non-cyclic permutations of the external legs. The expression above helps us disentangle the color structure (given by the trace factors) from the dynamical content. This property can be used to examine the kinematics of the all-loop planar gluon amplitude with a specific color ordering. Hence, in the rest of this section, we focus on the ordering with the 't-Hooft coupling a beign now the perturbative expansion parameter. In 2005, Bern, Dixon and Smirnov [71], based on previous work [74], were able to show that the sum (2.3) actually exponentiates to where A tree is the color-ordered MHV tree-level amplitude [75] containing all the polarizationdependent information, while the exponential factor carries the full infrared-divergent structure as well as a finite contribution. Kinematically, the exponent above depends only on the Mandelstam invariants s i,j = 2p i · p j , with the identification s n,n+1 ≡ s n,1 , where p i is the momentum of each external gluon. 1 It is interesting to note that, due to the planar limit, interactions among neighboring external particles are the only ones contributing with the divergent terms in the 1 We work with the metric η = diag(+, −, −, −) expression above. We will investigate the implications of this organization in the celestial basis in subsection 2.4. Since N = 4 super Yang-Mills is a UV-finite theory [76,77], the remaining divergences are soft and collinear singularities originating from integrations over internal loops. In order to regulate these divergences one needs ǫ < 0, a condition that will be used throughout this article.
In expression (2.4), double and single poles in ǫ are controlled by the coefficients γ (l) and G (l) , consisting in the l-loop contributions to the cusp and collinear anomalous dimensions respectively, i.e., (2.5) These quantities arise as coefficients in the renormalization group equations for certain observables, such as Wilson loops and form factors. The cusp anomalous dimension γ(a) shows up in the renormalization of the product of two semi-infinite Wilson lines [16,17,21,78], while G(a) makes its appearance in the RG flow that defines the Sudakov form factor [79]. What makes the formula (2.4) particularly appealing is the appearance of the full finite part of the amplitude, denoted by F n . For n = 4 and n = 5, F n is directly expressible in terms of one-loop data, however, starting at n = 6, it also acquires higher loop corrections depending on dual conformal invariant cross-ratios of the Mandelstam variables [80][81][82][83].
In what follows, we will make use of the celestial basis [38] to interpret the role played by the divergences in (2.4) as the vacuum expectation value of operators in a two-dimensional CFT. Furthermore, these operators will become a crucial ingredient when extracting the infrared-safe information from (2.4) and translating it to the celestial sphere.

IR divergences in celestial n-point correlators
The momentum of each external massless particle in an scattering process will be parameterized in terms of complex coordinates (z i ,z i ) and a real number ω i as where η i is ±1 for outgoing/incoming particles. In terms of the above decomposition, the Mandelstam variables read with z ij = z i − z j . We now consider the Mellin transform on each external leg of the BDS amplitudeM As it has been shown in [36][37][38], the main property of a Mellin transformed amplitude is its covariance under SL(2, C) transformations. Since this corresponds to the global part of the conformal group in two dimensions, one then expects that (2.8) can be expressed as a correlator of n insertions on the celestial spherẽ where one can identify the conformal primary with ℓ i being the helicity of each external particle. Note that for the MHV amplitudes considered here, all but two of the spins ℓ i are positive.
Whenever infrared divergences are present in an amplitude, it is clear that (2.9) is not a well-defined object as it will inherit these singularities when taking ǫ → 0. However, in order to define a celestial correlator free of divergences, we will show that soft and hard gluonic degrees of freedom decouple in the celestial basis. In fact, it has been observed in [69] that for gravity and scalar electrodynamics, celestial amplitudes can be arranged to be written as where the first factor controls the IR (soft) divergences through a correlator of n operators V σ i (z i ,z i ). These are primaries with zero spin and conformal dimension 2σ i that are both divergent and regulator-dependent. Based on this factorization, one defines the infrared-safe part of the celestial amplitude as the second factor in (2.10). This corresponds to the expectation value of n dressed operatorŝ O ∆ i ,ℓ i with conformal dimension ∆ i − σ i . It is then appealing to suggest, for gluon amplitudes, the normal-ordered relation that has, in fact, already appeared in [84,85] for electrons and gravitons, respectively. In the next section we focus on the factorization (2.10) and we reveal the role of (2.11) in the case of celestial BDS amplitudes.

BDS celestial amplitude
We would like to start by investigating the divergences of the celestial amplitude (2.8). To do so, let us expand the divergent part of the exponent in (2.4) in powers of ǫ with the understanding that γ 0 (a) = γ(a) and G 0 (a) = G(a) as defined in (2.5). To perform the Mellin integrals in (2.8), we first split the exponential in divergent and finite terms. We then promote the Mandelstam variables to operators acting on conformal primaries. As we will see below, the use of operators of this kind streamlines the expression for the hard celestial amplitude. LetP i be an operator such that [69] is the corresponding primary field defined on the celestial sphere. Acting on a celestial amplitude, one can useP i to define the Mandelstam operatorŝ (2.15) whose eigenvalues on plane waves are the Mandelstam variables.
Let us now go back to the Mellin integral (2.8). From (2.12), we see that the full celestial amplitude factorizes into a divergent part times a finite one, 2 i.e., where, using the cyclic property δ 1,n+1 ≡ δ 1,1 = 1, we have organized the divergent part as 4ǫ being a non-dynamical factor depending on the coupling only. The finite contribution becomes an operator acting on the tree-level celestial amplitudeÃ tree , namelỹ The hard amplitude in (2.18) has been defined in terms of shifted conformal dimensions ∆ ′ i = ∆ i + γ 1 4ǫ . Note that even though the dimensionful factor (µ 2 ) − γ 1 8ǫ blows up as ǫ → 0 − , the full object defined in (2.18) will be finite after using the explicit expression for F n (a,ŝ ij ). In fact, the exponential operator appearing in (2.18) simply enters through the integrals, transforming thê s i,i+1 operators into the Mandelstam variables s i,i+1 , thus allowing us to explicitly evaluate the full expression for the finite amplitudeM hard . This will be the focus of subsection 2.5. However, to properly define (2.18) as a celestial correlator in terms of dressed operators as in (2.11), we first need to examine in more detail the conformal structure of the divergent piece (2.17).

Divergent part as a conformal correlator
The first relevant property of (2.17) is that it transforms as a n-point celestial correlator in the sense that it is covariant under SL(2,C). From the transformation 20) thus yielding that the divergent part of the celestial amplitude behaves as a conformal correlator composed of spin-zero primary fields V σ (z i ,z i ) with weights Note that σ is positive due to the positivity of γ 1 (a). This follows from the fact that the cusp anomalous dimension γ is positive for arbitrary values of the Yang-Mills coupling α, [86] 3 . Then, one first concludes that γ 1 (a) is a monotonically increasing function of the coupling, but since γ 1 (0) = 0, the announced positivity is ensured. Another appealing feature is the possibility of finding an explicit representation for the operators V σ i (z i ,z i ). For such purpose, we introduce scalar fields Φ a (z,z) with an extra color index a labelling its m entries. We define the two-point function of these scalars as Let {e 1 , · · · , e n } be a set of n ≤ m vectors in R m satisfying the orthornomality condition e a i e a j ≡ e i · e j = δ i,j . We use this basis of vectors to represent the interactions between nearest neighbors appearing in the divergent part of the BDS amplitude. For each external particle located at the position (z k ,z k ), we define vectors τ a k = e a k − e a k+1 with the identification e n+1 ≡ e 1 . This is an overcomplete set, as one can verify that the total sum vanishes, i. e., n k=1 τ a k = 0 . (2.23) This constraint suggests a conservation law. We will later identify the vectors τ a k with the generators T a i in the adjoint, representing the color charge of each of the external gluons. The expression above thus simply reflects color conservation. We will briefly expand on this in subsection 3.3 where we will also identify the dimension m.
For the purpose of expressing the divergent amplitude as an explicit celestial correlator, it is useful to rewrite the sum in (2.17) in terms of the inner product 4 (2.24) One then considers vertex operators defined by Using the two-point function (2.22) and (2.24), one can show that the divergent part of the amplitude can be explicitly written as a correlator on the celestial sphere as Furthermore, the importance of identifying the explicit representation (2.25), is that we are now able to construct dressed operators controlling the hard celestial amplitude following the discussion of section (2.2). In particular, due to the OPE it is possible to invert relation (2.11) to obtain the corresponding hard operator, where we can clearly see that the role of V −τ i is to shift the conformal dimension by −σ i . Also notice that, except for the jet function, this expression coincides with previous results for celestial amplitudes in scalar QED [69]. Thus, the hard BDS celestial amplitude is defined through It is worth stressing here that, since J and the operators V τ i are regulator dependent, they are not measurable quantities. However, ratios of divergent correlators, like the hard amplitude defined above, are finite and give rise to IR-safe observables. In the next subsection, we show that this expression precisely matches with the computation of the Mellin integral (2.18) for four and five gluons.

Finite part: The UV and the positivity of the cusp anomalous dimension
We start by evaluating the finite part of celestial factorization in (2.18) for the case of four and five gluons. The fact that this is even possible is thanks to the main feature of the BDS formula, namely, that the infrared-finite part of the all-loop amplitude exponentiates as well. As a by-product, we will see how the positivity of the cusp anomalous dimension γ(a) is directly responsible for the existence of the corresponding hard CFT 2 correlation function on the celestial sphere.

Four gluons
Taking the Mellin transform on each external particle, the explicit form for the finite part of the four-gluon celestial amplitude is Here, r is the conformally invariant cross-ratio and C(a) is a function of the coupling only. 5 The integral above can be performed along the lines of [56] (see Section 2.2) where, for definiteness, we assume two ingoing (η 1 = η 2 = −1) and outgoing particles (η 3 = η 4 = 1). In terms of the shifted conformal weights we obtain a four-point CFT 2 correlation function with where f ∆ ′ i is a conformally invariant factor, having a distributional nature in the r variable 6 , given by After the change u = log(w/µ 2 ), one immediately sees that this is a Gaussian integral whose convergence in the UV depends crucially on the positivity of the cusp anomalous dimension γ(a), a fact that has been shown to hold to all orders in perturbation theory [86,87]. Therefore, turning the argument around, we can venture to say that the very existence of the celestial amplitude, corresponding to the planar gluon amplitude in N = 4 SYM, to all orders in the 't Hooft coupling, is directly related with the positivity of the cusp anomalous dimension. Evaluating the gaussian integral above yields As a consistency check, with ǫ held fixed and taking the zero-coupling limit a → 0, we should recover the tree level celestial amplitude. More precisely, we must obtain the same expression as in (2.33) but with bare weights (h i ,h i ) and with the conformally invariant function Indeed, note that (2.35) is a gaussian function since, in the a → 0 limit with ǫ held fixed, β becomes a pure imaginary number. From the fact that γ(a) = O(a) and using the representation of the Dirac delta 7 Since we now have an explicit expression for the full celestial correlation function corresponding to the all-loop four gluon (planar) amplitude in the bulk, we must make sure that, after computing its inverse Mellin transform, we arrive at the original BDS formula in momentum space. In particular, since this process involves integrating over β from −i∞ to +i∞, the correlator better have a good behavior in the large |β| region in order assure convergence. But since (2.35) is a gaussian function in iβ, the inverse Mellin transform, indeed, converges rapidly. 6 The precise form of the prefactor I ∆ ′ i (r,r) is somewhat irrelevant for our purposes here, but it is given by 6 γ(a)+C(a)+O(ǫ) 7 In taking the limit, we have also used that G(a) and γ 1 (a) are of O(a).

Five gluons
In the four-particle case, because the scattering takes place on a plane, the four-dimensional Dirac delta enforcing total-momentum conservation yields an overall factor δ(r −r) in front of the amplitude, i.e., the cross-ratio r in (2.31) is forced to be real. The rest of the factors from the delta function localize three of the four energies ω i , yielding a single Mellin integration left. For five-point celestial amplitudes, the 4D Dirac delta fully constrains four of the energies, which can be chosen to be ω * i for i = 1, . . . , 4. Then, one can write [46], where r 4 is one of the two conformal invariant cross-ratios defined by 8 In the case of five gluons, the four energies ω i , i = 1, . . . , 4, become localized at where the coefficients f i5 are functions that depend only on z ij and the cross-ratios r i . The full expression for the 5-gluon celestial amplitude, for all values of the 't Hooft coupling, isM where the function h 5 , entering in the exponent above, is given by 9 Before continuing, we would like to notice the following. All of the helicity structure of the original amplitude in momentum basis is contained in the tree-level factor M tree 5 . This, together 8 The relation between the cross-ratio r of the four-gluon (previous subsection) and the r 4 ratio defined here, is r 4 = r r−1 . 9 Here we used the expressions in section 4.C of [71].
with the rest of the |z i,j | dependence in the prefactors above, nicely combine to produce the fivepoint CFT 2 correlation function, with the correct covariance under SL(2,C) transformations. Therefore, since the energies ω * i do transform under SL(2,C), one may wonder whether the explicit dependence on the ω * i in the expression for h 5 in (2.43) could spoil the aforementioned covariance. The answer is of course not, since from (2.41) and the explicit expressions for f ij (see Appendix A of [46]) one notices that they transform in such a way that the final evaluated expression for h 5 ({ω * i }, ω 5 ) ends up being conformal invariant.
where the coefficients b(r i ) and c(r i ) are functions of the conformally invariant cross-ratios only. After making the change u = log(ω 5 /µ), one sees that the integral in (2.42) is again a gaussian whose convergence is assured by the positivity of γ(a).

Infrared divergences in gluon amplitudes
In order to extend the results of the previous section, we first describe the general structure of IR divergences of scattering amplitudes in SU(N) gauge theory. Throughout this section we use the color-space formalism of Catani and Seymour [72,90]. In this construction, contributions to amplitudes are vectors decomposed in the space of product of traces of n generators t a in the fundamental representation of su(N). Tree level n-gluon amplitudes are decomposed in terms of single traces. At loop level, apart from single traces, the decomposition also allows for the appearance of multi-trace contributions (for a review see e. g. [91]). It is convenient to introduce the color charge operator T i = {T a i }, which is a vector associated with the i-th external gluon, with N 2 − 1 entries. It is defined in the adjoint representation, therefore, it acts as a commutator on a generator in the fundamental representation Color charge is conserved in a scattering process. Therefore, when acting on amplitudes, the sum over all operators T a i vanishes We shall need the inner product between color charges T i · T j = T a i T a j . The Casimir operator C i = T 2 i is proportional to the identity and thus commutes with the color operators. By using dimensional regularization in D = 4 − 2ǫ dimensions, the all-loop n-gluon amplitude resummation exhibits the factorization property [12][13][14][15], where the path-ordered exponential Z contains all the information about infrared divergences. The hard part of the amplitude, M hard corresponds to a vector in color-space, and it is finite in the ǫ → 0 limit. The argument of the exponential in Z is integrated up to an arbitrary energy scale µ, and α(λ 2 ) is the coupling that is generated through the renormalization group equation where β(α) is the beta function of the theory. We consider the soft anomalous dimension matrix Γ n given by where the first term contains contributions from pairwise interactions. The proportionality constantγ(α) = γ(α)/C A is the cusp anomalous dimension divided by the Casimir in the adjoint representation. 10 The second contribution in (3.5) is a sum over the collinear anomalous dimensions associated with each gluon γ J i , and the last term ∆ n is a non-vanishing contribution for n ≥ 4 that starts at third order in the coupling α [94], The kinematic dependence of ∆ n is only through the cross-ratios ρ ijkl (3.7) From the second relation above one immediately sees that the cross-ratios ρ ijkl are also invariant under SL(2, C) transformations [94,95]. Hence, it is a conformally invariant function on the celestial sphere. Notice that it was recently observed that (3.5) receives extra contributions starting at four loops [70]. Unlike ∆ n , these extra terms are not constrained to depend on the cross-ratios ρ ijkl only, but rather depend explicitly on logarithms of the Mandelstam invariants s ij . Fortunately, precisely thanks to this logarithmic behavior (and color conservation), we will see that the effect of including these four-loop corrections are under control in terms of the celestial correlation functions (See Section 4 for more details).

Infrared divergences as a correlator of Wilson lines
The renormalization factor Z(α, {s ij }) in (3.3) can be expressed in terms of a correlator of the product of n light-like Wilson lines of the form [19,96] where A µ = A a µ T a i is the gluon field for each external particle with momentum p i . In order to see this, we first note that the renormalization factor can be rewritten as Notice that both γ ij and the exponential prefactor in (3.9) are ill-defined expressions since p 2 i = 0 for external gluons. Therefore, we shall think of them by first giving a small mass to each external gluon, and then taking the massless limit at the very end. When doing this, one obtains [14] Z(α(µ 2 ), {s ij }) = lim which is free of the aforementioned light-like divergences. However, notice that a different type of singularity will remain, namely, UV divergences in the Wilson lines arising from the presence of cusps. These infinities are in a one-to-one correspondence with the infrared structure of the full amplitude (3.3), and whose role on the celestial sphere will be the focus of the following subsection.

Celestial gluon amplitudes
For the purpose of finding the celestial amplitude of (3.3), it is convenient to write the Mandelstam variables as in (2.7) s ij = η i η j ω i ω j |z ij | 2 . Using this and the color conservation condition (3.2), one deduces that Γ n decomposes as (3.12) 11 In this expression we will actually need to introduce a damping factor in order to render the correlators finite. For a recent review on these developments see [97] and references therein.
The above separation of the angle-dependent part z ij and the energies of the external particles ω k has an important consequence on celestial amplitudes. More precisely, it factorizes the renormalization factor as Then, taking the Mellin transform of the amplitude (3.3), we obtaiñ whereM hard is the Mellin transform of the hard amplitude M hard with the µ 1 2K k C k prefactor from (3.13) included in its definition. Notice that this is the same type of factorization obtained in the planar limit N = 4 super-Yang-Mills in (2.16). However, this result is valid for any gauge theory (even without supersymmetry) because it follows from the general factorization displayed in (3.3).
The distinctive property of the celestial correlatorM n , is that it affects the hard factor by a mere shift in the conformal dimensions ∆ i of its primaries in the amount As the dimensional regulator ǫ is removed,K diverges, corresponding to an infinite shift, precisely matching with the one obtained for the case of planar N = 4 SYM theory.
Recall that the factor Z(α, µ 2 |z ij | 2 ) in (3.14) governs the infrared behavior of the full amplitude and it possesses the color structure of a SU(N) tensor product of Wilson lines in the adjoint representation. Among its appealing features are that it is solely expressed in terms of the celestial data z ij and transforms covariantly under SL(2, C). Furthermore, we can use the regulated version of (3.11) with p i replaced by µq i to express this factor in terms of purely celestial variables as where W q i (0, ∞) is a semi-infinite Wilson line in the regulated direction q i ≡ q(z i ) given by satisfying q 2 i = δ 2 with δ a dimensionless parameter. At this point, one may think that W q i can be regarded as a primary field on the celestial sphere. For this to be true, we would need to analyze the transformation properties of these Wilson line operators. At leading order in the regulator δ, the action of SL(2, C) on q i produces a simultaneous Lorentz transformation and a scaling, i.e., On the other hand, since Wilson lines on the half real line do not transform under rescalings q i → λq i , for all λ, we conclude that these operators are SL(2, C) invariant in the δ → 0 limit. Hence, they do not display the needed transformation laws accounting for the covariant properties of Z(α, µ 2 |z ij | 2 ) as a celestial correlator. However, due to the nature of the scattering process, we will see that W q i is a bare operator that suffers from UV divergences associated to the presence of cusp singularities. As such, its transformation properties after renormalization will become anomalous producing the desired covariance under SL(2, C). Let us consider the simplest observable containing only two semi-infinite Wilson lines 12 Notice that the second operator in the equation above has the reversed order due to color charge conservation T ≡ T 1 = −T 2 . This fact allows us to interpret this correlator as a Wilson loop that closes smoothly at infinity, but has a cusp at the origin. This type of singularity produces the following renormalization group flow [16,17] d log(W R ) d log(λ) = −Γ cusp (γ 12 , α(λ)) , (3.20) where Γ cusp is the anomalous dimension of the renormalized two-point function W R in the presence of a cusp singularity. To make contact with light-like Wilson lines, we study the properties of the equation above in the δ → 0 limit. Remarkably, it has been found in [21] that for all values of the YM coupling α, the cusp anomalous dimension admits the following expansion with C = T 2 being the Casimir of either colored state, and where the O(1) terms represent contributions that do not depend on the cusp angle, i.e., they are independent of the 'distance' |z 12 | on the celestial sphere. This fact, together with the logarithmic behavior above, yields that under the global conformal transformations (2.19), the correlator W R does transform covariantly. Therefore, we shall identify the renormalized version of W q i (0, ∞) as the primary operator of interest on a celestial CFT 2 . As a consequence, its anomalous dimension becomes the conformal weight σ i of the spinless operator Summarizing, the correlator of two semi-infinite light-like Wilson lines can be obtained by integrating (3.20) and using the limit (3.21), where J = J (α, ǫ) is, until now, an undetermined prefactor, while the Kronecker delta enforces color charge conservation. It is interesting to highlight the resemblance between the relation just found and the OPE involving operators V τ in the case of planar SYM amplitudes (2.27). More precisely, an examination of the soft anomalous dimension (3.5) in the case of two external gluons, [12,13] reveals that the prefactor J is given by where each external gluon contributes with the same collinear dimension γ J i = γ J . Interestingly enough, this quantity coincides with the jet factor defined for the first time in (2.17) in the case of SYM planar amplitudes.
Having identified the Wilson line primary field as the responsible for the infrared divergences, one can extend the application of (2.11) to the color space basis as that enables us to properly define the colored celestial hard amplitudeM hard as a correlator of IR-safe operatorsÔ ∆,ℓ (z,z).
In the next subsection we study the large N expansion of (3.14), arriving at some general results for gluon amplitudes in planar QCD. As a by-product, for the specific case of N = 4 SYM, we will reproduce all of our findings already obtained in section 2.

Large N contributions
We describe how the results obtained in section 2 can now be reproduced in the general framework of infrared divergences in non-abelian gauge theories. Let us first notice that, in the large N limit, the color structures appearing in (3.5) project the amplitude onto vectors in color space with single-traces only. In fact, the full amplitude in the planar limit can be reconstructed from n (square roots of) Sudakov wedges between adjacent external particles [79,98]. In particular, this argument implies that the higher loop order contributions to the soft anomalous dimension (3.5), must vanish, i.e., ∆ planar n = 0 . (3.26) Thus, the only object that needs to be studied in more detail is the action of T i · T j on single traces. This has been analyzed in appendix B of [15], and here we briefly review their argument.
When acting on single traces of n generators in the fundamental representation t a i , with fixed color structure, one has T i · T j tr[t a 1 · · · t an ] = tr[t a 1 · · · [t a i , t a ] · · · [t a j , t a ] · · · t an ] . (3.27) Expressions such as the one above simplify greatly with the use of the identity 13 involving the product of k generators. Then, one notices that the leading N contributions in (3.27) will only arise when the precise combination t a t a is present or, in other words, when i and j are nearest neighbors. The operator T i ·T j then effectively acts on single traces of n generators, for i = j, as We can therefore conclude that, on the celestial sphere, infrared divergences in the large N limit are fully encoded in the expression withK defined in (3.10) and J given by relation (3.24) where we have set γ J i = γ J , which is true for external particles of the same kind, such a gluons. It is worth noticing that the expression for Z planar applies to any non-abelian gauge theory in the large N limit. Therefore, the discussion presented in section 2.4 which allowed for the identification of the celestial operators V τ i , also holds here provided one recognizes 1 4 NK as the conformal weight σ of the primary field (2.25). Furthermore, the color operator T i is effectively replaced by the vector τ i in R m . Since T i is in the adjoint representation, we notice the number of entries m scales as m ∼ N 2 in the large N limit.
Finally, the N = 4 SYM case can be immediately recovered by recalling that beta function vanishes. Then, from (3.4), one has with α 0 is the coupling at the energy µ. We then perform explicitly the integration over the energy scale λ in (3.3) for ǫ < 0. Using a series expansion in the coupling α, the coefficientsK and J are then given by with the identification G (l) = 2γ (l) J . We have also used that, in the planar limit, Nγ m (α) ≡ γ m (a) and G m (α) ≡ G m (a) with γ m and G m defined in (2.13). Hence, the two quantities above show that infrared divergent contribution Z planar exactly coincides with the factorization given in section 2.3.

Final remarks
We have described the celestial representation of infrared corrected gluon amplitudes using the well-established results [12,14,70,71] in momentum basis. We have first analyzed the planar case and then we have discussed finite N contributions, where the crucial object is the soft anomalous dimension matrix Γ n , (3.5). However, it has recently been found that this object gets corrected at fourth loop order [70]. Therefore, one may ponder on the precise effect this correction has on the celestial side. This has been emphasized in [73] and here we would like to elaborate on it.
The four-loop correction to Γ n is [70] ∆Γ n (α(λ 2 ), {s ij }) = −g(α) i =j where d abcd is a symmetric invariant tensor given in terms of traces of symmetrized products of generators in the adjoint representation. That the kimematical dependence of this correction is logarithmic in the Mandelstam variables s ij is a very welcome feature, since it implies that its effect on the celestial CFT (3.14) consists in a correction to the conformal weight given by with (compare with (3.10)) K corrected (α) = 1 2 4) where N i is the dimension of the representation (for gluons, N i = N 2 − 1). The second term inside the parentheses is precisely the four-loop correction to the cusp anomalous dimension found in [70,93]. In this sense, our proposal for the celestial vertex operators (3.22) includes this correction after making the simple changê which shows the robustness of our results. Moreover, including these corrections, the analysis presented in this work can be straightforwardly extended to massless external fermions by simply using now the operator T i in the fundamental representation.
There are other possible extensions of this work that would be interesting to pursue in more detail in the future. Here we summarize some of them.
One of the most interesting aspects of the BDS formula is the all-loop exponentiation of the finite part, a feature due to an extra symmetry enjoyed by the theory: dual conformal invariance. This symmetry is particularly manifest for the case of n ≥ 6 particles, where the finite pieces kinematically depend, exclusively, on dual conformal cross-ratios. This is a consequence of the close connection between dual conformal symmetry and finite part of the BDS amplitude. It then seems worthwhile to investigate what dual conformal invariance can teach us about the deeper structures of the corresponding celestial CFT.
Another interesting aspect has been emphasized in [85] for graviton amplitudes. They have shown that there is a link between the resummation of IR divergences in graviton amplitudes and memory effects. Since these effects are a way to retrieve IR-safe information from the amplitude, it would be worth addressing this link in the present context under the light of the color memory [99]. the operator P is defined bŷ P = µ 2 r Recalling the operatorŝ ij = η i η jPiPj |z ij | 2 and using momentum conservation for two ingoing (η 1 = η 2 = −1) and two outgoing particles (η 3 = η 4 = 1), we noticê s 12M =ŝ 34M ≡ŝM ,ŝ 14M =ŝ 23M ≡ −tM , (A.11) that permits us to find Therefore, upon acting with the Mellin transforms on each external state in (2.4), one arrives at the exact same expression (A.12).