Celestial double copy from the worldsheet

Using the ambitwistor string, we compute tree-level celestial amplitudes for biadjoint scalars, Yang-Mills and gravity to all multiplicities. They are presented in compact CHY-like formulas with operator-valued scattering equations and numerators acting on a generalized hypergeometric function. With these we extend the celestial double copy to tree-level amplitudes with arbitrary number of external states. We also show how color-kinematics duality is implemented in celestial amplitudes and its interpretation in terms of a generalized twisted cohomology theory.


Introduction
The perturbative S-matrix of quantum field theories can have many interesting properties and relations. If we take these relations to reflect some fundamental property of the respective QFTs, we expect them to survive in some form after a change of basis in the Hilbert space. The S-matrix is usually computed in a basis of plane waves for external particles due to its simplicity, but other bases might be more useful in answering certain questions. In [1][2][3][4] it was argued that a basis of external particles transforming as conformal primaries on the sphere at the null boundary of Minkowski space is the most appropriate one to study a potential holographic CFT living on this sphere. Amplitudes with conformal wavefunctions as external states were dubbed celestial amplitudes. These can be obtained from amplitudes computed in the usual plane wave basis by a Mellin transform. This has the effect of mixing the UV and the IR. Yet universal properties of amplitudes, such as soft limits and IR factorization [1,[5][6][7][8][9][10][11][12][13][14][15][16], survive this change of basis albeit in a different guise. Recently, it was shown in [17] that relations between the S-matrices of Yang-Mills and gravity known as the double copy also survives this change of basis, at least up to four point amplitudes. 39]. But computations still relied heavily on the leftover momentum conservation and special properties of plane wave backgrounds. Ambitwistor string computations have also provided new formulas for amplitudes in AdS spacetimes [40,41] which have a structure very close to the one we find for celestial amplitudes. We would like to find a notion of double copy that can be more easily generalized to curved spaces and valid for a larger class of spacetimes including asymptotically AdS spacetimes.

Summary of results:
• We compute compact, n-point formulas for tree-level celestial amplitudes of biadjoint scalars, gluons and gravitons. They are given in a universal form as operator-valued numerators acting on the Mellin transform of the scalar contact vertex. Analogous to the case with external plane waves we can represent the amplitudes either as a sum over trivalent graphs (3.19), or as integrals over the moduli space of n-punctured Riemann surfaces localized to the operator-valued celestial scattering equations (3.10).
• As a consequence of the computation given above, we prove that the celestial double copy introduced in [17] is valid for amplitudes at all multiplicities.
• We introduce a generalization of twisted cohomology to operator-valued twisted forms which are the relevant objects in the CHY formulas for celestial amplitudes. We use this, together with a straightforward generalization of the results in [33], to define color-kinematics duality for celestial amplitudes and show how to obtain colorkinematical dual numerators from the ambitwistor string numerators.
This paper is organized as follows: we start in section 2 by reviewing some facts about celestial amplitudes, CHY formulas, double copy and twisted homology to make this paper more self-contained. The reader familiar with these topics can safely skip these sections. Section 3 contains the formulas for celestial amplitudes with explanations about its constituents. Next, in section 4, we introduce the ambitwistor model and go through the calculation used to obtain the formulas introduced previously. In section 5 we introduce a generalization of the twisted cohomology to operator-valued forms, give a definition of color-kinematics duality for celestial amplitudes and show how numerators obtained from the ambitwistor string naturally obey this duality. We finish with some discussion about this framework and possible further generalizations in section 6.

Review
In order to make this paper self-contained we quickly review in this section some background on recent technology used in the study of amplitudes. We keep the reviews short and focused on what is needed in the rest of the paper. In subsection 2.1 we review some facts about celestial amplitudes, in subsection 2.2 CHY formulas and the scattering equations are introduced. In subsection 2.3 we recall the basics of color-kinematics duality, and in subsection 2.4 we review how the CHY formulas and double copy are to be interpreted in light of twisted cohomology on the moduli space of punctured Riemann spheres.

Celestial amplitudes
Celestial amplitudes are obtained by Mellin transforming the usual momentum space amplitudes to make manifest their transformations under the conformal group of the celestial sphere at null infinity I + . This is essentially a change of basis [3] on the Fock space from the plane wave basis into a basis of conformal primary wavefunctions 1 Let X be a point in R 1,d+1 , with D := d + 2, and denote the massless scalar plane wave as e isωq·X where ω denotes its energy, q ∈ S d is a direction on the celestial sphere, and s = ±1 denotes if the particle is outgoing or incoming. The scalar conformal primary wavefunction is obtained by Mellin transforming the energy ω, giving an external wavefunction that lives on the celestial sphere with a new quantum number: the conformal dimension ∆. Conformal wavefunctions for external gluons and gravitons also exist but here we will work directly with the Mellin transform of the plane waves, a a,s which are given in terms of the scalar wavefunctions, the color generators T a , and polarizations µ ,˜ µ that only depend on q µ [4]. These wavefunctions transform as conformal primaries up to gauge transformations [3]. Wavefunctions for massive particles have also been worked out [42,43] but we won't have anything to say about massive particles. Amplitudes are in principle computed using (2.1) and (2.2) as external wavefunctions, but it is more convenient to simply Mellin transform the usual amplitudes computed with plane waves. A celestial amplitude is formally given by Explicitly carrying out the Mellin transforms can quickly become unwieldy as the number of external particles increases. In the following section we'll introduce the CHY representation of amplitudes which gives a very compact formula for massless n-point amplitudes.

CHY formulas
The CHY formulas present the D-dimensional, tree-level massless S-matrix of several quantum field theories [27,28,44] as an integral formula. A generic CHY formula is written as The integral is taken over the moduli space of n-punctured Riemann spheres M 0,n which carries an action of the group SL(2, C) denoted by the factor of 1/vol SL(2, C) in the measure. Concretely, this implies we can fix the position of any 3 punctures, for example {z 1 , z 2 , z 3 }, removing the associated differentials dz i and introducing the Jacobian gives the S-matrix for a cubic biadjoint scalar. Another interesting numerator is the reduced Pfaffian, with Pf the Pfaffian of the 2n × 2n matrix, with two lines and two columns removed, denoted by Ψ ij ij . Its components are given by the matrices Taking I = Pf Ψ n andĨ = PT n as numerators the CHY formula gives the S-matrix for external gluons in Yang-Mills. Using a Pfaffian for both, I = Pf Ψ(k, ) andĨ = Pf Ψ(k,˜ ), gives the S-matrix for gravitational amplitudes in NS-NS gravity. The universality of the CHY representation is due to the presence of the scattering equations E i [45], which are the last ingredient of (2.4) to be explained. They appear as the arguments of the delta functions in iδ (E i ). These delta functions are taken as holomorphic delta functions, that isδ(E i ) =∂(z) E i , effectively fixing the contour of integration to the solutions of the scattering equations. The symbol means that delta functions for three equations should be omitted from the product and another Jacobian of the form The scattering equations themselves are The delta functions impose them as constraints on the z i 's, completely localizing the integration over M 0,n . There are generically (n − 3)! points in M 0,n which solve the scattering equations. Denoting the solutions of the scattering equations by σ i , the amplitudes above can be written as with Φ = det ∂ j E i , the Jacobian coming from the delta functions.
In the CHY formulas the sum over trivalent graphs is replaced by the fully localized integral over M 0,n , making manifest properties which are hard to see from the Feynman diagrammatic expansion. CHY formulas also exist for many theories beyond the ones reviewed above [44,46]; many can be seen as originating from an unconventional string theory, called the ambitwistor string [26,47]. The advantage of having a worldsheet theory like the ambitwistor string is the ease of generalizations to other settings, e.g. higher loops [48][49][50][51][52][53] and curved backgrounds [35,[54][55][56][57][58], which in turns leads to new formulas for amplitudes in these settings. In section 4 we'll show how CHY-like formulas for celestial amplitudes originate from the ambitwistor string.

Color-kinematics duality and double copy
The original double copy [31] is a prescription to obtain gravitational amplitudes from suitable squares of Yang-Mills amplitudes. There are now many pairs of theories known to have amplitudes related by a double copy prescription, as well as several proposals for double copy of non-linear solutions. See [59] for a recent review of the field. Here we'll present the aspects of the original double copy which we'll need in the rest of the paper. Given a tree-level Yang-Mills amplitude it can be represented as a sum over trivalent graphs by opening up four point interactions in the Feynman diagrams, (2.11) Here Γ is the set of trivalent graphs, c γ are color numerators carrying the gauge group information, n γ are kinematical numerators which are polynomials in the external momenta and polarizations, and P e are the propagators associated to each edge e of the graph γ. Color numerators corresponding to graphs which differ only by BCJ moves on some subgraph, see figure 1, obey identities inherited from the usual Jacobi identity, If the four point vertices are opened up appropriately then, at tree-level, kinematical numerators can always be found such that they satisfy identities analogous to the ones satisfied by the color numerators. That is they satisfy the Jacobi-like relation, n γs + n γt + n γu = 0 . (2.13) These kind of numerators are called color-kinematical dual numerators [32]. Substituting in (2.11) the color numerators by a set of color-kinematical dual numeratorsñ γ (obtained γ s γ t γ u + + = 0 Figure 1. Graphs related by BCJ moves from n γ by replacing with˜ ) yields another amplitude, which turns out to be a gravitational amplitude. The requirement that numerators obey color-kinematics duality ensures that (2.14) is invariant under linear diffeomorphisms. This is the double copy prescription. The CHY formulas reviewed in the previous section also have a double copy structure. There, the analogue of color numerators is the Parke-Taylor factor PT and the analogue of kinematical numerators is the Pfaffian Pf Ψ. The sum over trivalent graphs is replaced by an integral over the moduli space and double copy becomes a simple substitution rule. Color-kinematics seems to be disconnected from the double copy in this case, but they turn out to be intimately related through twisted cohomology on M 0,n . We review this in the following subsection.

Twisted cohomology
The ingredients in the CHY formula have an interesting interpretation in terms of a cohomological theory on the moduli space of punctured Riemann spheres M 0,n . Here we'll quickly go over the relevant details of the constructions in [33,60] that we'll generalize in section 5. Proofs of the statements presented below, as well as details of the computations in the context of amplitudes can be found in [33,60,61]. Other mathematical details can be found in the original mathematical literature [62][63][64][65].
We start by defining a meromorphic one-form in M 0,n using the scattering equations: where the prime in denotes that three points have been fixed to account for the SL(2, C) invariance. Using momentum conservation one can show that the form (2.15) doesn't depend on which three points were fixed and that it only has simple poles along the boundaries of M 0,n . With this in hand we define the twisted de Rham operators (2.16) Both are flat connections [∇ ± , ∇ ± ] = 0 on certain line bundles L and its dual L ∨ over M 0,n called local systems. Sections of L (L ∨ ) are given by functions on M 0,n that are locally covariantly constant in ∇ + (∇ − ). The operator ∇ + acts naturally on Ω • (M 0,n , L), the complex of differential forms on M 0,n with coefficients in L, and squares to zero (∇ + ) 2 = 0. We use it to define twisted cohomology groups, Generically, i.e. for generic momenta, the only non-vanishing group is the middle-dimensional one H n−3 (M 0,n , L) with (n − 3)! independent generators. There is an analogous construction for cohomology groups with values on L ∨ which we omit.
Taking dµ a top holomorphic form on M 0,n , the numerators ϕ + = I dµ and ϕ − =Ĩ dµ can be interpreted as elements of the twisted cohomology groups ϕ + ∈ H n−3 (M 0,n , L) and ϕ − ∈ H n−3 (M 0,n , L ∨ ). Moreover, amplitudes in the CHY formalism can be interpreted as a bilinear pairing between these cohomology groups: giving the intersection number of ϕ + and ϕ − . An alternative evaluation of this pairing can be given where the contours around the scattering equations are deformed, picking up contributions only from the boundaries of M 0,n . A codimension 1 boundary is reached when the n-punctured sphere degenerates into a nodal surface given by two punctured spheres connected by a node. Higher codimension boundaries are reached with further dengenerations of these spheres into surfaces with more nodes. The deepest boundaries arise when all that is left is a nodal surface composed of spheres with three marked points. These spheres are connected through these marked points in such a way that no closed cycle can be drawn that passes through the nodal points. That is, the nodal surface resembles a tree-graph with only trivalent vertices. In fact, one can label these deepest boundaries by trivalent trees whose punctures are distributed along the external edges. Different boundaries correspond to different assignments of labels for the external punctures.
The end result is that the contour given by the scattering equations can be deformed to pick contributions from the deepest boundaries of the moduli space. This presents the amplitude as a sum over the boundaries of M 0,n labelled by trivalent graphs, In this expression, are residues of the top meromorphic forms ϕ ± along the boundary divisor v γ labelled by the trivalent graph γ. This representation is similar to the field theory one (2.11) and (2.14) but here the numerators are guaranteed to obey color-kinematics. This follows from global residue theorems on M 0,n as explained in [33]. The general argument can be summarized as follows: take a triple of trivalent graphs γ s , γ t , γ u differing only on a subgraph connecting four edges as shown in figure 1. These three divisors can be seen as arising from the same corner of M 0,n where a four-punctured sphere degenerates as one of its punctures, z, approaches one of the other three punctures, z s , z t , z u , generating the three boundary divisors associated to the trivalent graphs γ s , γ t , γ u , see figure 2. We model the neighborhood of these degenerations by a thrice-punctured sphere Σ 3 with coordinate z and the three marked points z s , z t , z u fixed to some values using the SL(2, C) symmetry. Near these degenerations the differential form ϕ + restricts to a 1-form on Σ 3 with poles along the marked points. Kinematical numerators are given by residues of this form a ∈ {s, t, u} . (2.21) and are related by a linear identity due to the global residue theorem on Σ 3 , While this argument would work for any regular differential forms on M 0,n , only forms which are elements of the twisted cohomology groups actually generate field theory amplitudes. Moreover, many simplifications occur by choosing good representatives for these cohomologies.

Celestial scattering equations
The CHY formulas reviewed in the previous section have a natural celestial analogue. These are most straightforwardly -if somewhat formally -derived by directly Mellin transforming the momentum space expressions. We start this section by performing these Mellin transforms and writing down explicit formulas for celestial amplitudes. These take the form of Gelfand A-hypergeometric functions in all dimensions, now governed by the celestial scattering equations. Some example computations are also provided.

Amplitude formulas
The n-point celestial amplitude can be written down as a Mellin transform of the momentum space formula (2.4) [8], The primary trick to simplify this is to perform as many ω j -integrals as possible against the momentum conserving delta function. We can already do this at the level of (3.1), but we will find it much more illuminating to first manipulate it a bit. For this we make use of the momentum generators in Mellin variables. The action of the momentum operator P αα on a function A(k) of null momentum k = s ω q is given by a trivial multiplication Its action on the Mellin transform of A(k) is then easily expressed as This dictates the definition of the celestial translation symmetry generators [66], which are operator-valued null vectors, K 2 i = 0. Momentum conservation is then equivalent to invariance under the diagonal translation generator, that is, n-point celestial amplitudes should be annihilated by Indeed, acting with this on (3.1) produces a factor of i s i ω i q αα i inside the Mellin transforms which vanishes by momentum conservation. So, at least whenever the Mellin transforms converge (or are understood distributionally), celestial amplitudes are invariant under diagonal translations.
Hence we can make the formal replacements, inside the scattering equations and the CHY integrands. By converting them into operators we can take these objects outside the Mellin integrals. Moreover, the various K i 's clearly commute with each other and there is no operator ordering ambiguity. We thus find the following expression, having defined the celestial scattering equations, (3.8) and the Mellin transformed n-point contact diagram, Rigorously speaking, replacements like (3.6) are supposed to be made inside analytic functions. In section 4, we will justify performing these replacements inside the delta functions δ(E i ) by deriving it from the ambitwistor worldsheet CFT. S n is the contribution of a φ n contact term to the n-point celestial amplitude of a theory of massless scalars φ. In appendix A we explicitly evaluate these contact diagrams for arbitrary multiplicity and dimension, citing here only the final expression for the celestial CHY formulas, In particular, the replacement (3.6) applied to the Parke-Taylor factor (2.5) and the Pfaffian (2.6) provides the celestial equivalent of kinematical numerators in (3.7) for biadjoint scalar, Yang-Mills and gravity celestial amplitudes. Using the celestial scattering equations, these take the form of operator-valued numerators acting on the same universal function The function F n ({∆, q, s}) takes different forms depending on whether n > D or n ≤ D.
n > D : In this case, there are more ω j integrals than momentum conserving delta functions. This allows us to perform D of the Mellin integrals -say those over ω 1 , . . . , ω D -and the rest simplify to Here, the various indices run over l = 1, 2, . . . , D while r, r , r = D + 1, . . . , n. The coefficients u and u lr in the integrand can be expressed in terms of D × D minors of the D × n matrix (q µ 1 , q µ 2 , . . . , q µ n ) of celestial positions. Define the determinants, where ε µ 1 µ 2 ...µ D is the D-dimensional Levi Civita symbol. Then we can express them quite compactly as (3.14) The Heaviside step functions constrain the u lr to be positive. The integrals over the ξ r 's produce Aomoto-Gelfand hypergeometric functions over the Grassmannian Gr(n − D, n) (or equivalently on Gr(D, n)) [67,68]. Such integrals were first identified in the context of celestial amplitudes in [18].
n ≤ D : In this case, the result is distributional with D − n + 1 leftover delta functions, Here, l = 1, 2, . . . , n − 1. The spacetime index µ has also been partitioned into two sets: Then we have localized the Mellin integrals on the delta functions imposing i s i ω i q a i = 0 by solving for the ω l in terms of ω n . The various coefficients in the above expression are where the determinants (i 1 i 2 . . . i n−1 ) are now (n−1)×(n−1) minors of the (n−1)×n matrix (q a 1 , q a 2 , . . . , q a n ). Again, these coefficients satisfy positivity constraints u ln > 0 for all l. The remaining delta functions impose momentum conservation in the "transverse" D − n + 1 dimensions. The result is still Lorentz invariant, the Mellin transforms themselves preserve the symmetry but in performing them explicitly, non-Lorentz covariant choices had to be made.
We also remark that the leftover Mellin integral over ω in (3.10) is generically divergent. An interpretation of such integrals was given in [69] in terms of "generalized delta functions", allowing one to declare Formally, we see for instance from (3.15) that we indeed need the condition i ∆ i = D to hold so that it has conformal weight ∆ n in q µ n . If the final amplitude is finite or marginally convergent, this awkwardness gets resolved by the application of the K i 's on S n as this shifts the conformal weights appropriately. We discuss such shifts below.
Let us now study some salient features of our formulas. First note that the various scattering equations E i trivially commute since the K i commute. Moreover, even though momentum conservation is absent, they are still SL(2, C) invariant when acting on translation invariant objects. This is a consequence of the fact that the residue of any E i at z i = ∞ is given by were we used K 2 i = 0. This residue is proportional to the diagonal translation generator which annihilates the contact diagram.
Having evaluated the Mellin transform of S n , we can use (3.10) to solve the scattering equations and obtain the celestial S-matrix of a variety of theories. Since the operators K i formally behave as commuting numbers, the scattering equations are solved by the same algebraic expressions coming from the (n−3)! solutions in momentum space. The worldsheet path integral implements a sum over all of these solutions and generates trivalent graphs γ contributing to the n-point amplitude. The resulting celestial amplitude has the following general structure, Each internal edge e of a graph γ now comes with an operator-valued propagator denominator, The numerators N γ andÑ γ are also theory dependent operators given by residues of the CHY numerators I andĨ analogous to eq. (2.20), see also section 5, which act on the celestial contact amplitude. For n > D, one can take these operators inside the Mellin integrals in S n . For i = r, one straightforwardly replaces e ∂ ∆r → ξ r , while for i = l one replaces e ∂ ∆ l → r u lr ξ r as expected. Simultaneously, one also needs to appropriately shift i ∆ i occurring in the distributional factor (3.17), and this is most easily seen by working through the examples given in the next subsection. The resulting integrals yield Gelfand A-hypergeometric functions [70,71]. Similarly, for n ≤ D, one can replace factors e ∂ ∆ l → u ln and e ∂ ∆n → 1, while again shifting i ∆ i appropriately.

Examples
To illustrate our methods, we work out some examples of biadjoint scalar amplitudes in different dimensions. Generalizing to gluons and gravitons is computationally tedious but straightforward.
4 points. In this case, we only have one scattering equation. We choose to fix the three points z 1 , z 2 , z 3 so that the remaining scattering equation corresponds to z 4 : Fixing the global conformal symmetry also introduces a Faddeev-Popov determinant (z 12 z 23 z 31 ) 2 into the worldsheet integrals, where z ij := z i − z j . Using 2πiδ(E i ) =∂ i (1/E i ), we then need to evaluate For convenience, we take the fixed points to be z 1 = 0, z 2 = 1, z 3 → ∞. As explained before, since the K i 's commute, we can formally solve the scattering equations as in momentum space. We find . (3.23) More precisely, this makes sense via its action z * 4 ·S 4 on the contact diagram which produces an ordinary number. To obtain particular amplitudes we can insert various different numerators. For example, we can compute the biadjoint scalar amplitude in the color ordering (1234|1234). This is done by substituting the (color-stripped) Parke-Taylor factors, for IĨ. Using 4 i=1 K i · S 4 = 0, one easily finds

(3.25)
Clearly, all that the derivatives do is shift some of the weights by −1. This phenomenon is also clear from the perspective of the plane wave basis. There the propagator denominators contain factors of ω i 's which induce precisely these shifts. Eq. (3.25) can be further simplified on each channel. When D < 4, using (3.12) and (3.17), the s-channel contribution becomes Other channels contribute similar terms.
Propagators at higher points. At n points, each Feynman diagram of the cubic biadjoint scalar theory comes with n − 3 propagators and trivial numerators. Clearly, each propagator denominator induces a shift i ∆ i → i ∆ i − 2 in the Mellin integral (3.17). Altogether, we find a factor, in the final amplitude. With D = d + 2, the on-shell phase space of massless particles in R 1,D−1 is spanned by conformal basis states on the principal continuous series ∆ i ∈ d 2 + i R [3]. So the distribution (3.28) is a true delta function precisely for d = 4, i.e., in six dimensions. It is curious to note that D = 6 is precisely the dimension in which the biadjoint scalar theory is a classical CFT. 2 In general, we will have propagator denominators acting on the celestial contact diagram. As a simple example, consider the action of a single propagator acting on S n for n > D, Due to the new quadratic polynomial in the integrand, the result now produces a Gelfand A-hypergeometric function [70,71]. Every propagator is taken care of by insertion of such quadratic polynomials. We may also find the insertion of higher degree polynomials through the numerators in gluon and graviton amplitudes, which are also handled by the same theory of hypergeometric integrals. Similar integrals have also recently occurred in the physics of loop level Feynman integrals [72] as well as stringy canonical forms [73].

Models
In this section, we derive our worldsheet formulas for celestial amplitudes using ambitwistor string theories [26]. These are chiral string theories with target the space of null geodesics in Minkowski space. The worldsheet action of a large class of such models takes the general form, Here, X µ and P µ are fields of conformal weight (0, 0) and (1, 0) respectively on the string worldsheet Σ, and S matter is the action of a pair of auxiliary worldsheet CFTs. T denotes the weight (2, 0) stress tensor, while µ and e are weight (−1, 1) Beltrami differentials gauging the constraints T = 0 and H := P 2 /2 = 0 respectively. The latter generates the gauge symmetry δX µ = αP µ , δP µ = 0 , δe =∂α , where α is a weight (1, 0) vector field on Σ. Thus, the constraint P 2 = 0 and associated gauge redundancy X ∼ X + αP reduce our target space to the space of null geodesics: ambitwistor space. The choice of matter action S matter gives rise to amplitudes of a variety of theories: Biadjoint scalar : S matter = S g + Sg , Yang-Mills : Gravity S g and Sg denote actions of current algebra CFTs corresponding to Lie algebras g andg. We will denote their weight (1, 0) currents as j a andã respectively. In the Yang-Mills and gravitational cases, the fields ψ µ ,ψ µ both denote weight ( 1 2 , 0) worldsheet fermions, while χ andχ are weight (− 1 2 , 1) fermionic Lagrange multipliers gauging fermionic worldsheet symmetries akin to supersymmetry. Worldsheet correlators of the current algebra systems give rise to the Parke-Taylor type numerators (2.5), while those of the fermions produce the Pfaffian type numerators (2.6).
Since both the XP system and matter actions are free CFTs, it is easy to find their fundamental OPEs. In the former case, one finds Similarly, the fermionic OPEs read, Lastly, the current algebra OPEs are the standard ones, where f abc and fãbc are the structure constants of g andg respectively. The levels of the current algebras can be non-zero, producing multi-trace amplitudes, but since we're only going to be interested in the single trace contributions we omit terms proportional to the level. As in standard string theory, on gauge fixing we will also add ghost fields for each of the gauge symmetries. For the bosonic symmetries generated by T and H, one adds two bc ghost systems consisting of fermionic ghosts b,b with weights (2, 0) and c,c with weights (−1, 0). For the supersymmetries generated by ψ · P andψ · P , one adds βγ systems with bosonic ghosts β,β of weights ( 3 2 , 0) and γ,γ of weights (− 1 2 , 0). These fields also have the well-known OPEs, and (4.10) The expression for the BRST operator will vary depending on the matter content of the ambitwistor string. For the tree models we discuss, biadjoint scalar, Yang-Mills and gravity, their BRST operators square to zero with appropriate choices of dimensions and current algebra level. When setting up the model we implicitly assume to be working with choices that render the BRST charge nilpotent, but our final expressions after all the worldsheet calculations have been performed are actually valid in any spacetime dimension. They are after all, tree-level amplitudes of field theories, which are the same independent of the spacetime dimension. More details on the BRST charges can be found in [26]. Next, we construct vertex operators for these theories in the conformal primary basis of states. It is easiest to start with the biadjoint scalar states. Their fixed vertex operators are given by where i is a particle label and φ i (X) ≡ φ s i ∆ i (X; q i ) denotes the scalar conformal primary wavefunctions of (2.1). Generators of the Lie algebras g andg are denoted by T a andTã respectively. Integrated vertex operators are obtained form the usual descent procedure.
Similarly, the fixed vertex operator of a conformal primary gluon external state in Yang-Mills is where the superscript −1 stands for picture number. To construct the corresponding picture number 0 vertex operator, we descend by computing the OPE of V −1 i with the picture changing operator Υ = δ(β) ψ · P . Using (4.6) and the scalar wavefunction (2.1), it is easily seen that Such OPEs are the means by which the celestial translation generators K i µ = s i q i µ e ∂ ∆ i of (3.4) will enter our analysis. They act on the wavefunctions to their right. This produces the picture number 0 vertex operators, (4.14) A similar analysis holds for gravitons, yielding the conformal basis vertex operators, at picture number −1, and at picture number 0. We reiterate that, just like ordinary momenta, the K i 's commute with each other and there is no ordering ambiguity here. BRST closure of these operators follows immediately from the conformal primary representatives being on-shell.

Worldsheet correlators
We start this section by computing the correlators of the biadjoint scalar vertex operators (4.11) in detail and deriving the celestial scattering equations. After this we give the generalization to Yang-Mills and gravity. To compute n-particle amplitudes we quantize these models on a n-punctured Riemann sphere with punctures located at z 1 , . . . , z n . We work in conformal gauge µ = 0 and also gauge fix e to be an element of the (n − 3)-dimensional Dolbeault cohomology group H 0,1 (Σ, T Σ (z 1 + · · · + z n )). The gauge freedom (4.2) is precisely the freedom in choosing a representative of its cohomology class. Explicitly, we pick the gauge fixing condition, (4.17) having chosen the punctures z 4 , . . . , z n as our moduli without loss of generality. The {e i } denote a standard basis of this cohomology group. They act against quadratic differentials like H by picking 2πi times their residues at the z i . So for instance we will find an insertion of inside the path integral for any correlator, coming from the action (4.1). The r i provide n − 3 moduli that are left to be integrated over after the gauge fixing. The gauge fixed action of the biadjoint scalar ambitwistor string is given by S = 1 2π Σ P ·∂X + b∂c +b∂c + S g + Sg . Using this action, the n-particle celestial amplitude is computed by the correlator, where the integral is performed over an appropriate middle-dimensional contour Γ in T * M 0,n [74]. 3 In this expression, we have inserted a product of n − 3 picture changing operators, needed to soak up fermionic zero modes and produce the measure on M 0.n . Their OPEs with fixed vertex operators produce integrated vertex operators corresponding to a choice of n − 3 z i 's. The other three puncture locations z 1 , z 2 , z 3 have been fixed using the residual SL(2, C) symmetry. One computes the products of B i andB i with V i using the OPEs (4.9) to find This strips off the c andc ghosts from n − 3 of the vertex operators, and correlators of the remaining ghost zero modes produce a factor of (z 12 z 23 z 31 ) 2 that can be accommodated by inserting a factor of 1/vol SL(2, C) 2 in the sense of Faddeev-Popov. Next, the current algebra correlators generate a pair of Parke-Taylor numerators PT n and PT n of the form (2.5) ( PT n contains traces over products ofTã i ). We throw out the multitrace terms since we're interested only in amplitudes where Yang-Mills states are exchanged. All these simplifications leave us with a correlator of the XP CFT, (4.23) We evaluate the last correlator by utilizing the OPE (4.13) computed before. It is easily shown that H(z) = 1 2 P 2 (z) acts on a product of scalar wavefunctions φ i (X) by the OPE, As a result, on performing the XP path integral, H(z) is frozen to its "classical" value, .
Note that the double poles dropped out in this computation due to K 2 i = 0. Finally, the exponential e − 1 2π Σ e H can be brought outside the correlator. Using (4.18), the r-integrals subsequently give rise to the n − 3 scattering equations, where Res This justifies the replacements ω i → e ∂ ∆ i done within the scattering equations in section 3. The result is a top-form integrated over M 0,n . The remaining correlator over the XP system contains only the insertion i φ i (X(z i )). Performing the path integral over P µ imposes its equation of motion following from the gauge fixed action (4.19),∂X µ = 0. On Σ = CP 1 , this reduces the path integral over X µ to an integral over its constant zero mode. Denoting this zero mode again by X µ , we find (4.28) Translation invariance of S n and commutativity of the K i 's guarantees permutation invariance of the scattering equations, so all choices of n − 3 scattering equations are equivalent.
An analogous story holds for Yang-Mills and gravitational amplitudes. For instance, in the former case we need to compute a correlator of the form, having chosen to insert picture number −1 vertex operators for particles 1 and 2 without loss of generality. The path integrals over the fermionic ghosts is done as usual and we pull out the insertion of e − 1 2π Σ e H from the correlator as before leading to the celestial scattering equations. The path integral over the currents j a generates a single Parke-Taylor factor plus multitrace terms we once again ignore. This leaves the correlators of the XP and ψ systems along with the βγ ghosts to be performed, (4.31) Following [26], these yield a Pfaffian-type numerator (2.6), now containing the operators K i in place of ordinary momenta k i . The final result is the formula, encoding gluon celestial amplitudes in arbitrary dimensions. For gravity, the only new ingredient is the replacement of the current algebra system with a second fermionic system. The correlator of interest is given by There are two unifying features of all these formulas. The first is the presence of the same contact amplitude S n and celestial scattering equations governing the three expressions. The second is the manifest double copy structure at the level of the CHY integrands given by a simple replacement rule. In the next section, we return to this point with the machinery of twisted cohomology as applied to these operator-valued integrands.

Celestial color-kinematics duality
We have seen in the previous section how the celestial versions of the CHY formulas still manifest a double copy structure for the operator-valued numerators. This can either be stated as a substitution rule for ambitwistor integrands or in terms of their residues and trivalent graphs. In the latter method, gauge invariance of the double copied amplitude is not manifest even though we know it is gauge invariant by construction. As reviewed in section 2 gauge invariance of the double copied amplitude is guaranteed if the kinematical numerators obey color-kinematics duality. We show below that an analogous requirement holds for the operator-valued numerators of celestial amplitudes obtained from the ambitwistor string. The proof will closely follow [33] by recasting the operator-valued numerators as elements of a generalized twisted cohomology.
Let C n denote the configuration space for the celestial data {∆ i , q i , s i } of n-particle celestial amplitudes. The operators K i defined in (3.4) map smooth functions on C n to themselves by shifting their conformal dimensions. These operators are not derivations as they don't satisfy Leibniz rule, but they are linear and commute with each other. Denote the C ∞ (C n )-algebra generated by the K i 's by C. It forms a subalgebra of Hom(C ∞ (C n ), C ∞ (C n )). Even though translation invariance is not manifest it is still a symmetry of celestial amplitudes, so our operators K i will only ever act on functions S satisfying Without loss of generality, we quotient C by the ideal generated by i K i to define as the reduced kinematical space where celestial amplitudes live. In what follows, we will treat the operator-algebras C and K as infinite-dimensional vector spaces spanned by monomials in K i 's. With some abuse of notation, we also denote by K the infinite-dimensional trivial vector bundle M 0,n × K. In analogy with the usual twisted cohomology, we define the twisted de Rham operators with d the exterior derivative on M 0,n and the connection one-form built out of the celestial scattering equations, The operator-valued connections ∇ + and ∇ − act via straightforward multiplication on forms valued in K, denoted by Ω • (M 0,n , K). 4 We also formally define the line bundles L and L ∨ whose sections are K-valued functions on M 0,n . 5 The connections ∇ ± square to zero so we use them to define twisted de Rham cohomologies with coefficients in L and L ∨ , denoted by H • (M 0,n , L) and H • (M 0,n , L ∨ ). 4 The images of i Ki under ∇± are again in i Ki since the Ki operators multiply as usual, commute with each other and with d. 5 More precisely, L and L ∨ are local systems given by an abelian operator-valued representation of π(M0,n). We leave the study of analytic aspects of these definitions to future work.
The operators K i , viewed as elements of K, obey the same algebraic identities as the usual momenta for plane waves k i . That is, they are null, commute with each other and obey a "momentum conservation" identity. Due to this, the analysis in [33,60] carries over to the celestial case with the appropriate substitutions. Let ϕ ± be two K-valued meromorphic top forms on M 0,n with singularities only along its boundary divisor. Take these to also obey ∇ ± ϕ ± = 0 so that they are representatives of cohomology classes in H n−3 (M 0,n , L) and H n−3 (M 0,n , L ∨ ). Explicit examples are the ambitwistor numerators (2.5) and (2.6) with the replacement k i → K i , multiplied by the top holomorphic form dµ on M 0,n .
We define intersection numbers ϕ + |ϕ − in analogy with [60] but in our case this pairing is not a number but an operator. There are two common ways to compute these intersection numbers: one is given by the CHY formulas, as a moduli space integral localized to the solutions of the scattering equations; the other by evaluating the paired forms ϕ ± near the highest codimension boundaries of M 0,n 6 labelled by trivalent graphs. The first method recovers the formula (3.10) for numerators ϕ + = I dµ and ϕ − =Ĩ dµ. The latter gives the amplitude as a sum over trivalent graphs γ, with edges denoted by e. Residues are taken along the boundary divisor v γ associated to the trivalent graph γ. Each edge has an associated denominator P e = ( i∈e K i ) 2 analogous to propagators acting as inverse operators. This is nothing other than eq. (3.19) which was obtained by deforming the contours prescribed by the scattering equations.
With the framework introduced above it is straightforward to adapt the arguments of [33] to show how color-kinematics is implemented for celestial amplitudes. The set up is analogous to section 2.4. Take a triple of trivalent graphs γ s , γ t , γ u differing only on a subgraph connecting four edges as shown in figure 1. The boundary divisors for these graphs arise from a four-punctured sphere degenerating as one of its punctures z approaches one of the other punctures z s , z t , z u . The neighborhood of this degeneration is modelled by a sphere with three punctures, Σ 3 , with coordinate z and fixed marked points z s , z t , z u , see figure 2. Take ϕ M as a K-valued 1-form on Σ 3 with poles along the marked points. The operator-valued numerators, are then related by a linear identity due to the global residue theorem on Σ 3 , valid in K, i.e. the sum of operators on the left vanishes precisely when it acts on translationally invariant functions. In C, that is, before taking the quotient, the statement of color-kinematics duality is for some operator O µ ∈ C.
As an example, we illustrate (5.8) for the four-point Yang-Mills amplitude. The operator-valued numerators are easily read off from their momentum space counterparts. We have for the s-channel. Other channels are obtained from permutations of this: N t by exchanging 2 ↔ 3 with an overall minus sign, and N u by the permutation (2, 3, 4) → (4, 2, 3). Making use of q 2 i = i · q i = 0, these three numerators are easily shown to obey As claimed, the kinematic Jacobi identity is violated only up to something proportional to the diagonal translation generators. With this we have shown that numerators given by Res vγ (ϕ ± ) satisfy the celestial color-kinematics duality and (5.5) makes manifest the the celestial double copy structure for any number of external particles extending the results of [17].

Discussion
Ambitwistor strings are important tools for the study of flat space holography as they provide a framework to study celestial amplitudes to any multiplicity and in several dimensions. This is exemplified by our all-multiplicity expressions (3.10), (3.12) and (3.15) which allowed us to generalize the computations of [18] to general dimensions and polarizations. In fact, we showed that for multiplicities with n > D, amplitudes in D dimensions can be uniformly expressed in terms of Gelfand A-hypergeometric functions. For n ≤ D, we found distributional expressions with (D − n + 1)-dimensional delta functions coming from residual momentum conservation. We also observed that the Mellin transforms of biadjoint scalar amplitudes are marginally convergent in D = 6, just as the Yang-Mills ones in D = 4 [4]. These happen to be the dimensions in which the classical theories are conformal.
There are also worldsheet models specialized to D = 4 with target space I that have been shown to compute Yang-Mills and gravity amplitudes in the plane wave basis [75,76].
These models are adapted to the spinor-helicity formalism and might be better adapted to the study of a conjectural 2d CFT on I . It would be interesting to understand how the computational methods we used above can be adapted to this case, and to compute the celestial amplitudes in these D = 4 models.
Beyond explicit expressions for the n-point celestial amplitudes, our worldsheet formulas have provided a new outlook on their double copy. Celestial color-kinematics duality and double copy depend on two properties of the amplitudes: kinematical numerators can be represented as operators acting on external kinematics, that is, acting only on quantities at the boundary of spacetime; and that total derivatives can also be characterized as operators acting only on the external kinematics. The latter is simply the statement that amplitudes are translation invariant in the celestial case even if not manifestly so.
These properties might be expected to hold only in the simpler case of flat spacetime. But we know of at least one case where they also hold for amplitudes in a curved background, namely, AdS. The works [40,41] used the ambitwistor string to write CHY-like formulas for amplitudes in AdS spacetimes. These have similar structure to the celestial ones, with kinematical numerators and the scattering equations acting as operators on the scalar contact vertex. These amplitudes are not translation invariant. Instead, the decoupling of total derivatives on the scalar contact diagram can be identified with a quotient by the ideal of diagonal conformal transformations in analogy with the celestial case. We then expect that a similar notion of color-kinematics duality holds for AdS with kinematical numerators obeying a relation like (5.8) up to some symmetry of the space of external data.
Double copy in AdS can also be expected to hold in a similar fashion to the celestial double copy. To show this we must first find the appropriate CHY numerators taking into account that some terms which naively look like non-zero contributions might decouple on top of the scattering equations. To characterize such terms a generalization of twisted cohomology analogous to the one we defined for celestial amplitudes would be very useful. We expect that several acceptable numerators could be found using insights from the ambitwistor string together with an interpretation of these numerators as a generalized twisted cohomology.
An interesting question is whether the framework introduced above holds for loop amplitudes. There are a couple of ambitwistor formulas for loop amplitudes [48,50], the most successful being the ones based on nodal surfaces [49,[51][52][53]. The latter can be generalized to the celestial case making use of our replacement rule (3.6) to extract numerators in front of the Mellin transform of a scalar quantity. The outstanding issue is that the loop-level scattering equations and the numerators can depend on the loop momentum which we'd rather not have in a purely celestial description. If one goes back to position space Feynman diagrams this problem is absent since we can leave all the internal propagators in their position space representation. Numerators that don't depend on loop momentum can still be pulled out in front as operators acting on a scalar loop diagram, but now there is a proliferation of different topologies coming from the loop diagrams. It would be interesting to see if there is a way to encode the effect of loop momentum in the numerators in terms of operators acting on the external variables. Perhaps, some insight could be gained from explicit Mellin transforms giving loop-level celestial amplitudes [16,19,77].

A Mellin transform of the contact diagram
In this appendix we evaluate the celestial n-point contact diagram There are two cases to consider: n > D and n ≤ D. In the former case, we will find the structure of Aomoto-Gelfand hypergeometric integrals which also appear in four dimensions [18]. In the latter case, we will be able to perform all the Mellin integrals and the result will be distributional due to leftover delta functions.
a) n > D : We begin by solving the delta functions for ω l , l = 1, 2, . . . , D, in terms of ω r , r = D + 1, . . . , n, and the other variables. To do this, we use Cramer's rule. Define the determinants where ε µ 1 µ 2 ...µ D stands for the D-dimensional Levi Civita symbol. Then first write l s l ω l q µ l = − r s r ω r q µ r . (A.3) Contracting both sides with ε µµ 1 µ 2 ...μ l ...µ D q µ 1 1 q µ 2 2 · · · q µ l l · · · q µ D D (where a hat denotes omission), we can solve for ω l to find The leftover integral has the standard form of an Euler-type integral.