Celestial IR divergences in general most-subleading-color gluon and gravity amplitudes

Gluon amplitudes at most-subleading order in the 1/N expansion share a remarkable simplicity with graviton amplitudes: collinear divergences are completely absent in both and, as a consequence, their full IR behavior arises from soft gluon/graviton exchange among the external states. In this paper we study the effect of all-loop IR divergences of celestial most-subleading color gluon amplitudes and their similarities with the celestial gravity case. In particular, a simple celestial exponentiation formula for the dipole part can be written. We also analize how this exponentiation is modified by non-dipole contributions. Finally we also show that, in the Regge limit, the soft factor satisfies the Knizhnik-Zamolodchikov equation hinting at the possibility that, in this limit, an effective Wess-Zumino-Witten model would describe the dynamics of the infrared sector.


Introduction
Celestial amplitudes offer an explicit realization of flat space holography. Momentum space amplitudes for massless external states are mapped to correlation functions on the two-dimensional celestial sphere at the null boundary of Minkowski space. After a Mellin transform over the energies of all the external particles participating in the scattering event, the Lorentz group SL(2,C) identifies the plane-wave asymptotic states with conformal primary fields defined on the celestial sphere at null infinity [1][2][3]. Celestial amplitudes thus exhibit the usual tranformation properties of correlation functions of a two-dimensional CFT, however, the precise underlying conformal field theory (or theories) describing these correlators is still largely unknown.
Momentum space amplitudes in gravity and gauge theories are plagued with infrared divergences. These infinities appear due to the interaction of external particles with an infinite amount of photons, gluons or gravitons carrying arbitrarily low energies. Remarkably, the entire effect of infrared divergences can be fully traced yielding a factorized form for the n-point amplitude, namely

JHEP01(2022)136
where the entirety of the infrared divergent part is contained in the universal soft factor Z n , whereas A hard is a process-dependent (and infrared-safe) n-particle amplitude. It has recently been shown that this factorization persists for celestial amplitudes [4][5][6], with Z being a function that depends solely on the celestial coordinates z i . For example, the Mandelstam invariants s ij get mapped to "squared distances" |z ij | 2 on the celestial sphere. Thus, being this factor a universal quantity depending only on celestial data, it appears to be a perfect candidate to describe general features of the precise but yet undiscovered celestial CFT (CCFT). In particular, one would like, at least approximately, to answer the question: what is the holographic two-dimensional CFT living on the celestial sphere at null infinity governing the infrared degrees of freedom of the four-dimensional gauge theory in the bulk?
Up until now, many steps have been taken in this direction. In [7] a theory of free scalars was put forward in order to describe the infrared sector of abelian gauge theory and graviton amplitudes. In [6,8] an analogous theory of free colored scalars was proposed in order to reproduce the celestial correlation function given by the celestial soft factor Z in nonabelian gauge theory. While succesfully describing the infrared sectors in terms of a 2-dimensional description of free scalar fields, we know that this is far from being the end of the story, since nonabelian gauge theory is much richer and subtle. For instance in [8], the 2dimensional theory takes into account the theory at finite N , but considering only the dipole contributions to the soft anomalous dimension matrix. In [6], all (known) contributions to the soft anomalous dimension beyond the dipole approximation were considered, but only in the large N limit. In both these cases, and in the abelian construction in [7], one may suspect that the absence of the effects from collinear divergences may be responsible for the apparent simplicity of the holographic description obtained [9]. Therefore, we believe that a richer 2-dimensional theory, beyond free colored-scalars, must be at play governing the dynamics of the infrared sector of nonabelian gauge theory. 1 In this article we take a step in this direction by focusing on the sub-sector of gluon amplitudes called most-subleading-color amplitudes [11][12][13][14][15][16]. These amplitudes are a specific (infinite) subset of loop contributions to the full amplitude that are most-subleading in the 1/N expansion at each order in the strong coupling constant. One of its remarkable features is that, similar to graviton amplitudes, these amplitudes are free of collinear divergences. This makes its infrared structure remarkably simple, while still being at finite N . In fact, within the dipole approximation, the infrared divergences are one-loop exact, and can be explicitly computed for a general SU(N ) gauge theory [16].
The similarity between graviton and most-subleading-color gluon amplitudes and their infrared structure is not a coincidence. There are a number of relations between them: exact relations for the 4-point functions at L = 0, 1, 2 loops, valid not only for the maximal supersymmetric cases (N = 8 vs. N = 4) [11], but also for arbitrary N [15], exact relations at one-loop 5-points extended from maximal supersymmetry [13] to general N [15], and relations between the first two IR-divergent terms in the maximal supersymmetry case [16].

JHEP01(2022)136
It is important to note that these relations among graviton and most-subleading-color gluon amplitudes hold for the full amplitudes, that is, soft and hard parts included. However, up until now, these connections are only available in a loop-by-loop analysis. This means that, contrary to the soft factor whose basic structure is known to all-loop orders, we do not yet have the luxury of making a statement between the celestial correlators for the full amplitudes of these theories. With this is mind, in this article we perform a loop order analysis for the celestial loop amplitudes similar to that of planar N = 4 SYM performed in [4].
We will also see that their infrared properties, being similar to those of gravity, makes writing down the celestial factorization of infrared divergences almost trivial in this case. Finally, we will focus on the celestial soft factor and show that, in the forward limit (Regge limit), it satisfies the Knizhnik-Zamolodchikov equation, signalling a WZW model as an effective theory of the infrared degrees of freedom in the Regge limit. This paper is organized as follows. In section 2 we review celestial loop amplitudes and their factorization properties. In section 3 we consider celestial graviton amplitudes and their infrared divergent structure, with a view towards their relation with most-subleading gluon amplitudes. In section 4 we focus on the case of most-subleading gluon amplitudes in SU(N ) gauge theory, studying their infrared factorization properties and their celestial avatar. We will restrict ourselves to dipole contributions only in the soft anomalous dimension matrix. Then, in section 5, we will see how the Knizhnik-Zamolodchikov equation emerges, in the Regge limit, as a consistency condition for the soft factor in the mostsubleading sector. Finally, in section 6 we analyze corrections to our results beyond the dipole approximation.

Celestial amplitudes, N = 4 SYM and IR divergences
In this section we review the general results obtained in [4] for the exact 4-point celestial amplitude with massless external states, and the structure of infrared divergences at arbitrary loop order in the planar limit of N = 4 SYM theory.

Celestial amplitudes
We consider scattering amplitudes in four dimensions 2 with external massless particles. Each external state is labelled by a momentum p µ i , a helicity i , and a sign η i = ±1 distinguishing incoming (−1) from (+1) outgoing states. We parameterize p µ i in terms of points (z i ,z i ) on the two-dimensional celestial sphere through the map with ω i ≥ 0.

Basis for 4-point amplitudes
With the external particles written in terms of their helicities i , momenta p i as in (2.1), and stripping off the overall momentum-conserving delta function, the momentum-space amplitude is denoted as It is also convenient to further split the stripped amplitude as where B(s, t) is a function of two of the Mandelstam variables s = (p 1 + p 2 ) 2 , t = (p 1 − p 4 ) 2 , and u = (p 1 − p 3 ) 2 (which is not necessarily an analytic function), and the full dependence on the helicities of the external particles is contained in the rational function R ( 1 , 2 , 3 , 4 ) . The latter can be written in terms of the usual spinor helicity products ij and [ij] which, rewritten in terms of the complex variables (z,z) on the celestial sphere and the energy ω, can also have an explicit dependence on the SL(2,C) invariant cross-ratio which, using total momentum conservation, is directly related to the scattering angle θ cm in the center-of-mass frame through from where we see that r is real and r > 1.
The Regge (or forward) limit corresponds to taking large r and we will particularly focus in this regime in section 5. Also, although from (2.5) it may seem that r is complex, total momentum conservation forces r to be real (see, for instance, equation (2.14)), which is simply a consequence of the fact that the four-point amplitude always occurs on a plane.
As all the information about the helicities i of the external states is encoded in the function R ( 1 , 2 , 3 , 4 ) , the function must be an eigenfunction of theˆ i , whereˆ i is the helicity operator for the i-th particle, We may express a solution to the eigenfunction constraint (2.7) in terms of powers of Lorentz invariant functions R i , defined by the relationˆ i R j = δ ij R j , as (2.11) Notice that no mention whatsoever has been made to perturbation theory, thus, the statements made so far apply to the exact (all-loop) four-point S-matrix element.

General 4-point celestial amplitude
Celestial amplitudes are defined through a Mellin transform over the energies ω i of each of the n external states in a momentum-space amplitude A ({ω i , i , z i ,z i }), maintaining the variables (z i ,z i ) as coordinates on the celestial sphere, and reinterpreting the helicities l i as spins J i of operators, (2.12) After performing the transform, ∆ j are the conformal dimensions and J j ≡ j the spins of conformal primary operators inserted at the points (z j ,z j ) ∈ S 2 . Indeed, under SL(2, C) Lorentz transformations on the conformal sphere in (z,z) coordinates, the celestial amplitude (2.12) transforms as [2,3,18] A ∆ j , J i ; (2.13) i.e., it has the correct transformation properties of a 2D CFT correlation function defined on the celestial sphere. Specializing to the 4-point case, inserting (2.2) with (2.3) into (2.12), and using the delta function expressed in terms of the celestial sphere,

JHEP01(2022)136
we obtain the four-point correlator (up to a numerical factor) Here r is the conformally invariant cross-ratio (2.5) which, a priori, is a complex number. However, because of the delta function δ(r −r), the cross-ratio becomes real and equal to r = −s/t due to total momentum conservation. The conformally invariant function f (r,r) is given by (2.15) are given in terms of scaling dimensions and conformal spins by where the scaling dimensions ∆ k are restricted to lie on the principal continuous series of the SL(2, C) Lorentz group [3], i.e., As an example, consider the celestial tree-level four-gluon amplitude in pure Yang-Mills 4 first discussed in [2,18]. The stripped MHV four-gluon amplitude is where in the last step, the constraints from momentum conservation (2.14) were used. For the helicities 3 = 4 = − 1 = − 2 = 1, from (2.11) we read off that α 1 = 1, α 2 = 0, and the arbitrary function B(s, t) = g 2 is simply the coupling constant squared. 5 After applying the Mellin transforms in (2.12), one obtains the four-point correlation function (2.15) with conformal weights (h k ,h k ) = ( i 2 λ k , 1 + i 2 λ k ) for the negative helicity gluons and (h k ,h k ) = (1+ i 2 λ k , i 2 λ k ) for the positive helicity ones. The conformally invariant factor f (r,r) is where with λ = 4 k=1 λ k being the imaginary part of ∆. Enforcing the delta function, (2.21) yields ∆ = 4 giving f tree (r,r) = 8πg 2 δ(r −r)Θ(r − 1)r  (2.22) 4 At tree level, the pure Yang-Mills gluon amplitude is the same amplitude as in any other N supersymmetric Yang-Mills theory, such as N = 4 SYM which is one of the cases of interest in this article. 5 Hereafter we omit an overall factor of i(2π) 4 .

Loop amplitudes in N = 4 SYM and IR divergences
Here we focus on gluon amplitudes in the planar limit of N = 4 SYM. Because this is a conformal field theory, scattering amplitudes would be ill-defined in D = 4 dimensions. However, the need for introducing a regulator to deal with infrared divergences breaks conformal invariance and renders the amplitudes well defined. More specifically, we consider loop momenta in D = 4 − 2 dimensions while the momentum and polarization of the external particles remain in D = 4 dimensions. Infrared divergences thus arise as poles in with 1/ 2L at L loops. 6 The MHV all-loop 4-gluon amplitudes can be written as the tree amplitude, times a scalar factor expanded in a loop expansion via the rescaled 't Hooft coupling a, The BDS ansatz for the full rescaled 4-point amplitude M can be written as [19,20] with λ = g 2 N , f (λ) the cusp anomalous dimension, and g(λ) the collinear anomalous dimension. Expanding in loops L = powers of λ, it is easy to see that we have where f L (s/t, ) has a pole of order 1/ 2L . Alternatively, as found in [4], expanding in a L (instead of λ L ), one finds in terms of r = −s/t that Then, at one-loop, the celestial amplitudeÃ ) has an f (r,r), from (2.16), of a similar form with the tree level one:

JHEP01(2022)136
and then note that, because of the (−t) − factor, leading to w − in the exponent, we shift the argument of the function I(λ) in the tree celestial amplitude. We also have the extra factor F 1 (r, ), for a one-loop result Finally then, the celestial one-loop amplitude is written as an operator acting on the tree amplitude,Ã one−loop = aF 1 (r, )(P) Ã tree , (2.29) where the translation operator iŝ This translation operator is related to the translation operator P +,k that shifts the conformal dimension of a single gluon by one unit, ∆ k → ∆ k + 1, but involves the product, instead of the sum, over all k. This P +,k can be also used for the description of IR divergences, but in the case of an IR momentum cut-off for loop integrals, as in [5].
Then at L-loops, we have similarly leading to I(λ) → I(λ + 2iL ), and finallỹ Then, if we take the actual F L obtained from the BDS formula, we see that the total amplitude is an exponential factor, that now acts as an operator, acting on the tree amplitude.

Celestial gravity amplitudes and IR divergences
For n-graviton amplitudes in a gravitational theory, similarly we obtain IR divergences expressed as poles in for D = 4 − 2 dimensions.
It was shown in [21] that we can separate the amplitude in a hard part H n that is finite in , times a soft function S n , that contains all the IR divergences, where the soft function is where λ GR ≡ (κ/2) 2 (4πe −γ ) . Note that, unlike nonabelian gauge theories, like the case of N = 4 SYM treated above, there are no jet functions, since we have no collinear singularities

JHEP01(2022)136
in gravity, only in nonabelian gauge theory. As a result, the IR divergent factors in the exponent contain only single 1/ poles, not double 1/ 2 poles. However, note that, since for (massless) gravitons i s ij = 0, S n can be rewritten as a formal sum with 1/ 2 poles, more like the nonabelian gauge theory case, The object in the exponent in S n can be rewritten as in a similar manner with what happened in the case of N = 4 SYM in the previous section. We want to consider the effect of this factorization of IR divergences on the celestial graviton amplitudes.
The tree level 4-amplitude (which is the same in all gravity theories, including N = 8 supergravity) is written as a square of gauge theory amplitudes, where K is the same kinematical factor, involving helicities and momenta, as in the nonabelian gauge theory, defined further below in (4.3). But we saw that the tree-level MHV gluon amplitude had the form (2.11), with a constant B(s, t) = g 2 , so we infer that the graviton one also has the same form (2.11), with a constant B(s, t) = (λ GR u) 2 (the effective gravitational coupling, up to a number). Now we need to consider the factorization of gravity amplitudes (3.1), and see from it, analogously to the N = 4 SYM case reviewed in the previous section, what that implies for the loop amplitudes, and their IR divergences. We have to write a form for the function B(s, t) as in (2.31), in terms of a standard IR object, µ 2 −t times a function of r = −s/t, which will then get outside the integral over w = −t, and thus remain unchanged by the Mellin transform of the celestial amplitude.
Before proceeding, let us first express each object in the soft factorization (3.1) in terms of the loop expasion, i.e., n is the tree level amplitude and H L n are all infrared-safe hard factors. Thus, for instance, at one-loop we have In the absence of mass terms, there will be an effective dimensionless coupling that appears by absorbing the required powers of momenta (the only dimensionful quantities).

JHEP01(2022)136
In the case of gravity, that is λ GR s ij . The dimensional transmutation scale µ can only appear in the IR divergent part, and to fix overall dimensions. Then in the hard (IR finite) part of the amplitude H n , we can only have the physical dimension of the amplitude, given by the corresponding power of s, times a function of s/t. More precisely then, in this case, we must have (3.8) or, for n = 4, (3.9) and the expansion in λ GR t gives us the loop expansion. Thus, in particular, at one-loop, we have and A 4,tree contains the overall dimension as a power of s, times yet another function of the ratio s/t. In the one-loop amplitude, from (3.7) we see that, besides this hard factor there is an extra contribution coming from the soft factor S 4 , which at one-loop is given by just the exponent in (3.4), Then the B(s, t) coming from the one-loop amplitude is the constant from the tree times a function of s/t, which means that the function f one−loop (r,r, ) contains (λ GR (−µ 2 ))µ 2 , a function of the cross ratio,F 1 (r, ), and the integral Similarly, again, at L loops we find and soÃ For the first two terms in the expansion in of the L-loop amplitude, we do not actually need to know the full hard part of the gravity amplitude, only the one-loop amplitude of the corresponding theory, since we have in any gravity theory [16] which is an operator acting on the one-loop amplitude A (1) n , but depending on L instead of . Thus in this case, knowing F 1 (r, ) is enough to calculate the first two terms in the expansion in of the L-loop celestial amplitude.

JHEP01(2022)136
Finally then, the IR divergences exponentiate also in the celestial amplitude, now in terms of the operatorP −1+ , and we therefore find also the same soft/hard factorization for the celestial case.
We note that the IR divergences of graviton celestial amplitudes were considered also in [5], however in the case of IR momentum cut-off regularization for loop integrals, Λ IR . Then the factorization (3.1), with with s ij = 2p i · p j , leads to a similar one for celestial amplitudes, This result is consistent with ours, though our methods are different.

IR divergences of celestial gluon amplitudes: most-subleading-color
In this section we will consider the IR divergences of celestial gluon amplitudes in a general gauge theory, and we will note that, as in the case of usual amplitudes, the most-subleadingcolor term has properties that are similar to the ones of gravity, making it easier to write a factorized expression.

Set-up: color and N expansion, factorization
The full amplitude is expanded in a trace basis {T λ } of single and multiple-traces of generators T a in the fundamental representation where the coefficients A λ are the color-ordered amplitudes. The A λ elements are organized into a ket |A , according to the standard notation in [22,23]. For 4-gluon amplitudes (n = 4), λ = 1, . . . , 9 (six single traces and three double traces), however, since only three out of the six single traces are independent, one can choose a basis using only three of them. A useful such basis is given by

JHEP01(2022)136
In this basis, the 4-gluon tree level amplitude is written as where K is a factor depending on the helicities and momenta of the external gluons and can be found in equation (7.4.42) in [24]. At loop level, for the gauge group SU(N ) (the relevant case for phenomenology), the amplitudes can be expanded in both L and 1/N , with coefficients A (L,k) , is the 't Hooft coupling and µ is the renormalization scale. Here k ≤ L and, anticipating our discussion of the following subsection, notice that for k = L there is no overall power of N multiplying A (L,k) . These (infinite) subset of amplitudes are the ones known as most-subleading color gluon amplitudes. At one-loop, we then only have A (1,0) and A (1,1) . However, unlike the tree level amplitude, these are theory-dependent functions. For the case of N = 4 SYM, we have  where the massless scalar box I 4 (s, t) has an exact expression to all orders in , given by Note that I 4 (s, t) can be written as t −2− times a function that only depends on the cross-ratio r = −s/t. The explicit expansion of I

(s, t) in powers of reads
From the first line above we also note that, when Mellin transforming these one-loop expressions, the first two terms will simply amount to the change I(λ) → I(λ + i ) (with I(λ) defined in (2.21)) plus constant terms since they are either pure numbers or they depend on ratios of Mandelstam invariants only.

JHEP01(2022)136
For N = 1, 2 SYM theories, the full one-loop amplitude is (along with its A (1,0) and A (1,1) components), to the best of our knowledge, not yet available in the literature. However, the general form of the amplitude being of the type t −2− times a function of r = −s/t is still valid here. Moreover, it is not only true at the one-loop level but at all-loop orders.
In order to see this, let us consider the following argument, which is actually even simpler than the one used for the case of gravity in the previous section. The theories we are considering here have no scales at all since the couplings are dimensionless. Therefore, in general, we must have that the IR-finite part (i.e., the hard factor) consists in an overall power of s ij , (maybe times a power of s ij /µ 2 for the anomalous dimension, if it is nonzero, although this is usually taken care of by renormalization and by the factoring of IR divergences), times another function depending purely on the ratios s ij /s kl and of the coupling, that is Note that the full loop expansion is obtained from the factor f (a(µ 2 ), r) which is just a function of the ratio r = −s/t.
In the above relations we have used a notation for the IR finite part that previews an important fact about general color-ordered gluon amplitudes: they can be factorized into the product of three objects: a jet function J, a soft factor S, and a hard IR-safe factor H. The jet factor J is just a scalar function (a number) while the soft function S is a matrix-valued object (defined by the soft anomalous dimension matrix Γ below). Both jet and soft factors are IR divergent, thus describing the large distance behaviour, while the hard function H, being IR-finite (as → 0), describes short distance physics. In the ket notation for the amplitudes, the factorization reads [25,26] A s ij µ 2 , a(µ 2 ), = J a(µ 2 ), S s ij µ 2 , a(µ 2 ), H s ij µ 2 , a(µ 2 ), . (4.12) The soft function, for any gauge theory, is generically where the central object Γ is known as the soft anomalous dimension matrix. 7 The loop expansion of Γ, 14)

JHEP01(2022)136
is fully known (for massless theories) up to three loops [28], while four-loop corrections have recently appeared in the literature [29]. The one-loop matrix is where T i is the color charge operator for the i-th gluon corresponding to the SU(N ) generators in the adjoint representation. Here T i · T j ≡ T a i T a j , T 2 i = C i , with C i being the quadratic Casimir where, in the adjoint representation, C i = N . The matrix element of T i in between states characterized by gluon adjoint indices b i , c i is Another way of characterizing these operators is by writing how they act upon the generators which defines the action on the color basis. Given this action of T i , in terms of the color-ordered expansion (4.1), the L-loop anomalous dimension matrix Γ (L) acts on a given element T λ of the trace basis (4.1) to yield a linear combination and it is the matrix Γ (L) κλ that then acts on the hard ket |H . The fact that the two-loop contribution to the anomalous dimension matrix is proportional to the one obtained at one-loop, together with arguments based on collinear limits and symmetry considerations regarding rescalings of the momenta of the external states, led to a fascinating conjecture coined the dipole conjecture, valid for any gauge theory, including QCD. This stated that all higher loop corrections in the soft anomalous dimension matrix Γ were proportional to the one-loop matrix Γ (1) [30][31][32]. However, as argued in [33] and later shown by an explicit 3-loop calculation [28], the dipole conjecture fails starting at three-loop order. 8 Within the dipole approximation, it is easy to see from (4.15) that [Γ(µ 1 ), Γ(µ 2 )] = 0, for any two different scales µ 1 = µ 2 . This makes the path ordering operator in (4.13) immaterial allowing us to directly compute the integral, yielding Note however that in this case, we still need to calculate all of the Γ (L) in order to have the full IR behaviour in the soft function S. We would indeed have an exponentiation due to the resummations, but in the exponent we still have to sum over loops. Moreover, one would also need to calculate the jet function J.

JHEP01(2022)136
In the rest of this section we will focus on the effect of the dipole structure and in section 6 we will explore corrections from non-dipole terms.
Under the dipole conjecture, and considering only the one-loop anomalous dimension contribution in the exponent in (4.19) to all orders in N , we obtain (4.20) with the understanding that, in the second equality, only the divergent and finite terms in the expansion are to be kept. Note that, due to color conservation n i=1 T i = 0, n j=1 i<j the second term in the exponent of (4.20) commutes with the first. Therefore, only the action of the first term upon the hard amplitude |H is non-trivial. In the four gluon case, (4.22) from where we see that the square bracket is a function of r = −s/t only, while the prefactor is the standard one that already appeared in N = 4 SYM and in gravity.
So in this case, yet again, when calculating the celestial amplitude, f one−loop (r,r, ) for the first term in the exponent of S contains aµ 2 , a function F 1 (r, ), that gets outside of the integral over w = −t, and the same I(λ + 2i ), since we have (4.23) and so involves the same translation operatorP as before. Its action is now not on the tree amplitude, but rather on the hard amplitude |H , which is a function of r (times the trivial dependence on t given by its dimension). The second term is just a constant. Then, at L loops, we obtain yielding I(λ) → I(λ + 2iL ) in the leading term, and so Finally, note that

JHEP01(2022)136
and α, δ are written in terms of S, T , U themselves, or their sums, where whereas β, γ are written in terms of differences of the same, so are functions of ratios of s, t, u. This implies that the leading (in N ) part of Γ (1) depends independently on s, t, u, whereas the subleading (in 1/N : the part with β and γ) part, Γ sub , depends only on the ratios.

IR divergences for the most-subleading-color amplitudes
We have seen that we can write a formula for the IR divergence of the celestial gluon amplitudes using the dipole conjecture, but considering only the one-loop term in the exponent.
However, and this is the crucial point, if we focus on the (still infinite set of) mostsubleading-color amplitudes at all loops, |A (L,L) , which contain only terms of order N 0 (the others contain positive powers of N ) we can completely determine the IR divergence [16]. Indeed, while in general Γ  where and For the celestial amplitude, we note that these functions are all functions of r = −s/t, so in the B(s, t) integral over w = −t, they become mere constants that we can take out of the integral. Therefore the IR divergence is a simple constant factor leaving only the integral of the tree amplitude, yielding I(λ). Therefore, the formula (4.28) continues to hold for the celestial amplitude.
Moreover, for the first two leading divergencies for each most-subleading amplitude A (L,L) , for L = 2l (even) or L = 2l + 1 (odd), we have the relations (first conjectured for JHEP01(2022)136 N = 4 SYM in [11] and then proven in [12], if and when the dipole conjecture holds, for a general gauge theory in [16]): Since, as we saw, X, Y, Z are functions of only ratios of Mandelstam variables, and get outside the w = −t integration in the celestial amplitude, equations (4.29), (4.33) and (4.32) for 4-graviton amplitudes are also valid in celestial amplitudes.
Here we should note that the behaviour of (4.33) and (4.32) for the most-subleading amplitude is similar to the gravity relation (3.15). In fact, the relation between the two can be made more precise, as was done in [16] for the usual amplitudes: (s log s + t log t + u log u) 2 2(log 2 (t/u) + log 2 (u/s) + log 2 (s/t)) 1 /u = −4iK/(stu) is symmetric in s, t, u). Now again X, Y, Z are only functions of the ratios of s, t, u, however s, t, u appear explicitly in the relations. We could note that, considering the form on the first line(s) of the equation(s), we can take outside an u L factor in both cases (L = 2l + 1 and L = 2l), obtaining the modified coupling factor (λt) L in front, and then only ratios of s, t, u, and the ratio of the amplitude to the tree ampliture, A (L,L) /A (0) 1 , just like on the left-hand side. Therefore, when going to the celestial amplitudes, the relation is the same, except for the extra w L = (−t) L in the integral over w giving I(λ), shifting its argument, which, as we saw, means an extra translation operatorP for each −t factor. So the (λt) L transforms into (λP) L , but otherwise the relation is unchanged.
For loop number L = 0, 1, 2, the relation between the most-subleading 4-point amplitude divided by the tree amplitude and gravity 4-point amplitudes divided by their tree amplitudes is actually exact, as shown in [11] for N = 4 SYM vs. N = 8 supergravity and JHEP01(2022)136 in [15] at any N : SYM,N (s, t) + cyclic perms of s, t, u , (4.36) where λ SYM = g 2 N and λ SG = (κ/2) 2 . Again the same observation holds: we can take a factor of w L = (−t) L common, so as to have (λ SG w) L and the rest just ratios of s, t, u and M SYM , and then we just shift the integration over w, so shift the argument of I(λ), or act with the translation operatorP L , giving (λP) L .
For the 5-point functions, there are also some relations. At one-loop, we have the exact relation (found for N = 4 SYM in [34] and proven for general N in [15]) between ratios of loop amplitudes and tree amplitudes, where β 12345 are one-loop numerators defined in [35]. We also have a relation between tree amplitudes involving the same numerators, found in [13] for N = 4 and used at general N in [15], Then, in the case of 5-point functions, going to the celestial amplitudes would result in more complicated relations between the two sides, because of the numerators β, so we will not describe them.

The soft factor and the Knizhnik-Zamolodchikov equation
In this section we consider an important relation of the celestial soft factor to the Knizhnik-Zamolodchikov equation. We start by reviewing the celestial analog of the soft factorization in nonabelian gauge theory, following [6] and using their notation. In dimensional regularization, the all-loop n-gluon amplitude in momentum space takes the factorized form (the jet function is now absorbed into the soft factor) with the soft factor (previously split into the -matrix-soft function S and the -scalar-jet function J) where α = 2e − γ E g 2 /(4π) 2− is the YM coupling and the (total, i.e. all-loop) soft anomalous dimension matrix Γ n is given by where γ(α) is the cusp anomalous dimension.

JHEP01(2022)136
The first term is the dipole term, described in detail in the previous section (see (4.15 for its one-loop part, the further loops correct just its coefficient, giving the cusp anomalous dimension γ(α) in front), and which accounts for pairwise interactions only. The second term is the sum over collinear anomalous dimensions of each of the external gluons (jet function contributions in the notation from the previous section). The non-dipole contributions ∆ n and ∆Γ n make their appearance starting at three [28] and four loops [29] respectively, and were ignored until now.
Using the momentum parametrization (2.1), the Mandelstam variables read s ij = η i η j ω i ω j |z ij | 2 . From here and the color conservation condition we see that Γ n reads (we have denotedγ(α) = γ(α)/C A , where C A is the Casimir in the adjoint representation, equal to N here, and C k refers to the Casimir of gluon k, obtained This separation of Γ n into a part logarithmically dependent on the energies ω k , plus the same function evaluated on the celestial data |z ij |, allows us to write whereK(α) = 1 2 µ 2 0 dλ 2 λ 2γ (α(λ 2 )), and therefore, from (2.12), we immediately see that the full n-gluon celestial amplitude becomes [6] i.e., the momentum space soft factorization persists after Mellin transforming the amplitudes into celestial correlators. Moreover, the full effect, at all loops and all N contributions (to compare with the previous section, there we considered only the most subleading component, as well as the one loop in all N components) of the infrared divergences onto the hard part of the amplitude is a renormalization of the conformal dimensions from their bare values ∆ i to the renormalized ones ∆ i − 1 2 C iK . When → 0,K diverges, giving an infinite shift, like it was explicitly found in the case of N = 4 SYM.
Here we will only concentrate on the contribution from the dipole term in Γ n , leaving the effects of the contributions ∆ n and ∆Γ n for future work, as we did in the previous section. In this case, it is easy to see that, for two arbitrary scales λ 1 = λ 2 , the commutator [Γ n (λ 2 1 ), Γ n (λ 2 2 )] vanishes, thus making the path-ordering operator P in (5.2) immaterial, as explained in the previous section. This allows us to write Z as

JHEP01(2022)136
where we have used color conservation n i=1 T i = 0 and the definitions We have also assumed that all external states are of the same species (e.g., gluons), that is, Writing then the infrared factor as the scalar jet function J n and the soft function e A , and taking the derivative of Z with respect to any of the positions z k on the celestial sphere, we have where ∂ k = ∂/∂z k . The first term is Thus, if the commutator expansion above can be thought of as a perturbative expansion in a small parameter with ∂ k A being its leading term, the soft celestial factor Z α, {µ 2 |z ij | 2 } indeed satisfies the Knizhnik-Zamolodchikov (KZ) equation in that limit, i.e., From here, one would identify the infrared divergent factor K( ) and N with the level κ and the Coxeter number g as For the particular case of the most-subleading four-gluon amplitude we will see that, at leading order in the Regge limit, one indeed obtains the Knizhnik-Zamolodchikov equation. The most-subleading amplitude has the form where the matrix M, for the 4-gluon case at most-subleading order, is given by (5.17) Here we have changed our momentum conventions, in order to perform the Regge limit in the physical region (r > 1). Notice that, unlike the general gluon amplitude, the dependence on the z k coordinates is through the cross-ratio r only. Thus, for large r, we have, Therefore, at leading order in the large r expansion all commutator terms in (5.11) vanish, yielding the KZ equation in the Regge limit r → ∞.
To conclude, note that the dipole contribution A (for finite N ) can be written as (1) . (5.19) Using this and the fact that A (0) is a diagonal matrix (see, e.g., the four-point case in (4.26)) the commutator expansion in (5.11) organizes into a 1/N expansion of the form thus, in the large N limit, we again see the appearance of the KZ equation. However, since only the diagonal part of appears at leading order, the CFT describing the soft celestial part seems to be a theory of non-interacting colored scalar fields, as previously hinted in [6,8].
We would like to close this section by mentioning that the appearance of the KZ equation also emerges in other kinematical limits, and for the general gluon amplitude, i.e., not only for the most-subleading-color amplitude (within the dipole aproximation). For instance, in collinear limits such as z 12 = z 34 = δ → 0, one also obtains that the commutators in (5.11) are subleading with respect to the first term in the expansion in powers of δ, although again the nonabelian structure seems to be too simple, similar to the large N case. It would certainly be interesting to explore these limits in more general situations, such as in higher-point amplitudes and in other multi-collinear regimes.
Finally, note that in [36], the authors obtained the Regge limit of the soft factor Z of the usual (non-celestial) amplitudes (within the dipole approximation). 9 Therefore, Z must reduce to the soft factor in (5.15) for the most-subleading-color case, which is indeed the case.

JHEP01(2022)136
Being so, we can import the calculation from there, and find An interesting observation is that, for the most-subleading color part of this first term in (6.3), a = d = 0, so which is the same block form that we have for Γ (1) most−sub . However, the form of the matrices b and c are different from the ones in Γ (1) most−sub , and the two matrices do not commute (they just act in the same way, switching between single traces and double traces in the trace basis). This implies, for instance, that these 3-loop corrections to the KZ equation are rather non-trivial since they would involve new extra terms in the commutator expansion in (5.11).
For our purposes, when computing the contributions of these corrections to celestial amplitudes, the important observation is that ∆ 4 is only a function of r = −s/t (and more generally, ∆ n is a function conformal invariant cross ratios only), so its contribution to the celestial amplitude is trivially found since it just contributes to the F L (r, ), that is, it becomes a simple constant that comes out of the integral over w = −t in (2.16).

JHEP01(2022)136
Note also that, as explained in [28], in the Regge limit r = −s/t → ∞, ∆ 4 becomes a constant matrix (although still not one that commutes with the dipole terms).

Conclusions
In this paper we have considered the IR divergences of celestial amplitudes in gravity and in gauge theories, in particular the most-subleading-color part of the latter. We have found that the simplifications that occur in the IR structure of the most-subleading gluon amplitudes, responsible for their close connection with gravity, also translate into a simple form of the celestial correlators yielding a simple exponentiation and factorization. The fact that the soft factor Z in the most-subleading-color case is a function of conformal cross-ratios only plays a crucial role in the simple form of the celestial soft factorization.
Gravity amplitudes were found to have an exponential IR factor acting on the hard celestial part. If the dipole conjecture is assumed to be true, (such as in the case of N = 4 SYM where it is true), we were able to compute exactly the IR divergent factor acting on the hard part for the most-subleading celestial gluon amplitudes, as well as to reduce the two most divergent in terms to the calculation of the one-loop amplitude A (1,1) , just like in the gravity case.
Also, under the dipole approximation, for the all-order in N soft factor Z in terms of the cusp anomalous dimension γ(α), we have found a generalization of the Knizhnik-Zamolodchikov (KZ) equation, reducing to the KZ equation in the case of the mostsubleading celestial gluon amplitude in the Regge limit, as well as for the full amplitude at large N . Finally, we have also considered the 3-loop corrections due to non-dipole terms, and how they impact the IR divergent terms.
To conclude, we would like to close with a few observations and further directions: • The KZ equation in the WZW model is a consistency relation as a consequence of the presence of a null state in the theory. Since the KZ equation is not satisfied here in general, suggesting the absence of null states, but it does arise in the large cross-ratio limit r → ∞, this seems to imply that a null state does emerge in this particular limit. Also, as mentioned at the end of section 4, there could be other kinematical regions where the KZ equation is satisfied and it would definitely be interesting to explore this further.
• Possible further directions include: (a) compute the subleading contributions to the Regge limit; (b) include ∆Γ n corrections to the dipole term appearing at four loops; (c) consider the Regge limit of the full gluon amplitude (not just most-subleading part), using the results of [36]; (d) search for a non-trivial simplification also in the large N limit; (e) consider the effect of the all-loop conjectured connection between N = 4 SYM and N = 8 supergravity in the Regge limit in the recent paper by S. Naculich [37].