Asymptotic Charges and Coherent States in QCD
Riccardo Gonzo Tristan McLoughlin Diego Medrano Anne Spiering

=1 
             

 
                     


 




























@width






@ignsp#1  #1empty#1@ignsp


















    a]Riccardo Gonzoa], Tristan Mc Loughlinb], Diego Medranoa]and Anne Spiering
    [a]School of Mathematics & Hamilton Mathematics Institute,
    
    	 Trinity College Dublin, Ireland[b]Instituto de Física Teórica UAM/CSIC & Universidad Autónoma de Madrid,
    
    	 C/Nicolás Cabrera 15, E-28049 Madrid, Spain 
    gonzo@maths.tcd.ietristan@maths.tcd.ied.medrano@csic.esspiering@maths.tcd.ie
                 We study the connection between asymptotic symmetries in non-Abelian gauge theories and the generalised coherent states following from the application to QCD of the Faddeev-Kulish approach to asymptotic dynamics. We compute the large gauge transformation properties of the soft evolution operators and use this to define the quantum corrected, non-linear contribution to the asymptotic charges. We then compute the  leading IR-divergent part of the one-loop correction to matrix elements of the charges and show that, with a specific ordering of soft limits, the  asymptotic charges continue to be conserved at this order.
                                  
    Asymptotic Charges and Coherent States in QCD
        [
        December 30, 2023
    =============================================
    
    
    
    § INTRODUCTION
    
    
    The study of the infrared behaviour of gauge theory scattering amplitudes has a long history and in the case of QED has essentially been understood since the work of Bloch and Nordsieck <cit.>, though it has been refined over the years <cit.>. 
    The standard approach involves the computation of amplitudes which are formally singular - they in fact vanish after exponentiation of the perturbative divergences. One then focuses on inclusive quantities involving arbitrary numbers of real soft-photons which cancel the IR-divergences from virtual photons in loops.  An alternative approach is to directly formulate infrared-finite S-matrix elements by choosing appropriate asymptotic states. For QED this approach where the asymptotic states are not eigenstates of the photon number operator but rather have the form of coherent states, was, starting from the work of Chung <cit.>[The idea of removing the IR-divergences from S-matrix elements by using appropriately dressed asymptotic states was independently considered in  <cit.>. We are grateful to Mario Greco for bringing this reference to our attention.], developed by a number of authors, e.g. <cit.>. Faddeev and Kulish, building on the work of Dollard for the Coulomb problem in non-relativistic quantum mechanics <cit.>, related the structure of these coherent states to the form of the large-time Hamiltonian <cit.>. The approach was partially extended to the much more complicated case of non-Abelian gauge theory <cit.>, <cit.> and more recently to perturbative gravity <cit.>. 
    
    
    The observation that asymptotic conservation laws, which follow from
    Noether's second theorem <cit.> for large gauge transformations, pave the way to understanding the infrared dynamics of gauge theories <cit.> (see <cit.> for a review and more complete references), has lead to a renewed interest in the study of coherent state operators and soft-dressing more generally <cit.>.   
    It has already been demonstrated that for QED <cit.> and perturbative gravity <cit.> the   coherent states relevant to the construction of an infrared finite S-matrix follow from the symmetry of asymptotic charges. 
    The existence of an infinite dimensional symmetry group has lead to the interpretation of the QED vacuum as being degenerate and that  scattering processes are accompanied by a shift in the vacua. Infrared divergences due to massless particles which result in the vanishing of S-matrix elements are thus connected with the “wrong" choice of the in- and out-vacua, and a cure can be sought in a systematic way by considering the asymptotic charges. An analogous statement can be made for perturbative gravity, using BMS supertranslation charges to find suitable asymptotic states for an infrared-finite gravity S-matrix. 
    
    
    The corresponding understanding of QCD infrared dynamics is significantly less complete. Compared to QED, the natural complication that arises is that gauge bosons self-interact in a non-trivial way. The persistence of these non-trivial self-interactions at early and late times is central to the behaviour of the infrared properties of QCD<cit.>. It is our aim here to understand the connection between asymptotic symmetries in non-Abelian gauge theories and generalized IR leading coherent states. We will follow the approach of  Catani,  Ciafaloni and  Marchesini,  <cit.> (see <cit.> for an introduction) which uses energy ordering in each interaction to systematically organise the divergences due to soft gluons. Non-Abelian gauge theories of course also have collinear divergences which have been treated in the coherent operator approach <cit.>, see also <cit.> and the recent interesting work <cit.>.  We will focus on the leading IR-divergence, though the method can be extended in principle to include sub-leading divergences. To be slightly more precise, when we extend the tree-level considerations to one-loop we will focus on the IR double pole terms arising from simultaneously soft and collinear regions of parameter space.  
    
    
    In the following a central goal is the definition of asymptotic charges at loop level. As we discuss below we start from the linear part of the charge that arises in the classical construction and which can be understood to act on Fock, or bare, asymptotic states at null-infinity. However, as the S-matrix is ill-defined for such states we must introduce the soft-evolution operators generated by the QCD Hamiltonian which is used to define dressed states. We then consider the transformation properties of these soft-evolution operators under transformations generated by the asymptotic charges and we show how the soft-evolution operators can be viewed as generating the non-linear part of the charge from the linear part. This implies the existence of an infinite sequence of quantum corrections to the non-linear charge.  This charge can then be used to define a Ward identity for scattering amplitudes. 
    The main result of this paper is the computation, at leading IR divergence, of the one-loop corrections to the Ward identity of the asymptotic charge which is related to the soft gluon theorem. We demonstrate explicitly that the contribution from the soft gluon dressing factor cancels the contribution from the one-loop soft factor in QCD (<cit.>), showing that the asymptotic charge produces a vacuum that is orthogonal to all scattering states built on the original vacuum as in <cit.>. We will see that this result depends on the precise prescription for the order of soft limits in the definition of the charge and dressing factor which is analogous to the order of limits discussed in <cit.>. There are two order-of-limit choices we must make, for the first, in the definition of the charge we take a order of limits used in <cit.>, which corresponds to the standard choice in dimensionally-regularized soft limits. This choice is motivated, as in <cit.>, by its use in problems such as the computation of physical cross-sections. For the second choice, in the definition of the dressed states, we compute the result using both orderings for comparison. As we discuss in sec:disc it is however natural to choose the prescription which preserves, where possible, the symmetries of the theory. We will see that there is indeed a particular prescription that preserves the asymptotic symmetries to the order we compute, as we find that with this choice the Ward identity receives no corrections at one-loop and leading IR-divergence. 
    
    
    
    
      
    Preliminaries
    
    Perturbative computations are relevant in QCD as they are related to experimental observables due to two important properties: asymptotic freedom and factorisation. Confinement naturally sets a scale Λ_QCD such that for partons to be well-defined objects we require the existence of a lower cutoff 
    
    
    
    
    
    on the lowest momentum transfer in a given process. Perturbative quantities are then related to physical cross-sections by convolution with non-pertubative but universal objects, e.g. parton distribution functions.
     In this work we will initially consider asymptotic states formed by acting on the Fock vacuum with creation/annihilation operators for the hard incoming or outgoing partons (these can be either massless gluons or massive quarks)
    
    
        ∏_ℓ  b^†_α_ℓ, σ_ℓ(p_ℓ)|0⟩=|{p_ℓ, α_ℓ, σ_ℓ}⟩
    
    which are labelled by momenta p_ℓ,  colour indices α_ℓ (corresponding to the fundamental representation for quarks and the adjoint for gluons) and helicity indices σ_ℓ as appropriate.
    
    We will consider n-particle IR divergent S-matrix elements between such asymptotic states
    
    
        ℳ_n({p_ℓ,α_ℓ,σ_ℓ})=   ⟨ 0| (∏_ℓ∈out  b_α_ℓ,σ_ℓ(p_ℓ))S(∏_ℓ∈in  b^†_α_ℓ,σ_ℓ(p_ℓ)) |0⟩
        
        =   ^n-2∑_L=0^∞^2Lℳ^(L)_n({p_ℓ,α_ℓ,σ_ℓ})
    
    which give the usual perturbative scattering amplitudes. The behaviour of such amplitudes as individual gluons become soft and the relation to the asymptotic charge plays a key role in our work. We will properly introduce and define the linearized asymptotic charge, Q^lin_ϵ,  in subsequent sections however for now it is sufficient to state that it involves soft gluon creation/annihilation operators a^a†_σ(ω_q)/a^a_σ(ω_q) with vanishing gluon energy ω_q≃ 0 and has matrix elements of the form
    
    
        ⟨out | Q^lin_ϵ S |in⟩∼⟨out | lim_ω_q → 0ω_q a^a_σ(ω_q) S |in⟩
    
    with a careful interpretation of the limit.
    
    At tree level it is known that the limit in the definition of the charge can be understood as
    
    
        ⟨out | lim_ω_q → 0ω_q a^a_σ(ω_q) S |in⟩ =   lim_ω_q → 0ω_q ⟨out | a^a_σ(ω_q) S |in⟩
         =   lim_ω_q → 0ω_q ℳ^(0)_n+1({q, a, σ}) .
    
    These terms can be computed using the tree-level soft-gluon theorem 
    
    
        lim_ω_q → 0ℳ^(0)_n+1({q, a,σ},{p_ℓ})= g_ YM  J^(0)a_σ(q) ℳ^(0)_n({p_ℓ})
    
    with the soft-current given by
    
    
        J^(0)a_σ(q)=[ ∑_ℓ∈outp_ℓ·ε^σ (q)/p_ℓ· q  t_ℓ^a -  ∑_ℓ∈inp_ℓ·ε^σ (q)/p_ℓ· q  t_ℓ^a] ,
    
    where the soft gluon of momentum q, colour a and helicity σ is taken to be outgoing. In this expression, and similar expressions below, as the limit does not strictly exist the notation lim_ω_q → 0 should be understood as referring to the leading term in an  expansion in small ω_q. If the gluon was incoming, there would be an overall minus sign and the corresponding polarisation vector would be ε̅^σ. 
    
     At loop level the issue of the soft limit is more subtle: after computing the perturbative terms defining the matrix elements one may attempt to take the soft limit at the level of the integrands before performing loop integrations or alternatively one may keep ω_q finite and take the limit only after performing the loop integrations. It is known, from the case of subleading IR behaviour of graviton amplitudes <cit.>, that the order has important consequences for the interpretation of quantum corrections to the Ward identities for asymptotic symmetries. In this work for the definition of matrix elements of the asymptotic charge  we always take the soft limit after the computation of matrix elements. That is we will define a regularised charge Q^lin_ϵ(ω_q) and then define
    
    
        ⟨out| Q^lin_ϵ(0) S |in⟩ := lim_ω_q → 0⟨out | Q^lin_ϵ(ω_q) S |in⟩
    
    so that for the insertion of the soft gluon operator we use the usual soft gluon theorem. The regularisation of the charge will be related to that of the soft-limits of amplitudes. When we discuss the soft-evolution operators we will follow <cit.> and use an energy cut-off. However for explicit computations we will make use of the one-loop soft-limits which were computed using dimensional regularisation, with parameter ϵ̂=(d-4)/2, in <cit.> (see also <cit.> for an earlier work in the case where the quarks are massless) and can be written as
    
    
        lim_ω_q→ 0   ℳ^(1)_n+1({q, a, σ}, {p_ℓ})=
        g_ YM J^(0)a_σ(q) ℳ^(1)_n( {p_ℓ})
        +g_ YM^3 J^(1)a_σ(q) ℳ^(0)_n({p_ℓ}) ,
    
    where the first term on the right-hand side is the iterated tree result involving the tree-level soft current (<ref>) while the second term is due to the one-loop soft current which is, to leading divergence, 
    
    
        J^(1)a_σ(q)=-C_A/16π^2 ϵ̂^2J^(0)a_σ(q) + 𝒪(1/ϵ̂) ,
    
    where C_A is the adjoint quadratic Casimir. It is important to note that this leading double pole is due to both collinear and soft divergences at leading logarithmic accuracy and as we will see the coherent state will deal with both of them at this order, as explained in <cit.>.
    
    
    However before computing such matrix elements we first understand in more detail the definition of the asymptotic charge and its expression in terms of free-field operators. In order to do this we must review both the classical construction of the charges, the soft-evolution operators in QCD and the transformation properties of the evolution operators under large gauge transformations. 
    
    
    
    § ASYMPTOTIC CHARGES FOR QCD
    
    
    We will be interested in the asymptotic charges of Yang-Mills theory related to large gauge transformations, which are those that are non-vanishing on the boundary of space-time. Quite generally, Noether's second theorem relates the existence of a local symmetry to a two-form κ^νμ which can be integrated over a codimension-two sphere to define a charge. That is, the local symmetry implies that there exists a conserved current
    
    
        j^μ(ϵ) = S^μ(ϵ) + ∂_νκ^νμ(ϵ) ,
    
    where the function ϵ(x) parametrizes local symmetry transformations and the current S^μ(ϵ) vanishes on-shell
         S^μ(ϵ) eom= 0   .
    
    If ϵ is constant, one recovers the usual conserved current that couples to the gauge field. If we consider the integral of j^μ(ϵ) over a manifold Σ
    
        ∫_Σ d Σ_μ j^μ(ϵ)    = ∫_Σ d Σ_μ[S^μ(ϵ) + ∂_νκ^νμ(ϵ)] 
           eom=∫_σ = ∂Σ d σ_μνκ^νμ(ϵ) ,
    
    this charge is non-vanishing only if the function ϵ has support at the manifold boundary ∂Σ, otherwise it is trivially zero.  This is equivalent to the statement that gauge symmetries are not really “symmetries" in the proper sense. 
    
    Let us review the construction of classical non-Abelian asymptotic charges at ℐ^± for QCD in four dimensional Minkowski space using the usual conformal compactification (for details see <cit.> and see Appendix <ref> for the notations used in this paper). The two-form, κ, is given by the 
    field strength associated to the gauge field 𝒜_μ = 𝒜_μ^a T^a
    
        ℱ_μν = ∂_μ𝒜_ν - ∂_ν𝒜_μ -i   [𝒜_μ, 𝒜_ν]
    
    which obeys the equations of motion
    
    
        ∇^μℱ_μν -i [𝒜^μ, ℱ_μν] =   j_ν^M ,
    
    where j_ν^M is the matter colour current. The classical charge for non-Abelian large gauge transformations with parameter ϵ is then given as, 
    
        Q_ϵ = ∫_σ∗tr[ϵℱ] ,
    
    see for example <cit.> and also <cit.>. 
    
    The relevant surfaces for our asymptotic charges are future and past null infinity, usually denoted ℐ^+ and ℐ^- (see fig:Penrose).
    The natural coordinates for discussing ℐ^+ are the so-called retarded Bondi coordinates
    
    
        r^2 =   ∑_i = 1^3 (x_i)^2 ,      u = t - r ,      z = x_1 + i x_2/r + x_3 ,
    
    such that the flat-space metric becomes
    
        d s^2 =    - d u^2 - 2 d u d r + 2 r^2 γ_z z̅ d z d z̅  with  γ_z z̅ = 2/(1 + z z̅)^2 .
    
    The boundary at r→∞, ℐ^+, is thus parametrized by the coordinates (u, z, z̅).
    The advanced coordinates, (r, v=t+r, z, z̅), are most convenient for ℐ^- and where necessary we will apply the antipodal matching conditions as in <cit.>.
    
    Furthermore, as it is convenient for making contact with perturbative computations,  relevant for the construction of the coherent states, we will focus on Lorenz gauge in this paper with the gauge fixing condition
    
    
        ∇^μ𝒜^a_μ = 0 .
    
    The same charge can be found in non-covariant gauges such as retarded radial or radiation gauge. In general the asymptotic symmetries may depend on the choice of 
    gauge fixing, which can also be seen at the path integral level (see <cit.>), but for the leading asymptotic charges one finds the same result for both covariant and physical gauge choices. 
    Finally we will be interested in large gauge transformations parametrised by a function ϵ = ϵ^a(z, z̅) T^a which labels the transformations at infinity, for example at ℐ^+,
    
    
        δ_ϵ A_z(u , z, z̅) = D_z ϵ(z, z̅) ,
    
    where A_z(u , z, z̅)=lim_r→∞𝒜_z(r, u , z, z̅). 
    
    In retarded Bondi coordinates on the celestial sphere 𝒮^2 defined by the u→ -∞ region of ℐ^+, i.e. ℐ^+_-, the charge (<ref>) becomes 
    
    
        Q_ϵ =  ∫_ℐ^+_- d^2 z γ_z z̅ tr[ϵ(z, z̅) F_r u] ,
    
    where F_r u=lim_r→∞ r^2 ℱ_ru. This charge is usually denoted Q^+_ϵ to distinguish it from the corresponding charge Q^-_ϵ on ℐ^-. In order to avoid a proliferation of superscripts we will mostly neglect this index and hopefully it is clear from context to which charge we are referring.
    
     The charge can now be rewritten, assuming the field strengths vanish at ℐ^+_+ (for details see appendix <ref>), as 
    
    
        Q_ϵ = Q_ϵ^lin + Q^non-lin_ϵ ,
    
    where
    
    
        Q^lin_ϵ = ∫_ℐ^+ d^2 z d u ϵ^a(z, z̅) [ ∂_u (∂_z A^a_z̅ + ∂_z̅ A^a_z) ]
    
    and
    
    
        Q^non-lin_ϵ   = ∫_ℐ^+ d^2 z d u   ϵ^a(z, z̅)  [ f^a b c (A^b_z ∂_u A^c_z̅ - A^c_z̅∂_u A^b_z) +   γ_z z̅ j_u^a(2)] .
    
    Here we have adopted the standard splitting of the asymptotic charge into a piece linear in the gauge field (usually called soft) and a piece non-linear in the gauge fields (usually called hard). 
    In the construction of asymptotic charges it is usually argued that it is the linearised theory which is used to define the two-form. This is the case in the framework of covariant charges <cit.> and in a related context was noted by Abbott and Deser <cit.> for Yang-Mills theory. 
    In the current context this would correspond to  linearising the theory around the background A^a_z = A^a_z̅ = 0 and would thus discard the non-linear gluon piece of the asymptotic charge. However,  
    it was argued that one should include the non-gluon piece inside the hard part of the charge <cit.>.  As we will explicitly discuss later, we start from the linear charge (<ref>) acting on asymptotic states however such a hard contribution in our case is indeed present as it will be generated by the time evolution given by the Mø ller operator. At ℐ^- there is an analogous construction where the linearised asymptotic charge is given by
    
    
        Q^lin_ϵ = ∫_ℐ^- d^2 z d v  ϵ^a(z, z̅) [ ∂_v (∂_z A^a_z̅ + ∂_z̅ A^a_z) ]
    
    in advanced Bondi coordinates.
    
    We will be interested in understanding the action of these charges on scattering states and in order to hew closely to the standard S-matrix formulation we will start with the usual equal-time commutators on a space-like hypersurface and by taking the t →±∞ limit of the usual free field expansion
    
        𝒜_μ^a(x)    =  ∑_σ = ±∫dq[ε̅_μ^σ(q) a^a_σ(q) e^iq· x + ε_μ^σ(q) a^† a_σ(q) e^-iq· x] .
    
    The sum is over helicities with polarisation vectors ε^σ_μ and the equal time commutation relations are given in Appendix <ref>.
    We can then evaluate the fields arbitrarily close to ℐ^+ by using the saddle point approximation in the r→∞ limit <cit.>  such that, 
    
    
        A_z^a   = -i/8 π^2√(2)/(1 + z z̅)∫_0^+∞ d ω_q  [a^a_+(ω_q x̂) e^-i ω_q u - a^a, †_-(ω_q x̂) e^i ω_q u] 
        
        A_z̅^a   = -i/8 π^2√(2)/(1 + z z̅)∫_0^+∞ d ω_q [a^a_-(ω_q x̂) e^-i ω_q u - a^a, †_+(ω_q x̂) e^i ω_q u] ,
    
    where q^0=ω_q and x̂=x⃗/r. 
    The asymptotic charges can now be expressed in terms of the creation and annihilation operators and for the linear, or soft, charge we have
    
    Q^lin_ϵ   =   lim_ω_q → 0ω_q/4 π∫ d^2 z √(γ_z z̅)[∂_z ϵ^a(z, z̅) (a^a_-(ω_q x̂) + a^a†_+(ω_q x̂)) 
    +h.c.] 
    
    which provides the connection with the action discussed in the introduction eq.(<ref>).
    The non-linear part of the charge is more complicated to interpret in the quantum theory as it requires normal ordering. For the gluon part we have
    
        Q^non-lin_ϵ|_gluon=    i f^a b c∫_ℐ^+d^2 z  /2(2π)^3γ_z z̅ϵ^a(z, z̅)  {∫_0^+∞ d ω_q ω_q ∑_σ=± a^b†_σ(ω_q x̂) a^c_σ(ω_q x̂) 
           -45pt  +  1/2∫_0^+∞∏_i = 1^2 d ω_q_iδ( ω_q_1 + ω_q_2) ( ω_q_1-ω_q_2)  [ a^b_+(ω_q_1x̂) a^c_-(ω_q_2x̂) - a^b†_-(ω_q_1x̂) a^c†_+(ω_q_2x̂) ]}
    
    where we have normal ordered the expression and dropped the resulting constant.  As the energy integral is over only positive values, the delta-function in the second term has no non-zero support in the integration region and it is therefore vanishing. Defining the gluon number density operator 
    ρ^a_g(q⃗)=-i f^abc∑_σ=± a^b†_σ(ω_q x̂) a^c_σ(ω_q x̂)
    
    we can write the first term in the particularly simple form
    
    Q^non-lin_ϵ|_ a^† a-gluon=   -∫dq ρ_g^a(q⃗) ϵ^a(x̂) .
    
    
    There is one important subtlety in relating the linear charge to large gauge transformations. Using the equal-time commutators and the expression (<ref>) one can show that
           [Q_ϵ^lin, a^a_+(q⃗)] =-(2π)^2δ(ω_q)√(γ_zz̅)∂_z ϵ^a(z, z̅)
          [Q_ϵ^lin, a^a †_+(q⃗)] =(2π)^2δ(ω_q)√(γ_zz̅)∂_z̅ϵ^a(z, z̅) 
           [Q_ϵ^lin, a^a_-(q⃗)] =-(2π)^2δ(ω_q)√(γ_zz̅)∂_z̅ϵ^a(z, z̅) 
           [Q_ϵ^lin, a^a †_-(q⃗)] =(2π)^2δ(ω_q)√(γ_zz̅)∂_ zϵ^a(z, z̅) .
    
    Thus the gauge transformation of the asymptotic field  A_z^a(u, z, z̅) is given by
           [Q^lin_ϵ, A_z^a(u, z, z̅)]=i2∂_z ϵ^a( z, z̅))  ∫_0^∞ dω_q  δ(ω_q)(e^-i ω_q u+e^+i ω_q u) .
    
    This will  give the incorrect large gauge transformations using the usual definition of the delta-function but it can be remedied by inserting a factor of two for the zero-mode contribution. This is related to the fact noted in the context of the Abelian theory <cit.> that
    the radiative phase space at ℐ^+ defined as Γ_+ = {F^a_u z, F^a_u z̅} will not give the usual linearised large gauge transformation. It is easy to check that
    
    
        [Q^lin_ϵ, A_z^a(u,z,z̅)] = i/2∂_z ϵ^a(z, z̅) ≠ i δ_ϵ A_z^a(u,z, z̅)
    
    using the standard Poisson brackets of the non-Abelian gauge theory
    
    
        [A_z^a(u,z,z̅), A_w^b(u^',w,w̅)] = - i /4δ^a bΘ(u - u^') δ^2(w -z) .
    
    The problem is that the zero modes at the boundary correspond to a single real scalar field and not to a complex one: in the zero frequency limit the two helicities are identified and there is a miscounting of the degrees of freedom. This problem can be solved as in <cit.> by imposing additional constraints at the boundaries of ℐ^+.  Alternatively one can add the factor of two for zero modes as was done by the authors of <cit.> used in the definition of the QED coherent state operator. We will follow a similar procedure however in our case we modify the charge, adding a factor of two for zero-modes, and leave the coherent state, and the amplitudes, unchanged. 
    
    
    
      
    Splitting of the charge contributions according to the energy scale E
    
    A key point is that once we pick a scale E we split the Fock space of free particles into the hard and the soft parts ℋ= ℋ^E_s ⊗ℋ^E_h. Indeed we could make a splitting of the energy integral to separate the two contributions
    
    
        ∫_λ^+∞ d ω_q →∫_λ^E d ω_q + ∫_E^+∞ d ω_q  ,
    
    where we have introduced λ as an infrared cutoff.
    As mentioned previously, we will in fact take the linearised charge in the asymptotic region as our starting point and see how the non-linearity emerges from commuting with the  evolution operators. As we will see this corresponds to effectively having
    
    
        Q^lin_ϵ   = Q^lin_ϵ,s⊗𝕀_h     and   
        Q^non-lin_ϵ = 𝕀_s ⊗ Q^non-lin_ϵ,h|_gluon+quarks ,
    
    where we have a non-linear hard term for both gluons and quarks. This is similar to the approach taken in <cit.>.
    
    
    
    § ASYMPTOTIC HAMILTONIAN AND SOFT EVOLUTION OPERATORS
    
    
    The starting point for the coherent state approach to IR divergences <cit.> is the choice of an appropriate asymptotic Hamiltonian describing the parton dynamics in the far future and far past. We will review,  following the arguments of <cit.>, how one can carry out the Faddeev-Kulish approach in the non-Abelian case at the leading order in the IR divergences. We start from the splitting of the standard QCD interaction Hamiltonian in the interaction representation into soft and hard parts:
    
    
        H^I (t) = H_h^E(t) + H_s^E(t) .
    
    This is done by introducing  at each interaction vertex the energy transfer ν
    
        ν = | ∑_iη_i ω_i |     ∑_iη_i q_i = 0 ,
    
    where ω_i are the energies of the interacting particles with η_i = +1 (resp. -1) for incoming (resp. outgoing) particles. We define the soft part of the Hamiltonian as containing only energies below a scale ν< E and we also introduce a lower cut-off λ < ν. The lower energy bound λ is not only required by our use of perturbative QCD but also to have a good definition of the FK states <cit.> however in the following we will sometimes leave it implicit. It is important to note that the region ν<E contains both soft and collinear subregions. For example <cit.> considered a cubic vertex with incoming gluon energy ω_1=|𝐪_1| and outgoing energies ω_i=|𝐪_i| for i=2,3. In particular ω_3=√(ω_1^2+ω_2^2 -2 ω_1ω_2 cosθ_12), where θ_12 is the angle between 𝐪_1 and 𝐪_2. The condition (<ref>) then defines a hyperbola in ω_2–cosθ_12 plane 
    
        cosθ_12 =-(1+E/2ω_1)E/ω_2+(1+E/ω_1 )
    
    which defines regions of small ω_2 and cosθ_12∼ 1. Alternatively one can define the asymptotic region by imposing both an angular cutoff θ_12<Θ and an energy cutoff ω_2<M. In this way either form of the cutoffs can be used to treat the soft and collinear regions for any of the particles. 
    
    The soft Mø ller operators are defined as the standard time-ordered product
    
    
        Ω^E_± = T exp[-i ∫_∓∞^0 H_s^E(t) d t ] ,
    
    but in order to isolate the leading IR singularities, which are of the form (^2log(E/λ))^k at k-loops, it is useful to transform to frequency space and so write the Hamiltonian as <cit.>
    
        H_s^E(t) = ∑_η = ±∫_λ^E d ν  h^η(ν) e^-i ην t .
    
    Using this expression we can write the soft-evolution operator as 
    
    
        Ω^E_± =∑_n=1^∞   ∑_η_i=±∫_λ^E dν_1 … dν_n h^η_n(ν_n)… h^η_1(ν_1)/(η_n ν_n+… +η_1 ν_1± i 0 )… (η_1 ν_1 ± i 0) .
    
    Infrared singularities come from vanishing energy denominators, and in particular leading logarithms come from the region specified by the strong ordering
    
    
        λ≪ν_1 ≪ν_2 ≪…ν_n ≪ E ,
    
    where we can thus approximate
    
    
        1/η_n ν_n + … + η_1 ν_1 ± i 0· .... ·1/η_1 ν_1 ± i 0≃∏_i = 1^n 1/η_i ν_iΘ(ν_n > …  > ν_1) .
    
    To leading order (double pole accuracy) in the IR divergences, for which we can neglect the distinction between Ω^E_+ and Ω^E_- due to the i 0 prescription at leading order, the Mø ller operator is then given by the frequency-ordered exponential
    
    
        Ω^E    = 𝒫_νexp(∫_λ^E  ∑_η = ±d ν/η ν h^η(ν)) .
    
    Furthermore, it is possible to show that at this order in the IR:
    
      *   In each three-gluon vertex we may assume there is always one gluon which is much softer than the others and so we can use the eikonal approximation,
    	
      *  Quarks interact only via eikonal vertices and pair production is neglected because the process g → q q̅ is IR finite,
    	
      *  Four-gluon vertices and ghost contributions can be neglected.
    
    With these assumptions it is possible to write the soft interaction Hamiltonian as a sum of two terms, one depending on quarks and a soft gluon and a second purely gluonic cubic term
    
    
        H_s^E(t) = H_ffg^E(t) + H_ggg^E(t) .
    
    These can be combined so that
    
    
        H_s^E (t) = -g_ YM∫_ω_qd p∫^E_λd q ρ^a (p) p̂· [a^a(q) e^i p̂· q t + h.c.] 
              with   ρ^a(p) =ρ_f^a(p)+∑_σ a^† b_σ(p) T_A^a_bc a^c_σ(p) ,
    
    where ρ now contains a contribution not only from fermionic matter ρ_f but also from the hard gluons with energies ω_p> ω_q.  It is important to note that while the density operator involves gluons which are harder than the soft gluon at that vertex it does not commute with all soft gluon operators. This is fundamentally different from the QED case and corresponds to the non-linear nature of the gauge symmetry.
    
    
    Using the Hamiltonian H_s^E(t) in (<ref>) the soft Mø ller operator becomes
    Ω^E
       =   𝒫_ωexp[ ∫_λ^E dq 𝒥_q·Π_q] ,
    
    where  Π^a_μ(q)=a_μ^a(q)-a^a †_μ(q) is the displacement operator and
    𝒥^a_q_μ= ∫_ω_qdp ρ^a(p)p_μ/p· q
    and the dot product denotes the contraction of both Lorentz and gauge indices where appropriate. The exponential is interpreted as being ordered in the soft gluon energies with smaller energies to the right. Here we have switched from the frequency ordering to an energy ordering. For a complete equivalence to the frequency ordering we should also include an ordering on angles and a cutoff on the angular region. However as we are only interested in the leading IR region - which is both soft and collinear- it is sufficient to only consider the energy cutoff. Where this leads to expressions appearing singular an angular cutoff can be mentally added. For our explicit computations in subsequent sections we in fact use dimensional regularisation which simultaneous treats both soft and collinear singularities but the frequency/energy ordering is useful for the physical picture it provides. 
    
    Finally, the soft evolution operator can be used to define an IR-finite S-matrix S^E by removing the IR singularities due to initial and final state interactions
    
    S^E=Ω^E_-  S Ω^E†_+
    
    as shown schematically in fig:SmatrixRel.
    
    
    
    § LARGE GAUGE TRANSFORMATIONS FOR SOFT EVOLUTION OPERATORS
    
    We are interested in understanding the transformation properties of the S-matrix under large gauge transformations and the corresponding Ward identity for amplitudes. As the S-matrix relates states in the far future and far past, where the theory is taken to be free, we will use the linearised charge, (<ref>), when we compute matrix elements of the commutator 
    
        [Q_ϵ^lin, S] := Q_ϵ^+,lin S - S Q_ϵ^-,lin ,
    
    where we use the appropriate linearised charge Q^± for incoming and outgoing states. The key to our approach is to use the soft-evolution operator to relate the free theory to the interacting theory and so relate the linearised charge to the non-linear contributions. As the soft-evolution operator can be defined in the quantum theory, this gives a method of defining the correct quantum non-linear corrections to the charge. This can be done by analysing the gauge transformation properties of the soft-evolution operator. 
    
    Before carrying out the computation, we can see how such non-linear terms imply a Ward identity for the S-matrix: If we had that the evolution-operator transformation was given by a relation of the form
    
    [Q_ϵ^lin, Ω^E]=Q_ϵ^hΩ^E
    
    for some non-linear charge Q_ϵ^h, then we would have the relation
    
        ⟨⟨{p_f, α_f} [Q_ϵ^lin, S] {p_i, α_i}⟩⟩=⟨{p_f, α_f}|[Q_ϵ^lin - Q_ϵ^h, S^E]|{p_i, α_i}⟩ .
    
    Here we have used dressed states, e.g.
    {p_i, α_i}⟩⟩=Ω^E†∏_i∈in b^†_α_i(p_i)|0⟩ ,
    
    in computing matrix elements. 
    If the linearised charge produces states orthogonal to scattering states constructed on the original vacuum, this becomes 
    
        ⟨⟨{p_f, α_f} [Q_ϵ^lin, S] {p_i, α_i}⟩⟩=-
        ⟨{p_f, α_f}|[Q_ϵ^h,S^E]|{p_i, α_i}⟩ .
    
    This relation, graphically represented in fig:TreeQ, is what we refer to as the Ward identity. 
    
    
    
    Our goal is to study this relation in the context of QCD using the soft-evolution operators and capturing the leading IR singularities. To this end we consider the explicit expression for the linearised large gauge transformations using the definition of the linearised charge (<ref>) with the  commutators (<ref>).
    It is convenient to first integrate by parts for the z and z̅ variables and start with the expression 
    
    Q_ϵ^lin=- ∫_ℐ^+ d^2z du [∂_z ϵ^a ∂_u  A^a_z̅+∂_z̅ϵ^a ∂_u  A^a_z ] .
    
     In principle we must carefully account for any boundary terms that occur in this step, however at the end of the computation  we will undo this integration by parts and so remove the boundary terms again. 
    Using the definition of the energy ordering, we  calculate
    
        [Q_ϵ^lin, Ω^E]=∑_n [Q_ϵ^lin, ∫_ω_n-1dq_n…∫_λ^E dq_1 𝒥_n·Π_n…𝒥_1·Π_1] .
    
    As the charge only involves zero modes, it acts only on the last term in each element of the sum, i.e. ∫_λ^E dq_1 𝒥_1·Π_1, and then only for ω_1=λ∼ 0. That is,  we are assuming that the charge acts on the mode with energy ω_1=λ∼ 0  but not on ω_2> ω_1. By straightforward computation one finds
    
    [Q_ϵ^lin, ∫_λ^E dq 𝒥·Π]
       =   𝒩(ϵ) ,
    
    where 𝒩(ϵ)=∫_λdp ρ^a(p) N^a(ϵ, p)
    with
    
        N^a(ϵ,p)=-(2π)^2 ∫dq/√(γ_z z̅)δ(ω_q)
        2p/p· q·(ε^-(q)∂_z ϵ^a(q)+ε^+(q) ∂_z̅ϵ^a(q)) .
    
    This expression for N^a depends on having performed an integration by parts and so ignores any potential boundary terms. For example it can be seen that such a term vanishes for the case ϵ^a being a constant. One can undo the integration by parts, and so recover the dropped boundary terms, 
    and write 
    
        N^a(ϵ,p)=/2π∫ dω_q δ(ω_q) ∫ d^2z γ_zz̅ -p^2/(p·q̂)^2ϵ^a(z,z̅) ,
    
    where we use the notation q̂^μ=q^μ/ω_q and which for the case of constant ϵ becomes N^a(ϵ,p)= ϵ^a ∫ dωδ(ω). 
    In these expressions we must interpret the δ-functions as taking the soft limit in the appropriate fashion. 
    
    We have thus found that 
    
    [Q_ϵ^lin, Ω^E]
       =   Ω^E 𝒩(ϵ) ,
    
    which is almost what we would have expected from the classical analysis of the non-linear charge. However as
    
    [𝒥_j ·Π_j, ∫_λdp ρ^b(p) N^b(ϵ, p) ]=0
    
    only if N^b(ϵ, p)=N^b(ϵ, q_j), we must correctly order the terms to put the expression in the form (<ref>). Performing the non-trivial commutators to one-loop, 𝒪(^2), we find 
    
    [Q_ϵ^lin, Ω^E]= [𝒩(ϵ)+𝒜(ϵ)]Ω^E ,
    
    where
    
        𝒜(ϵ)=
        i g_ YM f^abc∫_λ^E dq∫_ωdp ρ^a(p)
        p·Π^c(q)/p· q·[N^b(ϵ, p)-
        N^b(ϵ, q)] .
    
    Here we see that there are corrections to the tree-level expression for Q_ϵ^h which arise due to the non-Abelian nature of the theory. Indeed we can see that there will be further non-linear corrections at each higher loop orders. These corrections do vanish in the case of a constant ϵ, which is to say that the total colour charge does not receive any corrections. Moreover the correction involves a soft-gluon operator which we might expect to have vanishing contribution when computing matrix elements of the IR-finite S-matrix. We now turn to the computation of exactly such matrix elements. 
    
    
    
    
    
    § WARD IDENTITIES FOR DRESSED S-MATRIX
    
    To  compute matrix elements between dressed states, we make use of the non-trivial fact that, to leading order in the soft divergence, the dressing of external states factorises in colour space  <cit.>
        {p_i,α_i}⟩⟩≡Ω^E† |{p_i,α_i}⟩=∏_i∈in𝒰^p_iE_α_iβ_i(Π) b^†_β_i(p_i) |0⟩ ,
    
    where the coherent-state operator 𝒰^p_iE_α_iβ_i(Π) is a functional of the soft gluons only. The coherent-state operator for a parton in the gauge group representation with generators t_αβ^a is defined by the energy-ordered integral
    𝒰_αβ^pE=P_ωexp[-g_ YM∫^E_λdq p·Π^a_ω(q)/p· q t^a]_αβ ,
    
    where the dressed gluon field is similarly defined by 
    Π^a_ω(q)=𝒰_ab^qEΠ^b(q)
    
    and 𝒰_ab^pE is the adjoint coherent-state operator.  
    
    These non-linear equations can be solved iteratively so that to 𝒪(^2) we have
    
        𝒰^pE_αβ=    δ_αβ - 
        ∫_λdq p·Π^e(q)/p· q t^e_αβ
           +^2 ∫_λdq_1∫_ω_1dq_2
        (
        p·Π^e_2(q_2)/p· q_2 t^e_2_αγ)
        (
        p·Π^e_1(q_1)/p· q_1 t^e_1_γβ)
           -^2  ∫_λdq_1∫_ω_1dq_2
        (
        p·Π^e_2(q_2)/p· q_2)
        (
        q_2·Π^e_1(q_1)/q_2· q_1) · [t^e_2,t^e_1]_αβ .
    
    This dressing factor captures the leading-order effects of soft-gluon radiation of each of the hard partons. It includes all-order effects arising from arbitrary numbers of gluons being radiated, as well as loop effects which arise from normal ordering each of the terms.  
    
    
    
      
    Tree-level Ward Identity While the calculation is essentially identical to the Abelian case, to fix conventions we will start with the tree-level result. The relevant terms for computing matrix elements of the commutator of the linearised charge with the S-matrix are:
    
        ⟨⟨{p_f,α_f} [Q^lin_ϵ,S]  {p_i,α_i}⟩⟩= ⟨ 0|∏_f∈out  b_α_f(p_f)([Q_ϵ^lin,  S]  
           - g_ YM∫_λdq[∑_ℓ∈outp_ℓ·Π^e(q) /p_ℓ· q t^e_ℓ  [Q_ϵ^lin, S] + ∑_ℓ∈in[Q_ϵ^lin, S]p_ℓ·Π^e(q) /p_ℓ· qt^e_ℓ])∏_i∈in b^†_α_i(p_i) |0⟩ ,
    
    where the subscripts on the colour generators indicate the parton leg upon which they act and for convenience we denote, for example, t^a_α_jβ_j b^†_β_j =t^a_jb^†_α_j.
    
    
    The two contributions to the matrix element can be represented graphically as in fig:TreeSoft. The first term, corresponding to diagram (4.a), corresponds to the usual absorption or emission of the soft gluon in the linearised charge. We compute these terms by using the soft-gluon theorem, see (<ref>) and (<ref>),
    and due to the sign difference in the soft limits of incoming and outgoing gluons the two terms of the commutator add rather than cancel. 
    The next terms, corresponding to diagram (4.b) in fig:TreeSoft, arise from contracting oscillators in the linearised charge with those in the coherent states. These terms cancel due to the sign from the commutator. Thus we find 
    
        ⟨⟨{p_f,α_f}  [Q_ϵ^lin, S] {p_i,α_i}⟩⟩= -[∑_ℓ∈out Q_ϵ^h(p_ℓ)-∑_ℓ∈in Q_ϵ^h(p_ℓ)] ℳ^(0)_n ,
    
    where we have introduced the eigenvalue of the non-linear charge Q_ϵ^h(p)=Q^h,a_ϵ(p) t^a where
    
        Q^h,a_ϵ(p)
        =   -8π^2 ∫dqδ(ω) /√(γ_zz̅)·[
        ∂_z ϵ^a(q̂) ε^-· p /q· p
        +∂_z̅ϵ^a(q̂) ε^+· p/q· p]
    
    which, as expected,  is N^a(ϵ, p) as defined in (<ref>).
    
    As in the QED case the charge acts on the Fock vacuum to produce a state orthogonal to all scattering states built on the original vacuum.
    This can be seen by keeping only those terms where the charge corresponds to an emitted gluon, and for convenience considering only incoming hard partons. While we now only have half the terms from the commutator, we can still have contributions of the form represented in fig:TreeSoft, however due to a cancellation between diagrams (4.a) and (4.b), one finds
    ⟨⟨0 Q_ϵ^linS  {p_i,α_i}⟩⟩   =   0 
    at tree-level. 
    
    
    
      
    Finite One-loop S-matrixNow let us proceed to a one-loop calculation, and again for now considering only in-particles, we compute the one-loop matrix elements
    
        ⟨⟨ 0  S {p_i,α_i}⟩⟩    =⟨ 0|S ∏_i 𝒰^p_iE_α_iβ_i b^†_β_i(p_i) |0⟩
            = ⟨ 0| S∏_i  b^†_α_i(p_i) |0⟩
             - g_ YM∑_j t^e_1_j ∫_λdq⟨ 0| S p_j ·Π^e_1(q) p_j· q∏_i b^†_α_i(p_i) |0⟩
             +g_ YM^2/2∑_j≠ k  t^e_1_j t^e_2_k
        ∫_λdq_1∫_λdq_2 ⟨ 0| S
        p_j·Π^e_1(q_1)p_j· q_1p_k·Π^e_2(q_2)p_k· q_2∏_i b^†_α_i(p_i)|0⟩
             + one-parton terms .
    
    The first term is the usual IR-divergent S-matrix  element which arises from diagrams such as (5.a) in fig:IRdvg2pt, while the remaining terms are the compensating IR-divergent terms from the coherent state which are graphically represented in diagrams (5.b) and (5.c). Here we will focus on those diagrams which involve two external partons and neglect one-parton contributions such as those in fig:IRdvg1pt. These one-parton contributions are subleading in the IR divergences and so are not needed, but can in fact be reconstructed by insisting on gauge invariance. The second term (5.b) gives an 𝒪(g_ YM^2) contribution by using the tree-level soft limit of the S-matrix. 
    In particular it can be written as 
    
        I^(5.b) = g_ YM^2 ∑_j≠ k   t^e_j t^e_k ∫dq ∑_σ p_j·ε^σ(q) /p_j· qp_k ·ε̅^σ (q)/p_k· qℳ^(0)_n ,
    
    where we have again dropped the one-parton, i.e. j= k, terms.
    The next term, corresponding to (5.c), involves contracting two gluons from the dressing factors of two partons. We thus need
    
        Π^e_j_μ(q_1)Π^e_k_ν(q_2) =    (a_μ^e_j(q_1)-a^e_j†_μ(q_1))(a_ν^e_k(q_2)-a^e_k†_ν(q_2))
        
        =    -δ(q_1-q_2) δ^e_j e_k∑_σε̅^σ_με_ν^σ +terms with two oscillators
    
    so that we find
    
        I^(5.c) = -g_ YM^2/2∑_j≠ k t^e_jt^e_k ∫dq∑_σp_j·ε̅^σ(q) p_k·ε^σ(q)/p_j· q  p_k· qℳ^(0)_n
    
    which is the same as the (5.b) contribution up to an overall factor. 
    We can rewrite the product of polarisation vectors using (<ref>) which,
    if we include the one-parton contributions and impose total colour conservation to remove the c_μ dependent terms, becomes 
    ∑_σε̅^σ_με^σ_ν=η_μν. 
    Hence we find the S-matrix elements
    ⟨⟨ 0  S {p_i,α_i}⟩⟩   =   [1+ g_ YM^2/2∑_j, k  t^e_jt^e_k  I_jk] ℳ^(0)_n
    
    in terms of the IR-divergent loop integral
    
    I_jk=∫dqp_j· p_k/p_j· q  p_k· q .
    
    While we have previously considered an energy cut-off to regularise divergences, in order to compare with known results it is useful to instead use dimensional regularisation. Keeping the leading divergence in the parameter ϵ̂=d2-2, the loop integral is given by 
    
    I_jk=-1/2(2π)^2ϵ̂^2 ,
    
    which encodes the double-pole singularity due to both soft and collinear divergences at leading logarithmic accuracy. To one-loop order and to our accuracy the amplitude is known to be, <cit.>, 
    ℳ_n=[1+^2/16π^2ϵ̂^2∑_j≠ k  t^e_jt^e_k]ℳ_n^(0)
    and hence we see that, as expected, the singular parts cancel in the S-matrix elements. This can be repeated for generic outgoing states and again one finds that the one-loop leading singularities cancel 
    ⟨⟨{p_f,α_f}  S  {p_i,α_i}⟩⟩ |_𝒪(^2)=0+𝒪(1/ϵ̂)
    
    as expected. An analogous approach to removing the singularities would be the Wilson-line dressing which reproduces the same leading IR-divergences <cit.> (see also <cit.>).
    
    
    
      
    One-loop Ward Identity 
    We now want to compute the one-loop correction to the Ward identity using the dressed states.
    
    We thus consider the terms
    
        ⟨⟨ 0  [Q_ϵ^lin,S] {p_i,α_i}⟩⟩ |_𝒪(g_ YM^3)=⟨ 0| [Q_ϵ^lin,S] {1
        - g_ YM∑_j t^e_1_j ∫_λdq p_j ·Π^e_1(q) p_j· q
           +g_ YM^2/2∑_j_s∏_s=1^2 t^e_s_j_s∫_λdq_s 
        p_j_s·Π^e_s(q_s)p_j_s· q_s +g_ YM^2∑_j [t^e_1_j, t^e_2_j]
        ∫_λdq_1∫_ω_1dq_2 q_2·Π^e_1(q_2)q_2· q_1p_j·Π^e_2(q_2)p_j· q_2
            -g_ YM^3/3!∑_j_st^e_1_j_1 t^e_2_j_2t^e_3_j_3∏_s=1^3  ∫_λdq_s 
        p_j_s·Π^e_s(q_s)p_j_s· q_s
            -g_ YM^3∑_j_1,j_2t^e_1_j_1 [t^e_2_j_2, t^e_3_j_2]
        ∏_s=1^2 ∫_λdq_s p_j_s·Π^e_s(q_s)p_j_s· q_s∫_ω_2dq_3
        q_2·Π^e_3(q_3)q_2· q_3
            -g_ YM^3∑_j_1,j_2t^e_1_j_1 t^e_3_j_2 t^e_2_j_2∏_s=1^2 ∫_λdq_s p_j_s·Π^e_s(q_s)p_j_s· q_s∫_ω_2dq_3
        p_j_2·Π^e_3(q_3)p_j_2· q_3 +one-parton terms}|{p_i,α_i}⟩ ,
    
    where as before the subscripts on the generators denote the hard parton upon which they act. This expression simplifies significantly as
    the contributions where the gluon in the charge is contracted with a gluon in the dressing of the partons cancel between the two terms in the commutator. Thus we only need to keep the contributions where the charge contracts with the S-matrix as these terms add. The non-vanishing contributions are graphically represented in fig:OneLoopSoft where it can be seen that these contributions closely parallel those of fig:IRdvg2pt.  The result can be straightforwardly computed using the previous tree-level and one-loop results. For example diagram (7.b) gives
    
    I^(7.b) = g_ YM^2 ∑_j, k, rQ^h,a_ϵ(p_j) I_kr t^a_j t^e_k t^e_r ℳ_n^(0) ,
    
    where Q_ϵ^h,a is defined in (<ref>) and I_kr is the integral (<ref>). The diagram (7.c) gives the same contribution but with a coefficient of -12, so that I^(7.c) = -1/2 I^(7.b). These terms cancel against contributions such as those in diagram (7.a), which can be computed using the soft limit of one-loop amplitudes given in (<ref>).
    As can be seen in, for example, (<ref>) there are contributions with j ≠ k ≠ r that involve three external hard-parton legs. These contributions straightforwardly cancel against the appropriate terms in the first term of the one-loop soft limit involving the one-loop amplitude, that is the iterated tree term, 
    
        ∑_j Q^h,a_ϵ(p_j) t^a_j ℳ^(1)_n≃g_ YM^2/16π^2ϵ̂^2∑_j≠ k≠ r Q^h, a_ϵ(p_j) t^a_j t^e_k t^e_r ℳ^(0)_n
    
    using (<ref>) for the leading singular terms in the amplitude.
    
    
    Slightly more subtle are the contributions which involve only two external legs. 
    In particular for the diagram (7.b) there is an ambiguity regarding the order in which one takes the soft limits corresponding to the gluon in the charge and the lowest-energy gluon in the parton dressing. The two orderings are shown schematically in fig:OneLoopDSoft. Diagram (8.b.1) corresponds to first taking the soft-limit for the charge before computing the loop integral involving the soft gluon in the dressing factor.  In diagram (8.b.2) we instead take the soft limit for the dressing-factor gluon and then for the gluon in the charge. In particular this gives an ordering of the colour generators 
    
    I^(8.b.2) = - ^2 ∑_j≠ k Q^h,a_ϵ(p_j) I_kr t^e_j t^a_j  t^e_k ℳ_n^(0)
    which combines with the other two-parton contributions from (9.c) to cancel against the two-parton contributions from (7.a) which come from the soft limit of the one-loop amplitude.
    
    The final result is that using the ordering (8.b.2)  we find that the commutator vanishes at one-loop to leading singularity
    ⟨⟨ 0  [Q_ϵ^lin,S]  {p_i,α_i}⟩⟩ |_𝒪(g_ YM^3)   =   0+𝒪(1ϵ̂) .
    
    There are in principle contributions which only involve a single external parton momentum which we have not included as they are subleading in the IR expansion, however these can again be found by using gauge invariance. Thus we find that, using this prescription, there are no corrections to the Ward identity. 
    
    If we had chosen the alternative prescription (8.b.1) we would find a correction which is related to the one-loop soft current and which can be written as 
    
        ⟨⟨ 0  [Q_ϵ^lin,S]  {p_i,α_i}⟩⟩ |_𝒪(g_ YM^3)=-
        ^2 C_A/16π^2ϵ̂^2∑_ℓ∈in Q_ϵ^h(p_ℓ) ℳ^(0)_n .
    
    In this case we see that there appears to be a one-loop effect which in principle could be interpreted as a correction to the hard charge. Related to this we must also establish that the linearised soft charge acting on the vacuum produces states that are orthogonal to the dressed scattering states at one-loop, that is we must compute  
    ⟨⟨ 0  Q_ϵ^linS {p_i,α_i}⟩⟩ |_𝒪(g_ YM^3) .
    
    In order to do this, one needs to compute all 𝒪(g_ YM^3) contributions, a calculation which closely parallels the computation of the one-loop soft current <cit.>.  There are additional three-parton contributions, see fig:OneLoopOrtho3part, which combine, using essentially the manipulations from the tree-level calculation, to zero. That is the contributions from fig:OneLoopOrtho3part cancel those of fig:OneLoopSoft. 
    
    The more non-trivial calculation involves the two-parton contributions (see fig:OneLoopSoft2part and fig:OneLoopOrtho) and there are contributions with the same ordering ambiguity as in fig:OneLoopDSoft. Once again if we make the choice corresponding to (8.b.2)  we find that there is a cancellation of the two-parton contributions so that 
    ⟨⟨ 0  Q_ϵ^linS {p_i,α_i}⟩⟩ |_𝒪(g_ YM^3)=0
    
    to leading order in the IR divergences. Alternatively if we choose the ordering (8.b.1) of fig:OneLoopDSoft, we find that all the diagrams that arise from dressing the tree-level amplitude cancel amongst themselves. The diagram corresponding to dressing the one-loop amplitude, (10.a) of fig:OneLoopOrtho3part, cancels the contribution due to the  iterated tree term in the soft-limit of the loop amplitude (<ref>), so that the term arising from the one-loop soft current remains. Thus we find
    
        ⟨⟨ 0  Q_ϵ^linS {p_i,α_i}⟩⟩ |_𝒪(g_ YM^3)=-
        ^2 C_A/16π^2ϵ̂^2∑_ℓ∈in Q_ϵ^h(p_ℓ) ℳ^(0)_n ,
    
    which is the same as the one-loop contribution to the Ward identity for that choice. 
    
    While for simplicity we focused on the case with only incoming hard particles, it is possible to generalize our results to  generic in- and out-states by essentially using crossing symmetry. In order to check diagrammatically that crossing is satisfied for dressed states, one can argue that the non-Abelian gluon clouds at leading logarithmic order weakly commute with the S-matrix. The steps are essentially the same as in <cit.>, with the important difference that here we have non-Abelian generators inside the clouds. However even in the non-Abelian theory the coherent state-operators have the property 
    
    
        [ 𝒰^pE_αβ, 𝒰^p^'E_α^'β^'] = 0 .
    
    We have explicitly checked that for the ordering prescription (8.b.2) there are no corrections to the Ward identity for generic states, while for ordering (8.b.1) the correction term now has a sum over both outgoing and incoming partons with a corresponding sign. 
    
    
    
    § DISCUSSION AND CONCLUSIONS
    
    
    In this work  we have computed matrix elements of the commutator 
    
    [Q_ϵ^lin, S]=Q_ϵ^+,lin S -SQ_ϵ^-,lin
    between asymptotic states corresponding to hard partons dressed with soft-gluon coherent operators. At tree-level the result following from the soft-theorem is essentially the same as in the Abelian case
    
        ⟨⟨out[Q_ϵ^lin, S]in⟩⟩=-[∑_ℓ∈out N^a(ϵ, p_ℓ) t^a_ℓ -∑_ℓ∈in N^a(ϵ, p_ℓ) t^a_ℓ] ⟨⟨out Sin⟩⟩ .
    
    In the Abelian case the right-hand side of this relation corresponds to the contribution from the hard part of the asymptotic charge and so this expression is equivalent to the conservation law for the full charge. The non-Abelian theory is quite different as the correct definition of the non-linear charge is unclear and these  difficulties are likely only exacerbated by including quantum corrections. We choose a pragmatic approach of  defining our quantum non-linear hard charge by means of the soft evolution operators and then computing matrix elements to check whether the resulting Ward identities continue to hold. We have seen in sec:LGTSEv that the hard charge does receive modifications at loop order however they involve soft gluon operators which can be shown to not contribute at one-loop order if we use IR-finite S-matrix elements with a particular ordering prescription. 
    
    In the computation of matrix elements of the commutator we noted that there are several order-of-limits issues. First, we have assumed a “weak" definition of the soft charge insertion, where the soft limit for the gluon operator in the charge is taken after the evaluation of the matrix element.  This makes the connection between the asymptotic charge and the standard soft limit of amplitudes most direct. A similar issue arose in the case of sub-leading soft-theorems for graviton amplitudes and our approach is analogous to  that of <cit.>.
    For the case of the dressed states there remains a second ordering ambiguity, which can be seen in diagrams (8.b.1) and (8.b.2) in fig:OneLoopDSoft, as one must choose whether to take the soft-limit in the charge before—ordering O_1—or after—ordering O_2—the contribution from the coherent state operator.
    
    These two orderings imply:
    
      *  Ordering O_1: The orthogonality relation
    	
        . ⟨⟨ 0 Q_ϵ^lin S in⟩⟩|_𝒪(^3)=-
        	^2 C_A/16π^2ϵ̂^2∑_ℓ∈in Q_ϵ^h(p_ℓ)
        	⟨⟨ 0Sin⟩⟩ .
    
    	is broken by pure one loop effects (i.e. no tree-level iterated terms) and the soft-charge no longer produces states orthogonal to scattering states.  The Ward identity receives the exact same correction and one may attempt to correct the hard charge at one loop in  fashion similar to <cit.> to preserve the conservation law. 
    	
      *  Ordering O_2: Using this definition the orthogonality relation is still valid at one-loop and leading singularity
    	. ⟨⟨ 0 Q_ϵ^lin S in⟩⟩|_𝒪(^3)   =   0+𝒪(1ϵ̂)
    	
    	and moreover the Ward identity is preserved
    	. ⟨⟨out[Q_ϵ^lin,S] in⟩⟩|_𝒪(^3 )   =   0+𝒪(1ϵ̂) .
    	
    	In this case the interpretation is clearer and we can see that loop effects at this order do not affect the asymptotic symmetries of the S-matrix. 
    
    It is quite natural to choose a prescription that preserves the symmetries of the theory where it does not lead to inconsistencies which argues in favor of the second choice. That this second choice is to be preferred on symmetry grounds can be further seen as follows. 
    In this work we have used the dressing following from the evolution operator (<ref>). However one could, following <cit.>, modify the dressing so that
    𝒥^a_q, f_μ= ∫_ω_qdp ρ^a(p)(p_μ/p· q-f_μ/ω_q) ,
    
    where f_μ is a q-dependent vector and whose value could be determined by insisting on physical asymptotic states. Relatedly one can choose a modified vacuum
    0⟩⟩_Λ=P_ωexp(∫^E_λdqΛ^a(q)Π^a(q) )|0 ⟩
    which is an eigenstate of the linearised charge
    
    Q^lin_ϵ 0 ⟩⟩_Λ=Λ̅_ϵ 0 ⟩⟩_Λ
    with 
    Λ̅_ϵ=-4π∫dq δ(ω_q) /√(γ_zz̅)Λ^a (q)
     (ε^- ∂_z ϵ^a  + ε^+ ∂_z̅ϵ^a)
    
    if we also assume that the Fock vacuum is an eigenstate of the linearised charge with eigenvalue 0. Strictly speaking this corresponds to the identification of the Fock vacuum 
    with the dressed vacuum |0⟩≡ 0 ⟩⟩_0 but this is in fact what we have used as we have set Ω^E|0⟩=|0⟩ in our regularised, perturbative calculations. In this sense we may interpret the failure of the orthogonality relation in the case of ordering O_1 as being related to the coherent states no longer being eigenstates of the linearised charge at one-loop and so requiring a one-loop modification. Thus if we make the physically sensible choice that the asymptotic states remain eigenstates of the linearised asymptotic charges we are forced to use ordering O_2 and we thus find that there are no corrections to the Ward identity. 
    
    While our explicit computations have been at one-loop, the leading IR divergence coherent state construction is valid to all-loop orders and so with the correct ordering prescription there should be no leading IR divergent quantum corrections to the Ward identity at any loop order. More non-trivially it should also be possible to repeat the analysis at subleading order in the IR divergences.  The construction of the dressed states is more complicated since single-parton coherent states can no longer be defined and one needs non-factorisable coherent states which take into account two-parton correlations <cit.>.  At sub-leading order the issue of collinear IR-divergences will also become more involved.  In this work we have included collinear effects to the order we consider, that is in the soft-region, which gives rise to the double pole in ϵ̂, and we see that the asymptotic charges continue to be conserved with the appropriate ordering prescription. At subleading order this will need to be reconsidered as there are also collinear singularities in the non-soft region<cit.>. Such collinear divergences play an important role in understanding tree-level asymptotic symmetries where they are related to the algebra of charges, <cit.>. Though at loop-level they are less well studied, the coherent state approach can treat such divergences, see <cit.> and the recent work <cit.> which proposed using the soft-collinear effective theory Hamiltonian to generate the asymptotic time evolution. It can in fact be shown quite generally <cit.> that the dressed S-matrix is completely IR finite and as the proof makes no use of an explicit expression of the asymptotic Hamiltonian it includes any collinear divergences. 
    
    The Faddeev-Kulish dressing of asymptotic particles with soft gauge boson clouds is related to the dressing of hard external particle fields with Wilson lines. This method was first introduced by Mandelstam <cit.> in order to achieve a gauge-invariant formulation of QED. Nowadays, it is, amongst other applications, an important tool to efficiently compute the infrared-divergence structure of scattering amplitudes. See e.g. <cit.> for a discussion of the virtual soft gluon contribution to multi-leg amplitudes in non-Abelian gauge theories obtained in this formalism. It is also known that the Wilson line approach to soft radiation produces the Faddeev-Kulish dressing when taking the Wilson lines along time-like paths <cit.>. Furthermore, it was recently used in <cit.> to study soft photon hair of black holes. It would be interesting to further study the connection between these formulations and asymptotic symmetries in the context of QCD. 
    
    Finally in this work we have focused on single insertions of the asymptotic charge but it would be of interest to consider the generalisation to multiple insertions to study the algebra of asymptotic charges.  He et al  in  <cit.>, see also <cit.>, showed that the tree-level double-soft limit of two positive helicity gluons could be rewritten as the level zero Kac-Moody algebra.  This can in fact even be used to construct a stress-energy tensor for gluons by the Sugawara method and to derive a KZ-like equation for MHV amplitudes <cit.>. By considering more general tree-level double-soft limits <cit.> the algebra of currents can be extended to include mixed-helicity gluons and sub-leading, in the soft expansion, currents <cit.>. The coherent state approach may make it possible to understand quantum corrections to the algebra of such currents. 
    
    
    We would like to thank  Einan Gardi, Mario Greco, Lorenzo Magnea, Dhiritman Nandan and Gim Seng Ng for useful discussions. We thank Dhritiman Nandan for comments on an draft version of the current work. This project has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No. 764850 "SAGEX". T.McL. and A.S. are supported by the SFI grant 15/CDA/3472. D.M. acknowledges support from the Spanish Government grant FPA2016-78022-P and Spanish MINECO Centro de Excelencia Severo Ochoa Programme (SEV-2016-0597).
    
    
    § NOTATIONS
    
      *  We use a mostly positive signature g_μν=diag(-1, +1, +1, +1) so that x^0=ct and p^0=E.
    	
      *  Retarded coordinates are defined as
    	
    	x^μ=(u+r, r(z+z̅)/1+ z z̅, ir(z-z̅)/1+ z z̅, r1-z z̅/1+ z z̅)
    	
    	and polarizations are 
    	ε^μ _-=1/√(2)(z, 1, i , -z) ,   ε^μ _+=1/√(2)(z̅, 1, -i , -z̅) .
    	
    	 The polarisation vectors satisfy
    ∑_σε̅_μ^σ(q)ε_ν^σ(q) =  η_μν+c_μ q_ν+c_ν q_μ ,
    
    where c_μ is a fixed vector which depends on the choice of polarisation vectors. In our conventions for retarded Bondi coordinates we have c_μ=(1+z z̅)2ω(-1,0,0,1).
    	
    	
      *  
    	The free mode expansion of the gluon field is 
    	𝒜^a_μ(x)= ∫dq[ε̅_μ^σ(q⃗) a_σ^a(q⃗)e^iq· x+ε_μ^σ(q⃗)
    	a_σ^a†(q⃗) e^-iq· x] ,
    	
    	where 
    	dq=d^3 q/(2π)^3 (2ω)
    	and we use the commutator 
    	
    	[a_σ^a(q⃗), a_σ'^† b(q⃗')]=δ̃(q⃗-q⃗')δ_σσ'δ^ab
    with δ̃(q⃗-q⃗')=(2π)^3 (2ω)δ^(3)(q⃗-q⃗').
    
      *  The gauge generators of 𝔰𝔲(N) are defined by the relations
    
    [t^a, t^b]=i f^abct^c
    
    and are normalised such that Tr(t^at^b)=12δ^ab. Quarks transform in the fundamental representation t^a=T^a and gluons in the adjoint (t^a)_b c = (T_A^a)_b c=-i f^a_bc.
    	
      *  The covariant derivative is defined as
    	 D_μ=∂_μ-ig_ YM[𝒜_μ,  ].
    
    § ASYMPTOTIC EXPANSIONS
    
    In this appendix we review the construction of the asymptotic charges for non-Abelian gauge theory. We follow <cit.> and the review <cit.>. 
    In order to understand the asymptotic symmetries, one must impose fall-off conditions compatible with the equations of motion and which allow for relevant solutions. We consider Lorenz gauge, which in our coordinates is
    
    
        -∂_u (r^2 𝒜_r) - ∂_r (r^2 𝒜_u - r^2 𝒜_r) + γ^z z̅ (∂_z 𝒜_z̅ + ∂_z̅𝒜_z) = 0
     
    and impose the conditions at null infinity (r, t→∞, u=t-r=const)
          𝒜_z=A_z(u, z, z̅)+𝒪(r^-1) ,   𝒜_z̅= A_z̅(u, z, z̅)+𝒪(r^-1)
          𝒜_r=1/r^2A_r(u, z, z̅)+𝒪(r^-3) ,   𝒜_u=1/rA_u(u, z, z̅)+𝒪(r^-2) .
    
    The corresponding condition on the field strength component is 
    ℱ_ru=r^-2 F_ru+𝒪(r^-3)    with    F_ru=-(A_u+∂_u A_r)
    
    which using the u-component of the field equations, satisfies
    
        ∂_u(A_u+∂_u A_r)- γ^zz̅   ∂_u(∂_z A_z̅+∂_z̅A_z) = -i γ^zz̅([A_z̅, ∂_u A_z]+[A_z, ∂_u A_z̅])+ j_u^(2) ,
    
    where j_u=j_u^(2)/r^2+𝒪(r^-3). Hence one finds
    
    Q_ϵ=∫_ℐ^+_-r^2 d^2 z  γ_zz̅ϵ^a ℱ^a_ru
    giving the expressions (<ref>) and (<ref>). 
    
    Alternatively, one can impose the following conditions on the falloffs of the non-Abelian gauge fields at large-r<cit.>
    
        𝒜_u   =1/r^2 A_u^(2)(u,z,z̅)+log(r)/rÃ_u^(1)(u,z,z̅) +log(r)/r^2Ã_u^(2)(u,z,z̅) +𝒪(log(r)/r^3) ,
        𝒜_r   =1/r^2A_r^(2)(u,z,z̅)+log(r)/r^2Ã_r^(2)(u,z,z̅) +𝒪(log(r)/r^3) ,
        𝒜_z   =A_z^(0)(u,z,z̅)+1/rA_z^(1)(u,z,z̅)+log(r)/rÃ_z^(1)(u,z,z̅)+𝒪(log(r)/r^2) ,
    
    and use the freedom in the residual gauge transformations at order 1/r to set A_u(u,z,z̅) = 0. 
    From (<ref>) we get the leading components for the field strengths
    
    
        ℱ_ur^a   = 1/r^2[∂_u A_r^a,(2) - Ã_u^a,(1)] +log(r)/r^2[∂_u Ã^a,(2)_r + Ã^a,(1)_u ] + 𝒪(log(r)/r^3)
        ℱ_uz^a   =∂_u A_z^a,(0)+𝒪(1/r) , 
        ℱ_rz^a   = 1/r^2[-A_z^a,(2) + Ã_z^a,(1) - ∂_z A^a,(2)_r + g_ YM f^a b c A_r^b,(2) A_z^c,(0)]
           +log(r)/r^2(-Ã^a,(1)_z - ∂_z Ã^a,(2)_r + f^a b cÃ^b,(2)_r A_z^c,(0)) +𝒪(log(r)/r^3) ,  
        ℱ_zz̅^a   =∂_zA_z̅^a,(0)-∂_z̅A_z^a,(0)+ g_ YMf^a b c A_z^b,(0) A_z̅^c,(0) +𝒪(1/r) ,
    
    so that the radiation flux is non-zero and finite on ℐ^+ as required (<cit.>).
    The equations of motion imply the following constraint equations on ℐ^+ at 𝒪(1) in the large-r expansion
    
        -∂_u Ã^a,(1)_u + ∂_u^2 A_r^a,(2)
           = - γ^z z̅[∂_u( ∂_z̅ A^a,(0)_z + ∂_z A^a,(0)_z̅ )+ g_ YM f^a b c(A_z̅^b,(0)∂_u A_z^c,(0) - A_z^c,(0)∂_u A_z̅^b,(0)) ]+  j_u^a,(2)
    
    and 
    
        2 (∂_u A_z^a,(1) - ∂_u Ã_z^a,(1)) + ∂_z (∂_u A_r^a,(2) + Ã_u^a,(1)) - g_ YM f^a b c∂_u (A_r^b,(2) A_z^c,(0))- g_ YM f^a b cÃ_u^b,(1) A_z^c,(0)
           - ∂_z (γ^z z̅ (∂_z A^a,(0)_z̅ - ∂_z̅ A^a,(0)_z + g_ YM f^a b cÃ_z^b,(0) A_z̅^c,(0)))
           +g_ YM f^a b c[-A_r^b,(2)∂_u A_z^c,(0) + γ^z z̅ A_z^b,(0)(∂_z̅ A_z^c,(0) - ∂_z A_z̅^c,(0) + g_ YM f^c d eÃ_z̅^d,(0) A_z^e,(0))) ]=  j_z^a,(2)
    
    and
    
        -∂_u A_r^a,(2) - Ã_u^a,(1) + γ^z z̅ (∂_z A_z̅^a,(0) + ∂_z̅ A_z^a,(0)) = 0
    
    whereas at order 𝒪(log(r))
        2 ∂_u Ã^a,(1)_z - 2 ∂_z Ã^a,(1)_u - g_ YM f^a b c(∂_u Ã_r^b,(2) A_z^c,(0) + Ã^b,(1)_u A_z^c,(0) +Ã^b,(2)_u A_z^c,(0) + Ã^b,(2)_r ∂_u A_z^c,(0))  = 0.
    
    In these expressions we have used that the currents have the decay properties
    
        j_u = 𝒪(1/r^2)      j_r = 𝒪(1/r^3)      j_z , j_z̅ = 𝒪(1/r^2) .
    
    Just as before (<ref>) can be used to rewrite the charge as in (<ref>) and (<ref>). One can also use retarded radial gauge or radiation gauge and at least at leading order one ultimately finds the same expressions
    for the charge. 
    
    
    
    
    *JHEP-2
