Descendants in celestial CFT and emergent multi-collinear factorization

Multi-collinear factorization limits provide a window to study how locality and unitarity of scattering amplitudes can emerge dynamically from celestial CFT, the conjectured holographic dual to gauge and gravitational theories in flat space. To this end, we first use asymptotic symmetries to commence a systematic study of conformal and Kac-Moody descendants in the OPE of celestial gluons. Recursive application of these OPEs then equips us with a novel holographic method of computing the multi-collinear limits of gluon amplitudes. We perform this computation for some of the simplest helicity assignments of the collinear particles. The prediction from the OPE matches with Mellin transforms of the expressions in the literature to all orders in conformal descendants. In a similar vein, we conclude by studying multi-collinear limits of graviton amplitudes in the leading approximation of sequential double-collinear limits, again finding a consistency check against the leading order OPE of celestial gravitons.


Introduction
Recent decades have seen many discoveries of alternative mathematical structures from which the standard principles of perturbative QFT emerge as derived consequences. One of the primary motivations of such investigations has been to find a holographic description for scattering amplitudes in flat space, akin to the highly successful AdS/CFT paradigm. In this context, celestial conformal field theory (CCFT) is a recent proposal that claims to identify Yang-Mills and gravitational amplitudes in R 1,3 with correlators of a putative 2d CFT living on the celestial sphere at null infinity. And even though no explicit candidate or stringy construction for such a holographic dual has been found yet, great progress is being made in understanding the abstract structures and symmetries that such a CFT could possess. Some of the notable advances include the main work on celestial amplitudes [1][2][3][4][5][6][7][8][9][10][11], on asymptotic symmetries and soft theorems [12][13][14][15][16][17][18][19][20][21][22][23][24][25][26][27], and on the CCFT operator algebra [28][29][30][31].
Due to the absence of an actual candidate CFT, most of the work has been kinematical and one-sided: trying to understand the CFT side by studying properties of the amplitudes. This begs the question: how do we begin to test this holographic proposal? An interesting direction was taken in [29], which provided a purely holographic derivation of the CCFT operator product expansions (OPE) via imposing asymptotic symmetry constraints. In turn, this gave a new holographic foundation for the universality of the well-known collinear limits of gluon and graviton amplitudes [32][33][34]. The work of [30,31] took this further by showing that even subleading terms in the collinear expansions of low-multiplicity graviton amplitudes can be ascribed to BMS descendants in the gravitational CCFT. Coupled with the idea that a CFT is in principle completely determined from its CFT data, i.e., its operator content and OPE algebra, such computations are enough to allow us to come up with some simple tests of the duality.
One of the hallmarks of scattering amplitudes is the structure of their factorization poles and residues. These are completely fixed by the principles of locality and unitarity. An important test of any holographic dual would then be to discover them as emergent properties of the corresponding CFT correlators. From the viewpoint of the OPE, the most natural object to study in this regard are multi-collinear limits of the amplitudes. These are maximally singular limits that recursively probe all possible factorization poles and residues (see [35] and references therein). In this work, we show that these can indeed be holographically determined by the symmetries and OPE of the dual CCFTs. This provides an example of a calculation that utilizes the CCFT to essentially "bootstrap" the physics of amplitudes.
Such a reconstruction of the bulk physics requires us to understand the contributions of descendants to the celestial OPE, and we will mostly focus on the gluon OPE for sake of simplicity. After a review of some standard material in §2 and §3.1, we begin with this task in §3.2. Global supertranslation symmetry is used to fix the OPE coefficients of all the (global) conformal descendants of celestial gluon operators. For completeness, in §3. 3 we also compute examples of Kac-Moody descendants contributing to the OPE of two positive helicity gluons. This is done by imposing Poincaré as well as Kac-Moody invariance. In fact, naively the symmetries overdetermine the OPE coefficients, but the results are beautifully mutually consistent. In appendix A, we directly verify that these descendants are exchanged in the 4-gluon amplitude with precisely the predicted OPE coefficients.
In §4, we look at multi-gluon collinear limits. In the language of conformal correlators, this corresponds to bringing multiple operators close together. To get a feel for the idea, consider a correlator of N operators in a CFT 2 , with (z i ,z i ) denoting complex coordinates. The various operators have operator products taking the generic form, with C ijk ≡ C ijk (z ij ,z ij , ∂ j ,∂ j ) denoting Wilson coefficients that depend on the CFT data and z ij ≡ z i − z j . Now, for operators in a given ordering, say |z 12 | < |z 23 | < · · · < |z n−1,n | (where n ≤ N ), we can replace the n-fold product n i=1 O i (z i ,z i ) by using the OPE to perform n sequential Wick contractions: Then the product of Wilson coefficients gives the equivalent of a multi-collinear splitting function (or more appropriately "splitting operator") for CFT correlators. It is automatically universal since it does not depend on the other N − n operators inserted in the correlator.
Thus, having determined the CCFT Wilson coefficients to the desired accuracy, we can approximate celestial amplitudes with such recursive celestial OPEs in the limit of small z ij . As the main utility of these coefficients, we will find the leading multi-collinear factorization behavior of the usual momentum space amplitudes without any input whatsoever from Feynman rules or the usual techniques of 4d QFT. Finally, the existence of a CCFT interpretation can guarantee the universality of these limits. Some integration techniques relevant to these computations are described in appendix B.
In §6, we attempt a similar calculation for graviton amplitudes. Due to a lack of literature to compare with on the gravitational side, we will only be able to outline a leading order computation for the simplest multi-graviton collinear limits. This nevertheless provides a nice consistency check of the formalism and a motivation for further work.

Background
In this section, we collect some standard conventions about celestial amplitudes and results for celestial gluon and graviton operator product expansions. Then we review the basics of multi-collinear limits that will come to use later, noting relevant results from the literature.

Celestial amplitudes and OPE
The null 4-momentum k αα of a typical massless particle can be decomposed as Here, the sign s = ±1 denotes whether the particle is outgoing or incoming, while ω ∈ R + denotes its energy. The remaining null vector q αα stands for the embedding of the celestial sphere CS 2 as the projective light cone of any point in flat space, with z,z giving complex coordinates on the sphere. By convention, the corresponding spinor-helicity variables are taken to be With Lorentzian signature (− + + +), one chooses the standard reality condition on these:z = z * (complex conjugation). We will stick to this, except for using split signature (− − + +) in appendix A.
"Celestial amplitudes" are the scattering amplitudes of conformal primary wavepackets of gluons and gravitons. In short, they can be defined in terms of a change of basis implemented by Mellin transforms: Here i = 1, . . . , n are particle labels, with k i = s i ω i q i , q i ≡ q(z i ,z i ) as described above. A n denotes the usual momentum space amplitude. For gluons, it will be further augmented with color indices. i = ±1 or ±2 are the gluon/graviton helicities. Under a Möbius transformation of CS 2 coordinates, the celestial amplitudes A n transform conformally covariantly with weights [2,3] in the (z i ,z i ). Consequently, these are conjectured to be the correlators of certain conformal primary operators in a 2d CFT living on CS 2 , called a celestial CFT. Such conformal primaries dual to gluons and gravitons are referred to as celestial gluon/graviton operators. Celestial gluons are denoted by O a, s ∆ (z,z), where a is an adjoint index. Celestial gravitons are commonly denoted by G ,s ∆ (z,z). For most of what follows, we will focus only on outgoing particles for which s = +1, so we drop this superscript to avoid cluttering notation.
Collinear singularities z ij → 0 in momentum space amplitudes are interpreted as OPE singularities of the correlators obtained from the Mellin transform. This has allowed a determination of the leading celestial OPE of gluons and gravitons [28,29]. For instance, the OPEs of two outgoing gluons read

Multi-collinear factorization
Gluon and graviton amplitudes in flat space are completely characterized by on-shell methods like BCFW recursion and MHV diagrams [36][37][38]. Such relations follow from their factorization poles and residues as a consequence of locality and unitarity. However, unlike collinear limits, general factorization limits don't seem to have an obvious analog for celestial amplitudes. One possibility to make a connection between the two formalisms is thus to study more involved collinear singularities, namely multi-collinear limits. They correspond to simultaneously taking all possible factorization limits involving a certain subset of the scattering particles. Such limits have been well-studied in the literature on QCD, with important progress originating from the usage of CSW rules and MHV diagrams [35]. For gravity, there has been no such progress beyond the 2-graviton collinear limit, more or less due to a lack of strong theoretical foundations for similar recursive methods. 1 The kinematic configuration probing a multi-collinear singularity corresponds to the null momenta k 1 , . . . , k n of a subset of the particles becoming collinear. For simplicity, we take all of these to be outgoing. Then all the propagators of the form (k i 1 + · · · + k ir ) −2 , i 1 , . . . , i r ∈ {1, . . . , n}, diverge. This leads to a maximally singular sub-amplitude to bubble off, yielding a universal factor called a splitting function. From a 4d perspective, its universality in the gluon case again follows from MHV diagrams. From the holographic viewpoint, our claim is that its universality is a consequence of the celestial OPE -an argument which might also extend to gravity.
To make the limits precise and set up some notation, note that we can always express the sum of multiple momenta in terms of two auxiliary null momenta, k 1 + k 2 + · · · + k n = p + n , (2.9) where for instance we can somewhat canonically choose n to be the null generator of I + . It follows that since ω i = n · k i , this choice results in where ω p = n·p is the energy of p. Define the longitudinal-momentum fractions ξ i := ω i /ω p . The collinear regime is defined by along with taking to be infinitesimal. In this regime, an N -point momentum space amplitude (with n ≤ N ) factorizes as where the superscripts are particle helicities. The universal splitting functions split (· · · ) are neatly organized by the number of negative helicity gluons participating in the collinearity.
If 1, 2, . . . , k are negative helicity among the n collinear particles, then the corresponding splitting function is denoted by (2.13) These will take center-stage in the latter half of this work. For gluons, we will content ourselves with considering the k = 0, 1 cases. In this case, the collinear gluons must be adjacent for the splitting function to be maximally singular. For these, the results for the splitting functions come in fairly compact expressions found in [35]. Using the convention (2.2), we can write them in variables adapted to the celestial sphere. The simplest of these occur in the case when all the collinear gluons have positive helicity (in this case we denote them by split (n) without any arguments), , (2.14) having stripped off color factors (these will be reinstated for comparison with OPE later). When the first gluon is negative helicity, one finds along with the relatively much more interesting expression, split (n) where s 1j is the generalized Mandelstam variable, At the level of the first three among these, one only observes two-particle factorization poles.
In the language of MHV diagrams, this is a consequence of the fact that only MHV subamplitudes happen to blow up for these configurations. The collinearity 1 − 2 + · · · n + → p + is the first case where NMHV sub-amplitudes can blow up. It will be much more novel for celestial CFT to make contact with multi-particle factorization poles of the sort showing up in (2.16), even if only in leading order approximations in some of the variables. This will be our goal in §4. For gravitons, as mentioned above, there is a distinct lack of data for multi-collinear limits beyond n = 2. The double-collinear limits that we need are given by [33] split (2) The graviton collinear limits are not singular in general and will even depend on the order in which the momenta are made collinear. So we will have to restrict our analysis to the simplest case: using the double-collinear limit to sequentially compute leading order approximations to the multi-collinear splitting functions.

Symmetry algebra
In a conformal field theory with an extended symmetry, the states and their corresponding local operators arrange themselves in representations of the symmetry algebra. In the 2d CCFT dual to 4d Yang-Mills, the conjectured symmetry algebra is Poincaré plus a holomorphic Kac-Moody symmetry [40][41][42]. 2,3 The representation multiplets are then organized into primaries of the Kac-Moody symmetry and their global conformal descendants and Kac-Moody descendants.
The 4d Lorentz group acts as the global conformal group of CS 2 . Its generators can be denoted by the standard combinations {L 0 ,L 0 , L ±1 ,L ±1 } of SL(2, C) dilatations and rotations. The generators of global supertranslations P a,b are identified with translation generators (momenta) in the bulk, with the following convenient arrangement: Their algebra is given by [17,24,30] [L m , Celestial gluons are primaries of the full Poincaré group, with transformation laws: Notice that the action of P a,b induces a flow, (h,h) → (h + 1 2 ,h + 1 2 ), in the conformal dimensions.
As for the Kac-Moody generators, the holomorphic current is identified with the conformally soft limit of the outgoing positive helicity celestial gluon [16,19,28], Using the OPEs (2.5) and (2.6), one can show that celestial gluons also transform as Kac-Moody primaries in the adjoint representation, Taking = +1 and ∆ → 1 again in this OPE gives the usual OPE bewteen Kac-Moody currents at level 0. Expanding the current in its holomorphic modes, one can straightforwardly find their action on O a ∆ from (3.7), Including the Kac-Moody modes, we obtain the following extended algebra: The last of these can be justified by taking a conformally soft limit of (3.5), or equivalently by showing that the OPE of [P a,b , j a (z)] with an arbitrary celestial gluon is non-singular. These form the leading 4 global symmetry algebra of the Yang-Mills CCFT. A typical descendant occurring as a subleading term in the OPEs (2.5), (2.6) is of the form, with each k a ≥ 1. Adding these to the OPEs with unknown OPE coefficients, one can apply various symmetry generators to generate constraints on the coefficients. However, the constraints coming from conformal and Kac-Moody actions become increasingly cumbersome very quickly. But quite luckily, since P a,b commutes with the Kac-Moody generators, its action does not mix descendants of different Kac-Moody weights. In §3.2, this allows us to use just translation invariance to fix the OPE coefficients accompanying all the purely conformal descendants (L −1 ) m (L −1 ) n O c ∆ without having to worry about the Kac-Moody descendants. For completeness, we also give an example of how OPE coefficients of some j a −1 -descendants can be determined from this data in §3. 3.
Note however that this idea hinges on the assumption that there are no further global symmetries of the Yang-Mills CCFT whose descendants might mix with these purely conformal descendants under translations. Thus, our calculations in §4 and appendix A also bolster our confidence that these are all the kinds of descendants that occur in the gluon OPE. However, more work along the lines of [31] might be needed to confirm this.

Conformal descendants
We will use the action of P 0,−1 and P −1,0 to determine the coefficients of conformal descendants in celestial gluon OPEs. To do this, we need the following easily derived relations, and evaluated at the origin for simplicity. Let us denote the contribution of a celestial gluon of helicity to the OPE of two gluons of helicity 1 , 2 by (3.14) where (h,h) are the conformal weights of O c ∆ 1 +∆ 2 −1 , and C 1 2 ∆ 1 ,∆ 2 encodes contributions from this primary and its conformal descendants: Here, we are viewing L −1 ,L −1 respectively as the holomorphic and antiholomorphic derivatives ∂ 2 ,∂ 2 when taken against the primary O c ∆ 1 +∆ 2 −1 (z 2 ,z 2 ). To find recursion relations on the coefficients C 1 2 m,n (∆ 1 , ∆ 2 ), we set z 2 = 0 =z 2 and act with P 0,−1 and P −1,0 on (3.14). Applying (3.5), (3.12) and (3.13), this process generates From the leading OPE (2.5) and (2.6), we have the following boundary conditions for these recursion relations: Thus, solving (3.16), we readily discover the entire series of OPE coefficients, where (a) q := Γ(a + q)/Γ(a) are Pochhammer symbols. Hence, we observe that translation symmetry is a very powerful constraint on the structure of the CCFT which would be absent from a garden-variety CFT. Performing such an all order calculation using Virasoro and Kac-Moody symmetry constraints is almost inconceivable, and it is almost always more useful to work out 4-point conformal blocks rather than the descendants' OPE coefficients. However, knowing all order contributions to the OPE as we do here will help us make contact with interesting universal statements about scattering amplitudes at arbitrary multiplicity.

Kac-Moody descendants
As an example, we analyze the first Kac-Moody descendants contributing to the like-helicity OPE (2.5). In particular, we need to justify that the OPE data generated above using translation invariance is consistent with the other symmetries in the problem.
To begin with, let us take the following ansatz containing j a −1 -descendants: . Now, we observe that c 1 can be read off from (3.18) evaluated for m = 1, n = 0, To fix c abcd 2 , we set z 2 = 0 =z 2 and act with j e 1 on this OPE. Using (3.9) and (3.10), this yields the relation, To solve for the remaining coefficient, we further guess an ansatz of the form One can in principle add other group-invariant "tensor structures" to this ansatz, like higher degree polynomials in structure constants, but we won't need them at this level. Courtesy of the Jacobi identity, this already satisfies (3.23) for the values, Hence, we have So, we find not one but two descendants contributing at this level: j a −1 O b ∆ 1 +∆ 2 −1 as well as its permuted partner j b −1 O a ∆ 1 +∆ 2 −1 . Next, we work out the constraint coming from an application of L 1 . The resulting relation is Substituting (3.22), (3.24) into this, we find the excess condition, which a priori overdetermines the system of equations. But this is already beautifully satisfied by our solution (3.26). This demonstrates, at least by way of an example, that the enormous amount of symmetry in a CCFT can indeed allow for non-trivial CFT data consistent with all of it.
For the details of how these descendants can be extracted from the OPE limit of an actual 4-gluon amplitude, the reader is directed to appendix A.

Multi-gluon collinear limits
The OPE of two celestial gluons was derived in [28] by Mellin transforming the doublecollinear splitting functions. The goal of this section is to generalize this computation to obtain multi-gluon OPE from Mellin transforms of multi-collinear limits. But the former can also be computed holographically by recursively applying the 2-point OPE. Most importantly, in §4.2.1, we match the contributions of conformal descendants derived in §3.2 to all orders with the triple-collinear limit. This provides a mechanism for factorization poles and residues of 4d amplitudes to emerge from the CCFT, generating the footprints of locality and unitarity.
A celestial amplitude of N gluons, the first n of which are outgoing and will be taken collinear, can be written as a CCFT correlation function, We can use the OPE to expand this around the multi-collinear regime. Fixing an ordering of points |z 12 | < |z 23 | < · · · < |z n−1,n |, the OPE of the collinear gluons can be accessed by sequentially applying the 2-gluon OPE. It takes the general form, where the quantity "ope" (in general a differential operator) is the celestial analog of the color-stripped splitting function in (2.12) and we have suppressed its dependence on the ∆ i 's. Inserting this in (4.1) leads to universal asymptotics, But (4.3) can also be obtained from Mellin transforming (2.12) (after reinstating color factors). Hence, to leading order in the collinear kinematics, holographic duality will relate the multi-gluon OPE with the splitting functions via Mellin transform, where both sides are to be viewed as "integration kernels" acting on the remaining nonsingular momentum space amplitude. Matching these will be our primary "test" of celestial holography. For brevity, we will use the same notation ope (n) (1 − 2 − . . . k − ) for the ope factor that we introduced for the splitting functions in (2.13). Recall that this corresponds to the case when gluons 1, 2, . . . , k among the collinear ones have negative helicity. Our focus will be on the k = 0, 1 cases, with the notation being just ope (n) when k = 0.
Next we perform a direct calculation in the bulk. A celestial amplitude involving these gluons factorizes as with ω p = n i=1 ω i and q p ∼ q n , and the relevant splitting function given in (2.14). The Mellin integrals over ω 1 , . . . , ω n can be easily done via a change of variables to the longitudinal-momentum fractions, ξ i = ω i /ω p ∈ (0, 1), and the total energy ω p ∈ (0, ∞). Straightforward manipulations bring (4.7) to the general form (4.3). Evaluating the left side of (4.4), one finds ∆ p = n i=1 ∆ i − n + 1 as expected, along with the ope factor, ope (n) where B(a 1 , . . . , a n ) = Γ(a 1 ) · · · Γ(a n )/Γ(a 1 +· · ·+a n ) is the multivariate beta function, and the integral has been performed by recognizing that its integrand is a standard Dirichlet distribution ( [44], chapter 49). 5 To see that this evaluation of the OPE singularity matches with (4.6), simply rewrite the beta functions in (4.6) in terms of gamma functions. The product of gamma functions collapses telescopically, yielding a match.
4.2 1 − 2 + . . . n + → p ± Next we come to the case when one negative and n − 1 positive helicity gluons become collinear. Using sequential Wick contractions, this time we arrive at with ∆ p = n i=1 ∆ i − n + 1 the same as before. Here, at leading order, while ope (n) which we have simplified by collapsing individual beta functions to multivariate beta functions. However, we will see that we will also need subleading contributions from descendants for this case. As we did in (4.8), we want to perform a Mellin transform on the splitting functions and check whether they match with the OPE. Mellin transforming split = B(∆ 1 + 1, ∆ 2 − 1, . . . , ∆ n − 1) z 12 z 23 · · · z n−1,n , (4.12) which again matches with (4.10).
To venture beyond such basic consistency checks, we finally come to the Mellin transform of the second splitting function split where S 1j := 1≤k<l≤j ξ k ξ l |z kl | 2 . Firstly notice that this has the same number of terms as (4.11). Each term in (4.13) comes from a particular MHV diagram [35], so that this counting points to a plausible 1:1 correspondence between each term in the multi-gluon OPE and an MHV diagram. We will begin with a detailed exploration of the various terms in this integral in the "toy model" of a triple-collinear limit. Subsequently, we will briefly explain how to scale this up to general n, in particular holographically recovering the j = 2 term in the second line of (4.13). The means of recovering the j ≥ 3 terms and the last term are still being investigated.

n = 3
In the triple-collinear limit, we need to evaluate ope (3) where S 12 = ξ 1 ξ 2 |z 12 | 2 and S 13 = ξ 1 ξ 2 |z 12 | 2 + ξ 1 ξ 3 |z 13 | 2 + ξ 2 ξ 3 |z 23 | 2 . Let's look at the first term in (4.14). Substituting for S 12 and writing z 13 = z 12 + z 23 , it can be converted into The integration has been performed using the integral formula (B.2). The beta functions are the ones we anticipated in (4.11). But we can go even further and holographically predict the entire Gauss 2 F 1 appearing in (4.15). To do this, we compute the 3-gluon OPE by incorporating the contributions from conformal descendants discussed in §3.2. We compute the coefficient of the 1/z 12 z 23 term in the 3-point OPE O − O + O + by using (3.14). This is precisely the term, We only need the leading term in the second Wick contraction because other terms would come with further descendants of O + e ∆ 1 +∆ 2 +∆ 3 −2 and give subleading contributions to the triple-collinear limit. Applying (3.15), (3.20), we see that the terms in C − + + ∆ 1 ,∆ 2 that contain positive powers of∂ 2 cannot act on 1/z 23 , 6 while the action of the rest of the terms produces the coefficient, having used the series expansion given in (B.2). As promised, the expected Gauss hypergeometric function is generated dynamically in the CCFT.
Let us also try to evaluate the second term in (4.14), though we will only partially succeed at matching this with the OPE. The calculation of the OPE in (4.17) motivates us to use the new variables w := z 12 /z 32 and its conjugate as expansion parameters around the leading singularity. These allow us to reexpress the second term of (4.14) as . To manifest some of the hidden structure in this integral, we Taylor expand inw. This results in Evaluating the ξ 2 -integral, each term in this series takes the form of a type-(3, 6) Aomoto-Gelfand hypergeometric function (see [45], section 3.3.5), Such functions were already encountered in the context of celestial amplitudes in [5]. One can perform these integrals explicitly by computing a double series expansion in w,w and hope to match the coefficients with the 3-gluon OPE.
For instance at O(w 0 ), setting w = 0 while formally keepingw non-zero directly in (4.18), this time we arrive at a 2 F 1 inw : This can be recovered by computing the coefficient of the 1/z 12z23 term in the O − O + O + OPE, and again keeping only conformal descendants, Using (3.19), we find This agrees with (4.21). The most plausible origin of the non-trivial functional dependence of (4.20) on w lies in contributions coming from Kac-Moody descendants. This will require finding an all order understanding of these analogous to our analysis of conformal descendants in §3.2.

General n
We can now describe some methods that may help scale up these computations in the future.
Denote the j th term, 2 ≤ j ≤ n − 1, in the expression (4.13) by Observe that, except for the delta function, the integrand of I n,j can be factorized into a product of two functions, one depending on ξ 1 , . . . , ξ j and the other on ξ j+1 , . . . , ξ n . In fact, Euler integrals like I n,j satisfy an elegant factorization property whereby factorization of the integrand also breaks the integral into two smaller such integrals. This is discussed in some detail in appendix B. Using this, the Euler sub-integral involving ξ j+1 , . . . , ξ n is found to be a simple Dirichlet integral. Then (B.9) yields one of the anticipated multivariate beta functions occurring in the various terms of (4.11), leaving us with removing the redundancies in the calculation.
With the help of this simplification, the j = 2 case for instance reduces to a product of beta and hypergeometric functions, This is the all multiplicity generalization of the calculation (4.15) of the same term in the triple-collinear case. Scaling up the Wick contractions in (4.16) by adding n − 3 further positive helicity gluon operators, and keeping conformal descendants for just the first O − O + contraction, CCFT predicts precisely this singularity.
To end this section, let us describe a systematic way to recover the rest of the leading singularities entering the expression (4.11) of ope (n) + (1 − ). Inspired by the variable w = z 12 /z 32 that showed up in the triple-collinear limit, we define a set of new variables in terms of ratios of consecutive distances, To leading order in the ratios w i , Similarly, we find the numerator factor, j l=1 ξ l z 1l 30) and the denominator factors, The main signifance of these expansions lies in the fact that one can also systematically keep subleading terms of O(w i ) to probe descendants exchanged in the OPE. With these leading order results in w i , along with judicious use of the constraint j k=1 ξ k = 1, (4.25) simplifies to Again using the factorization techniques of appendix B, specifically (B.6), this evaluates to the desired leading ope singularity, And finally, one can similarly show that the last term of (4.13) produces the last singularity in (4.11).

Multi-graviton collinear limits
One can hope to attempt similar computations for perturbative gravity. However, the collinear regime of graviton amplitudes is generically non-singular (i.e., not meromorphic). Consequently, the question of multi-collinear limits also becomes much less precise. In fact, generally it even depends on the order in which the gravitons are made collinear. Hence, in this section, we will restrict ourselves to an analysis of sequential double-collinear limits. These are the natural objects that we can expect to get mapped to sequential applications of the OPE under a Mellin transform. They would also have to act as leading order approximations to more precise notions of multi-collinear limits for sake of consistency. 8 For simplicity, we will only consider the case where all helicities are positive and the case where only one helicity is negative. In the first case, on applying the 2-point OPE between two positive-helicity gravitons (2.7) recursively, we have n−1z 12z23 · · ·z n−1,n z 12 z 23 · · · z n−1,n In the case of multi-gluon collinear limits split  − (1 − ), this approximation happens to yield the exact splitting functions as they only possess 2-particle factorization singularities. These also generated OPE coefficients that were permutation symmetric in the positive helicity gluons. So we did not need to discuss these issues earlier.
where ∆ p = n i=1 ∆ i . Unlike the n-point OPE singularity (4.8) for gluons, here the coefficient does not collapse to a multivariate beta function symmetric under permutations. This means that the order in which we perform the recursive OPE is important, and from the bulk point of view it corresponds to obtaining the multi-collinear limit from repeatedly taking double-collinear limits in the same order.
In momentum space, repeated double-collinear limits give the following multi-collinear splitting function, split (n) where ω p = n i=1 ω i , having concatenated the double-collinear splitting function split (2) + given in (2.18) n − 1 times. Along the lines of (4.4), Mellin transforming this gives the conformal weight ∆ p = n i=1 ∆ i and the ope factor, ope (n) where the integral can be performed by noticing that the integrand is an example of a generalized Dirichlet distribution ( [44], chapter 49). This matches the celestial CFT result (5.1). The non-trivial consistency check here is the fact that Mellin transforms do map concatenated splitting functions to sequential OPEs, which is a basic requirement before one can embark on an analysis of BMS descendants.
Similarly, the leading contribution to ope n−1z 12z23 . . .z n−1,n z 12 z 23 · · · z n−1,n with the same weight ∆ p = n i=1 ∆ i as before. In this case, the corresponding momentum space splitting function is found to be split (n) which again matches the celestial CFT result (5.4).

Conclusions
The emergence of bulk physics from celestial CFT is an important subject of much ongoing research. Ideally speaking, given the OPE algebra of the holographic dual, one should be able to work out all its correlators recursively. We have shown that even in the absence of this, we can determine many interesting limits of celestial amplitudes already with the leading order OPE. Our focus has been on finding an understanding of emergent locality and unitarity through the lens of asymptotic symmetries and the celestial operator algebra.
Our methods aim to utilize the operator spectrum of the CCFT to all orders, and hint at interesting organizational principles that could generate multi-particle factorization behavior from the CFT data. Moreover, they also open the doors to many interesting directions of speculation. The operator spectrum of the Yang-Mills CCFT clearly contains much more information than we have been able to find from just translation invariance. Even though we managed to fix the contributions of all global conformal descendants and some leading examples of Kac-Moody descendants to the OPE of celestial gluons, the absence of the remaining Kac-Moody descendants is still a big gap that needs to be filled. We suspect that these extra descendants will help in finding a truly holographic derivation of all the remaining terms in the multi-collinear splitting functions discussed above.
On a similar note, we also need to find how the subleading soft gluon symmetry of [43] fits into this paradigm of primaries and descendants. This should be an interesting representation theoretic problem in its own right. In fact, the subleading soft gluon symmetry will in general impose non-trivial constraints on the multi-gluon OPE, just as it constrained the 2-gluon OPE in [29]. These constraints could take the form of differentialrecurrence equations like the well-known Gauss contiguous relations and give an alternative way of discovering the hypergeometric functions occurring in the Mellin transformed splitting functions. Such a method would also be easier to scale up to higher multiplicity in contrast to our technique of summing up infinite series of descendants.
Another interesting route for future work is the study of the CCFT spectrum dual to general relativity and possibly quantum gravity. Initial steps in this direction have been taken in [30,31] where the OPE coefficients of some of the BMS and other descendants were computed by using symmetry constraints on the celestial graviton OPEs. However, here the set of symmetries that form a non-trivial algebra with translations is much larger, obstructing an all order computation analogous to that in §3.2. But the recent work on a double copy for celestial amplitudes [10] holds promise to overcome these issues. It should be possible to find a notion of color-kinematics duality that acts as an algebra homomorphism on the CCFT operator algebra and maps celestial gluon OPEs to those of celestial gravitons. Hints of this are already present in our results from §3.3. There, if one maps j a −1 O b ∆ 1 +∆ 2 −1 and j b −1 O a ∆ 1 +∆ 2 −1 to P −2,−1 G + ∆ 1 +∆ 2 −1 and −P −2,−1 G + ∆ 1 +∆ 2 −1 respectively in the notation of [30], then the OPE coefficient for the graviton supertranslation descendant P −2,−1 G + ∆ 1 +∆ 2 −1 is found to be α − β. This is given in (3.28) and matches with the result found in section 8.3 of [30]. Also, the vanishing of the OPE coefficient of its antiholomorphic partner P −1,−2 G + ∆ 1 +∆ 2 −1 seems to go hand in hand with the absence of antiholomorphic Kac-Moody descendants in the O + O + gluon OPE. Similar lines of research could also be explored in the full Einstein-Yang-Mills theory.
Above, we also saw that we only possess limited knowledge of the multi-collinear behavior of gravity. Celestial CFT might also be able to help with this by giving a concrete foundation for the gravitational MHV formalism of [46]. We saw indications of this when working out the multi-gluon OPEs in §4. There, various terms in the multi-collinear splitting functions were in 1:1 correspondence with both terms in the CSW recursion relations and the multi-gluon OPE singularities. This might have a straightforward generalization to gravity and help in making novel universal statements.
Finally, we would like to mention that one of our original hopes in deriving subleading terms in the OPEs was to find conformal block expansions for 4-point celestial amplitudes. There has been some work on partial wave expansions in [7,11], but relating them to the operators flowing in the celestial OPE is still an open question of great import. This will bring us a step closer to viewing scattering amplitudes as conformal correlators.

B Factorization of Euler integrals
We first note some useful integral representations of special functions. The Euler beta function is given by The Gauss hypergeometric function can be represented by a similar Euler integral formula as well as by a series expansion, Such integrals occur numerous times in §4. Next, we provide some standard tricks that can be applied to recursively simplify the Euler-type integrals occurring in this work. Suppose we start with an integral of the form ξ j   f (ξ 1 , . . . , ξ k ) g(ξ k+1 , . . . , ξ n ) , (B.3) for a pair of integrable functions f and g. Also assume that g is a homogeneous function of degree β under a diagonal rescaling, g(t ξ k+1 , . . . , t ξ n ) = t β g(ξ k+1 , . . . , ξ n ) , t ∈ R * .

(B.4)
This will be a ubiquitous property in our splitting functions. Since the dependence on ξ i 's is factorized in the integrand, we now show that this is also enough to factorize the integral into two smaller Euler integrals. Simply insert identity in the form, Since 0 < n i=k+1 ξ i < 1 due to the delta function constraint in (B.3), we only need to integrate ξ 0 over this range. Inserting this in I[f · g] and rescaling ξ i → ξ 0 ξ i for i = k + 1, . . . , n produces  Moreover, if f is also a homogeneous function of degree α under a diagonal rescaling, say f (t ξ 1 , . . . , t ξ k ) = t α f (ξ 1 . . . , ξ k ) , t ∈ R * , (B where, as expected, This factorization will be very useful in extracting beta functions from complicated integrals. Such properties may also be helpful in studying more general factorization behaviors of celestial amplitudes in the future.