Colour Evolution and Infrared Physics

We give a complete account of how soft gluon, massless quark, evolution equations in colour space originate, from a factorization into a hard cross section density operator and a soft function encoding measurements and the projection on definite colours. We detail this formalism up to the two loop level and we demonstrate how the evolution kernels relate to infrared subtractions, and how the resolution of infrared singular regions conspires with the structure of observables the algorithm should be able to predict. The latter allows us to address evolution in different kinematic variables, including energy ordering and angular cutoffs in non-global observables. The soft factor and its evolution resembles a hadronization model including effects such as colour reconnection, and could give insight into the structure of power corrections in observables which require soft gluon evolution.

1 Introduction The reliable prediction of realistic final states produced in high energy collisions is at the heart of both, event generator simulation as well as, in a less versatile but typically more precise fashion, analytic resummation.Recent work has been highlighting the basic fact that both approaches need to rely on analyzing amplitudes with many external particles in order to fully account how the final state builds up in detail.Cross sections are then obtained by squaring such an amplitude and in the presence of a (certain class of) measurements leading contributions can then be accounted for.A simplification is possible by using QCD coherence for the resummation of global observables, or dipole algorithms suited for non-global observables in the large-N limit, though in general more complicated evolution equations remain to be addressed.
Parton branching at the amplitude level [1][2][3] has become an important theoretical framework to construct and analyze parton shower algorithms.These approaches generalize the soft gluon evolution algorithm discussed in [4], which has recently been used to resum non-global logarithms beyond the leading-N limit [5,6].Understanding the structures in the soft limit unveils most of the complexity we need to account for in improved algorithms since collinear physics is typically colour diagonal and only soft physics can guide us to what structures new parton showers need to build on [7,8].In particular, we have recently been pointing out that there exist colour-diagonal, though formally 1/N suppressed contributions in the two-loop and one-loop/one-emission contributions [9].In the same work, we have also laid the ground to analyse virtual corrections differential in phase-space type integrals without resorting to unitarity; this allows us to extract imaginary parts as well and gain insight into the distributional structure of the integrands, as has recently also been pointed out in [10] in the context of an effective field theory analysis [11].The physics of non-global observables has also been addressed starting from a duality to small-x evolution in [12], and the resummation of non-global logarithms at the NLL level has also been achieved by [13,14].The purpose of this work is to approach this resummation program, at the same level, from the colour evolution point of view.We will derive the algorithm we have previously been studying at the LL level [4] that is implemented for full colour evolution in the CVolver library [5,6], however we go beyond this in three aspects: We will formulate the structure of the evolution at the next order to incorporate the two loop ingredients [9] and double soft currents [15,16]; we will point out how hadronization can enter the calculation of cross sections, making firmer links to our previous observation that part of the physics modeled by colour reconnection and cluster fission might actually stem from a perturbative interface to colour evolution [17,18]; and we will analyse how the new evolution factor impacts the accuracy and attempt to relate our analysis to factorization theorems to make contact with first principle analyses and other approaches.Once our soft algorithm is extended to the hard collinear case, a task which mainly consists of a more complicated book keeping, we will be able to address more general algorithms which then would account for the splitting functions we have identified in [19].Definitions of singular phase space regions become crucial to derive the evolution equations and generalize the usual parton shower infrared cutoff.We show that variations of the partonic shower cutoff necessarily force the soft function (and thus a hadronization model) to compensate for changes in the hard evolution.This might connect to previous and ongoing work on the interpretation of the top mass parameter [20], as well as similar EFT approaches to hadronization corrections following the concepts set out in [21].While we have now gained a significant amount of understanding for the QCD case, electroweak physics only now is recognized as in need for a treatment beyond the customary high energy evolution addressing the quasi-collinear limit.In recent work we have pointed out out that an amplitude factorization similar to QCD can serve as a starting point for evolution [22] and part of the present work's motivation is to address electroweak evolution.This is particularly tied to the question of what role the measurement and the projection on the physical final state will play (also see [23] for a new conceptual approach to this question that highlights the similarity with the need for hadronization in QCD).
Figure 1: Graphical representation of the present approach.A cross section is calculated by considering a partonic amplitude and its conjugate as an operator in colour space of intermediate partonic configurations.We sum over all configurations weighted by a projection on final state colour structures and a measurement, integrated over the finally observed objects.The details of the measurement are not of concern to the present work but we will investigate to what extent this factorization dictates evolution equations for both M n and U n .

Aim of the present work and guide to the paper
Our starting point will be a cross section definition which factors hard scattering matrix elements from possibly non-perturbative factors in the amplitude which describe the formation of hadrons and the implementation of the observable of interest.We shall perform this analysis on the level of "density operators" in colour space, as we have previously studied for partonic cross sections in earlier work (see e.g.[3,4,9] for the general formalism and to to set the notation).Graphically we represent our approach in Fig. 1.We assume that we will be able to calculate the scattering matrix element and its conjugate, conveniently assembled into the cross section density operator M n = |M n M n | in perturbation theory and as a vector |M n in the space of colour structures for n additionally emitted partons on top of a hard scattering process.The cross section density operator, in a given basis, would thus be represented as a quadratic matrix [4].The combined effect of projecting the partonic final state on certain colour configurations, and on performing a measurement described by an observable function on the final state, will then be encoded in an operator U n , which contains the truly quantum mechanical external states and operators which overlap with these states in the sense that we can create and annihilate (hard) partons from these states.The form of the resulting cross section will be (schematically) where dφ n indicates an integration over the partonic momenta (effectively, the partonic phase space, though this is already a very specific assumption).U n might be structured such that it allows or prevents cancellation of infrared singularities.In the latter case we need to consider U n as a genuinely non-perturbative object such as a fragmentation function or parton distribution function is.The aim of the present paper is not to work out the complete definition of U n and the conditions under which we can actually write a cross section in the form above (we will give an outlook on this topic), but rather what the implications of the shear presence of a non-trivial U n are: we analyse to what extent we can use the anticipated factorization to consistently construct evolution equations -and this mainly means finite, iterative algorithms -to build up the hard scattering operator and the measurement operator (which we can think of as anything in between a jet measurement, fragmentation function, or a completely exclusive hadronization model), subject to a given class of observables.Any finite algorithm, if it shall be based on events of fixed multiplicity, necessarily involves an infrared resolution.The presence of this resolution, e.g. through an infrared cutoff scale µ S , shall be our starting point.We will devise re-definitions (or renormalization transformations) from M n and U n , onto finite density operators A(µ S ) = (A 0 (µ S ), A 1 (µ S ), ...) and effective measurement operators S(µ S ) = (S 0 (µ S ), S 1 (µ S ), ...) such that the cross section is invariant in the sense that where α 0 = α S Z g relates the bare and renormalized couplings, and µ S is a resolution scale (in fact, as we demonstrate below, a collection of resolutions scales), of which the final cross sections needs to be independent in the same way as it is of the renormalization scale (which we have supressed for simplicity in these introductory notes).Jet cross sections would be recovered if S n are the unit operators in colour space times a scalar measurement function, and would then be represented by a scalar product matrix of overlaps of colour states, see [4] for more details.The independence of the "bare" objects M n and U n of the scale µ S will result in evolution equations for the A n (µ S ) and S n (µ S ), and also implies the fact that the cross section is independent of the chosen resolution scale.Had we chosen to work with vectors of density operators A rather their components A n for each partonic multiplicity, then the action of Z and X would in fact be matrices of colour charge operators which are inverse to each other.While the former notation might look more transparent on conceptual grounds, we here aim at a practical analysis, which is better based on the density operators for individual partonic multiplicities and the cross feed between them, mediated by the emission of additional partons.We should thus warn the reader that the implementations of the re-definitions Z and X in the main text will appear in a more complicated fashion than by representing them as simple linear operators acting on all the density operators.However we exactly exploit the fact that X is inverse to Z, and we spell out and prove this relation in very detail in the main text: We will actually determine Z from X to be its inverse, possibly order-by-order in the strong coupling α S .X itself needs to be determined to provide infrared subtractions to the hard density operator, after ultraviolet renormalization, and its factorization properties.This step very much resembles what happens in a fixed-order calculation in which infrared divergences are subtracted either by explicit subtraction terms or by re-defining objects like a PDF or fragmentation function.In particular this procedure needs to demonstrate how iterated, or uncorrelated, unresolved emissions are removed from simultaneous unresolved limits of two or more emissions, such as e.g.discussed in the context of the work presented in [13,14].We explicitly demonstrate that our formalism is able to achieve this and it thus gives rise to a finite A, which is given by the correlated matrix elements at the relevant orders without contaminating these by iterated contributions.This rest of this work is organized as follows: In Sec. 2 we will start from the postulate of a factorization theorem and show how the renormalization of a soft factor gives rise to systematic infrared subtraction terms, thus recovering the structure of a fixed-order calculation.From the very fact that we want to eventually calculate a finite cross section we deduce in Sec. 3, how the hard cross section density matrix is factored into a finite contribution which accounts for large logarithms, and renormalization factors which absorb the universal infrared divergences exhibited by the subtraction terms.The former object will satisfy the evolution equation recently proposed in [4] and generalized to the hard collinear case in [3].This evolution equation, as well as the evolution equation for the soft function will be discussed in detail in Sec. 4, where we provide all necessary ingredients up to the second order in α S .To show how we will combine real and virtual corrections subject to a resolution criterion driven by the class of observables we aim at, we give explicit one-loop results in Sec. 5, while a full application of the two loop results will be saved for a follow-up publication.In Sec. 6 we comment on the structure of the second order evolution, and highlight how the colour structure and the resolution criteria are constrained by a given observable.In Sec. 7 we finally comment on how the soft function can also be viewed as a hadronization model including colour reconnection, effectively generalizing the role fragmentation or parton distribution functions would play in canceling out infrared divergences.

Infrared Subtractions and Measurement Operators
We will only focus on soft gluon evolution, and generalize the starting point of the algorithm presented in [4] to the following Ansatz by writing a cross section as where we consider the emission of n additional soft partons on top of a set of hard momenta collectively referred to as Q.We assume that the additional partons are ordered (in the simplest case in energy), such that a full cross section might need to sum over different permutations of orderings, but this does not make a difference for our arguments below.Effectively, the ordering reflects a hierarchy upon which we can factorize the amplitudes.Specifically it does so in a way that if a certain number of last emissions with respect to this ordering cannot be singular through some phase space restriction, that any emission before will not lead to divergences of the same kind.Virtuality, transverse momenta or similar variables [1] might be used to achieve this, but also more complicated combinations of energy and angles can serve the same purpose.θ functions implementing the ordering are thus always implicit when we refer to additional emissions and the ordering corresponds to the labeling of the momenta.α 0 is the bare QCD coupling, and µ the 't Hooft mass as we are working in dimensional regularization with d = 4 − 2 dimensions. [ denote the integration over final state momenta, but we will use the same notation to denote integration over loop momenta, and δ is the mass-shell restriction of the physical phase space. 1 In writing this formula we have assumed a factorization of the cross section in a sense that the measurement definition, and in particular infrared and possibly nonperturbative phenomena are contained in U n .This operator can encode jet cross sections when we choose it to be 1 n × u(Q; p 1 , ..., p n ), the identity operator in colour space times a scalar observable function.If we were to include hard-collinear contributions we could also address cross sections for identified hadrons in case of which the observable function would constrain partonic momenta to be aligned along a certain direction.In the latter case we would have assumed that there are further, unobserved, hadrons for which the accompanying partonic activity would provide the colour charges needed to neutralize colour of observed hadronic final states.It is this mechanism which will break down if we ask for a completely exclusive final state, or even for an identified hadron with a specific quark content.In both of the latter cases projections of colour into certain singlet systems is needed and U, at any fixed partonic multiplicity, cannot be the identity operator in colour space.This shall be our starting point.In the conclusions, Sec 8, we will sketch how this picture can arise from a first principle analysis, work which we will defer in detail to an upcoming paper.

Infrared subtractions for the hard density operator
The hard density operator M n is calculable in perturbation theory as where each term M (l) n refers to the product of an amplitude and its conjugate (see [3,4,9,19] for more details and aspects of this formalism) for n emissions, and l loops.The sum of the amplitude's and conjugate amplitude's number of loops equals l and it is understood that they are distributed symmetrically.For example a one-loop term would read in terms of tree-level and one-loop amplitudes n , respectively.Using the renormalized coupling at a renormalization scale µ R , we will remove the ultraviolet divergences from M n2 .However, subject to the structure of the observable, infrared divergences might partially or fully cancel, or remain all together.We will therefore start to individually remove these divergences with a subtraction formalism which we express by re-defining the observable operator as3 The dependence on µ S appears in all of the quantities S n , X n and F (s) † n+s S n+s F n+s on the right hand side and has only been supressed for readability.In the redefinition, will be used as infrared counter-terms for virtual corrections in the n-parton amplitude, X n , and for s unresolved partons emitted from an n parton state, F n .Here, the • denotes that the operator is understood to act on the left, and on the right of a given operator in colour space: Subject to fixing a particular basis in colour space (or possibly including other quantum numbers as well [22]), the X n operators mix colour basis structures (and would hence be represented as quadratic matrices acting on the colour space for n partons), while emission operators such as F s n do change to a larger parton ensemble, hence a larger basis of colour states and as such would be represented as rectangular matrices.In the case of only soft gluons, or other means to entirely work at the level of factorizing amplitudes [22], there is in principle no reason to link both operators, though we choose to do so here.We will explicitly keep X n and n to not be related by unitarity, since in presence of a definite projection on final states virtual and real contributions might not relate one-to-one.In Sec.6.2 we discuss one example where the unitarity assumption might lead to a problematic algorithm.In particular we expect that this will be the case when considering electroweak evolution when cuts through exchanges of vector bosons involve finite width effects, see [22] for more details.In general, we can only demand S to define the physical observable for which we will be able to calculate a finite cross section.However, if it is possible to demand that also U reflects a physical measurement, then we can use Eq.2.6 to assess the accuracy of the subtraction we perform.Eq. 2.6 reflects by what amount the original definition of the measurement is altered with respect to the re-defined one after subtraction, something which we will be discussing in more detail in Sec. 6.
We will in particular address this case in relation to hadronization corrections for colour evolution and their connection to hadronization models in Sec. 7, though a generalization to fragmentation functions and parton distribution functions will become obvious by then, as well.
To see how the re-definition of the measurement Eq. 2.6 facilitates the task of providing infrared subtractions we expand the cross section Eq. 2.1 in the renormalized coupling Eq. 2.5.This can of course be done at a fixed, minimum, multiplicity, in case of which we reproduce the structure of a fixed-order calculation.Anticipating that we do want to construct an all-emission, all-orders cross section it is more instructive to consider fixed order expansions at each multiplicity, i.e. expanding up to including O(α n+nc S ) at each, fixed n ≥ 0. Once these contributions are rendered finite by the subtraction, the all-order expression which we will derive in the next section will be finite at all orders which iterate O(α nc S ) building blocks.We therefore consider to re-write the cross section as where Mn = l≥0 α l S M(l) n .The individual M(l) n can then be found by comparing fixed orders at fixed multiplicities, where we only consider l ≤ n c and s+l ≤ n c for the expansion of the operators X (l) n and F (s,l) n , respectively.In particular we find for n c = 2 M( 1) tree M (1) tree M (2) where n,R = M (1)  n + nZ (1)  g M (0) n (2.12) and define the renormalized density operators n,R with UV finite virtual corrections up to two loops, and it is understood that M (l) n is zero whenever n < 0 or l < 0. The operators S nc L (with L = tree, 1-loop, ...) act on specific contributions n,R to the hard density operator with n legs and l loops to provide infrared subtractions required up to an expansion to O(α nc S ) in the sense of Eq. 2.8.The result of S is the one pertaining to a density operator of n partons.Had we chosen to collect density operators per partonic multiplicity in an infinite vector, then S would be a matrix which mixes different multiplicities though we refrain from such a presentation since the change between different partonic multiplicities is an important i j respectively.For further details see the main text.The reader is also referred to [9,15,16] for further reference in the context of soft gluons, and to [19] for emission operators beyond this limit.
concept in the problems we analyze in the present work.It is also important to note that the Mn are now finite through the subtractions at each multiplicity, but only up to the additional orders α nc S we have considered in the present analysis.Information about higher orders will only be available once we have facilitated the re-definition of M n to A n thus completing the factorization and renormalization program for the cross section.To this end we have distinguished Mn from A n , though they will coincide if the Mn have effectively been build up by the evolution equation to be put in place for A n later.We in particular find that S tree provides subtractions for unresolved limits stemming from tree-level emission diagrams, for which we can choose F (1,0) n to be given by the single soft emission current D(1,0) n (see Ap.A for more details).Notice that this is an over-subtraction in the double unresolved limit, where H (l) n−1 admits factorization of another emission, as explicitly analyzed in Ap. A. We therefore need to take (2.15) in terms of the true unresolved behaviour D(2,0) n as e.g.given by the soft gluon currents [15].Diagramatically, these contributions have been depicted in Fig. 2, among some of their real/virtual or purely virtual counter parts, which will be dicsussed later in the text.We also require a resolution criterion Θ which projects onto singular regions: Θ n,k is supposed to be one for any singular configuration of k out of n partons (and zero outside), and this in particular implies that S (2) tree [H (l) n ] vanishes whenever n partons are fully resolved.Let us illustrate the requirements of these resolution criteria for example when cutting on energies and angles of the emitted partons.We would then choose or, for two emissions, where in both cases Θn,1 and Θn,2 indicate resolved angular separations: they will vanish on those phase space regions which allow for collinear divergences, e.g. one could choose however all of these choices need to be carefully subjected to analysing the classes of observables we intent to reliably predict.This is further detailed in Sec. 6.Also notice that we strictly have no need to use step functions.Had we chosen resolution in virtualities, then one choice of resolution function could be for one emission (where s ij = (p i + p j ) 2 ), and possible generalizations for more than one emission.Again, the resolution function would be unity whenever any propagator would become singular.Cutting in virtuality for more than one emission is somewhat more complicated, though the requirement that Θ n,2 becomes one whenever any of the singular limits is reached, indicates that also here a product over cuts on all virtualities of the emissions would be sufficient.In what we have presented so far, the Θ n,s functions have been applied as global factors to the emission operators, or as global insertions into the integrands of virtual matrix elements (see below).There is, however, no need to limit this to this case and one can, as well, understand this prescription to act differently for different contributions (even individual graphs like those depicted in Fig. 2).As an example, in the case of virtuality cuts, we could group the contributions in terms of divergencies involving a certain set of invariants, say those resembling the behaviour of the fifth graph in 2 would cut on Θ had we chosen in this case to attach the gluon n − 1 to the line k on the right hand side of the diagram.Generally, D(s,l) n originates from factorizing soft gluon emissions from the renormalized hard density matrix operators and is thus free of UV divergences (see App.A for a precise definition of what factorization properties we assume here).For loop integrals this in particular implies that their genuine UV divergences have been removed and the corresponding integrations are cut off accordingly (more aspects of this will be highlighted in future work).We further discuss the removal of over-subtraction in App.A, and stress the fact the subtraction terms are constructed such that the subtracted quantities are entirely projected on those phase space regions which are free of singularities.We might split up the action of the emission operators into different terms and use different Θ functions e.g.depending on the dipole in between a gluon is exchanged, though we will not make this explicit here for the sake of readability.Furthermore provides one-loop subtractions, and provides two-loop subtractions.Similarly to what happened in the emission case, we need to adjust X n to be given by the singular (or unresolved) part of the one-loop exchange V (1) n [Ξ n,1 ], where Ξ n,k refers to a restriction on k unresolved loop momentum modes in presence of n other resolved momenta, with the shorthand V(l) n ≡ V(l) n [1].ξ = 0, 1 labels two alternatives of the subtraction, and both yield a finite algorithm.Their difference amounts to whether the matrix element from which virtual exchanges have been factors is rendered finite for emissions, as well (ξ = 1) or not (ξ = 0).We can therefore think of ξ = 0 amounting to an inclusive algorithm, for which we would allow the observable function to be infrared safe and allow for unresolved radiation, and ξ = 1 would refer to a fully exclusive algorithm in which we will not rely on any cancellation in between different partonic multiplicities.The latter effectively means that we can achieve a merging of different fixed-order calculations, though we will not explore the consequences further in the present work.More details are given in App. A. Following [9], we use a sum over cuts α and write the virtual exchanges on external lines as where Ξ n,l is the analogue of Θ n,k for the virtual corrections (though we, again, might split this up into more terms, or combine cuts with each other).Since any more complicated situation follows by generalization we again keep things simple and refer to one function Ξ n,l in the following.I (l) n,α (p 1 , ..., p n ; k 1 , ..., k l ) denotes the integrand, which might include δ-functions which pin down the loop momenta to particular kinematic regions.They thus are distributions in general, including the i0 prescriptions.For the two-loop case we use in terms of the singular parts of the two-loop exchanges V(2) n (note that the X n • X (1) † n subtractions automatically adjust for the two loop case).These virtual corrections again stem from factorization of the renormalized matrix elements, and it is thus understood that suitable UV counter terms have been added to their definition, or that the genuine UV divergent integrals are cut off accordingly (artificial UV divergences might arise from the soft gluon approximation and are understood to be cut off by the resolution functions).Some examples are again given in Fig. 2. E.g. the one-loop self energy insertion on a one-loop soft gluon exchange would require a UV counter term which resembles Z (1) g upon integration (such counter terms can readily be constructed by using the methods developed in e.g.[24,25]).No other genuine UV divergences arise in the two-loop case.
Constructing counter terms for the one-loop/one-emission contributions which are free of over-subtractions we find in terms of the full unresolved behaviour D(1,1) n (and similarly for F ), which does not exactly coincide with the one-loop soft gluon current of [16] which in turn has some iterated contributions already removed.However it can easily be constructed from this result.On the very same note V(2) n and D(2,0) n are understood to contain the full sum of contributing diagrams, see App.A for details, and Fig. 2 for sample diagrams.

Resummation
In the regions of unresolved emissions and exchanges we expect the hard density operator to factorize as where It is important to stress that this re-definition is to be thought of in the sense that we renormalize the divergent density operator M n and trade it off for a finite, though possibly logarithmically enhanced, object A n from which we can hope to build up an algorithm for the calculation of physical, infrared finite, cross sections.This operation needs to constitute the inverse of the transformation providing infrared subtractions such that the cross section stays invariant under the combined transformation.We thus will be able to infer the operators Z n and E n from their counterparts X n and F n .Also note that the left hand side of Eq. 3.1 is ultraviolet finite, such that we solely accumulate infrared divergences in the renormalization factors on the right hand side.We can then perform fixed-order expansions of M using Eq.2.3, and The cross section reads where ∆ n encodes the effects beyond the exact cancellation of the combined redefinitions Eq. 2.6 and Eq.3.1, and Eq.2.5.We can enforce ∆ n = 0, if holds to all orders in α S .This is the previously mentioned invariance of the cross section under the combined transformation of the hard scattering and measurement operators: transforming M n into A n is precisely inverse to transforming U n into S n .To this extend note that, since the trace is cyclic, the cross section involves a sum over all mutliplicities, and phase space factorizes, the action of re-defining U n in terms of S n can be recast into an action on the M n as where equality in the above holds upon performing the trace, the summation over all mutliplicities and the integration over phase space.Had we chosen to arrange all M n and all U n each into a vector with elements indexed by partonic multiplicity, the above would then be a linear operation, though several of the key features of the evolution will then again be hidden.Inserting the re-definition of M n Z n g directly shows that the operations are inverse to each other so long Eq.3.5 holds.In the below we will determine Z and E in this way order by order, however the formula above might serve as a more general starting point and proves to be especially relevant in the electroweak case, in which F might involve decay processes only (corresponding to observed final states), while E might involve both, emission of unstable particles which subsequently count as hard lines, as well as their decays, and we need to match these identities at the level of different (final) states rather than at fixed order.
In QCD, at fixed order, we simply adjust the renormalization constants as then ∆ n will contribute at O(α n+2 S ) at fixed S n , and with Z (2) ∆ n will contribute at O(α n+3 S ) at fixed S n .Together with the removal of the oversubtractions we then obtain Note that we have left expressions for Θ n,1 (1 − Θ n,1 ) or Θ 2 n,1 to fully keep track of the distributional nature, noting that e.g.θ(x)(1 − θ(x)) = 0 everywhere except x = 0, while still (d/dx)θ(x)(1 − θ(x)) = 0, and also of course 1 − θ(x) + θ(x) = 1.Other choices of resolution functions, however, might thus also be possible to use.
Let us stress the importance of the over-subtraction removal, e.g. in the case of the real emission.Iterating the redefinition of M for two emissions from a hard configuration with no emissions, we would have + D (2,0) 0 as well as where we have put 1 is manifestly finite in the singly-unresolved limit where . We also have In the doubly unresolved limit both

Hence the above gives rise to
Thus, both the presence and the choice of resolution function on the blue contribution are needed such that A (0) 2 is finite (even in a singly unresolved limit when Θ 2,2 → 1 and Θ 2,1 → 1) and solely -without double counting -given by the regularized doubleemission factors and finite in all doubly unresolved limits.No finite algorithm could be achieved had we left out the blue contribution.It is in particular also the presence of the blue contribution which guarantees that in the ordered phase space regions iterated contributions are removed, while we need to consider the full matrix elements elsewhere (to be discussed in more details in the next section).Similar remarks apply to the other contributions, and in Ap.A we give a detailed analysis of this construction.

Evolution Equations
We will now study how the dependence on the different infrared resolution scales as well as the renormalization scale µ R gives rise to evolution equations for the hard and soft factors.Notice that the running of α S is given by the d-dimensional β-function as and it is understood that β S,k = 0 if S = R indicates that we differentiate to a different scale µ S in place of the renormalization scale µ R .From Eq. 2.5 we then have at second order β R,0 = 2 Z (1)

Hard density operator
If we differentiate Eq. 3.1 to any of the resolution scales µ S contained in Θ n,s and Ξ n,s , or the renormalization scale, we then obtain (∂ S ≡ ∂/∂ log µ S as above) where .
Explicit expressions up to second order are obtained in App.B.1.Note that the appearance of β S,−1 will mostly be canceled by derivatives acting on the µ R factors, thus no explicit account of their effect will be needed.For the cases where an explicit account of it is needed we employ (in a distributional sense for functions with support x ∈ (0, 1]) where the + prescription is to be taken with a subtraction at x = 0. Notice however, that most of such cases will actually be dropping out due to the d-dimensional β function.

Soft function
The transformed soft function is subject to an evolution where we now find the anomalous dimensions and evolution kernels to be given by and Explicit expressions up to second order are again given in App.B.2.A final remark is in order with respect to the choices of ξ.For ξ = 1, we will obtain the anomalous dimensions we would have used for ξ = 0, accompanied by additional factors which guarantee that the emissions in the state which the anomalous dimension acts on, are all resolved.The second contribution then is a purely finite virtual contribution acting by pinning down the emissions to the resolution cutoff by means of a derivative of the additional resolution function.The latter contribution thus gives a specific contribution to configurations at the infrared cutoff, i.e. those which do not radiate any further.We will investigate these specific aspects in future work, but in the following only consider the case ξ = 0, which also leads to a finite evolution.As a particular example, consider the one loop anomalous dimension (the real emission at this order is not affected by the different choices of ξ).Here we find, considering a soft evolution variable with β S = 0 only, Γ

Examples for Evolution at Leading Order
In this section we consider the resummation of non-global logarithms in an ordering with respect to energy.Specifically we require virtual counter terms which can be inferred from the cuts obtained with the methods in [9].The one-loop virtual correction is given by and admits the cuts given in [9], out of which the double cut vanishes in dimensional regularization and a real-valued and imaginary single cut remains, if i and j are both incoming or both outgoing, and the absorptive cut is canceled in the case of an incoming/outgoing pair.These give rise to the subtractions (notice that cutting introduces a minus sign in the use of the Feynman tree theorem) X n,abs = i<j In this case we have assumed that ξ = 0, and indicated that we have evaluated these integrals in the dipole rest frame, which introduces additional complications when combining with the real emission.In the dipole frame, and for the radiative cut, we have x = (p i − p j ) • p n+1 / 2p i • p j and E = (p i + p j ) • p n+1 / 2p i • p j with p 2 n+1 = 0, and writing the integral in the above form requires Ξ(ij) n,1,rad to be symmetric around x = 0.For the absorptive part, Ex = p j •p n+1 / 2p i • p j and p i •p n+1 = 0.A boost back to lab frame needs to separately map the dipole's legs into each other, Λ ij ni,j = n i,j , where p i,j = E i,j n i,j , such that the cancellation with the real emission can be ensured.In the case of a cutoff on energies we tag the exchange as resolved, whenever its energy is within the real emission phase space (here indicated by the limit on Q), or it has a minimum energy or a minimum opening angle to one dipole leg, Then we find for the anomalous dimensions from the radiative cut by differentiating with respect to µ S ∂/∂µ S and ∂/∂λ Notice that the anomalous dimension for the angular cutoff would exhibit a pole in if it was not for the upper bound on the energy.We strongly suggest that this upper bound should not only be inferred by considerations about the real cancellation, but by amplitudelevel considerations about the ordering variable as advocated in [4,26,27].The role of the boost and frame dependence of the collinear cutoff will be discussed in a separate work, notice only that allows to translate between the lab frame and the frame in which we consider regularizing the virtuals.The UV running at first order is given by confirming the role of the hard and soft density operators.The absorptive cut is regularized through the i0 terms which do not cancel, and is integrable for any resolution which is symmetric around, and not vanishing at x = 0 (in fact, an angular cutoff is not needed), provided that the energy is less than another 'Glauber' scale µ G : 3) . (5.9) The real emission subtraction would be given by the square of the single soft gluon current, and accordingly be supplemented by a resolution criterion which we take to be and the derivatives analogously give the contributions to the evolution from the real emissions.The evolution in the energy scale then provides us with the soft gluon algorithm studied in [4]; the evolution direction with the cutoff has not yet been considered completely but will likely cancel for observables in which collinear divergences cancel.We will leave this to future work. 4 Accuracy Considerations

Evolution at the Second Order and Comparison to Shower Approaches
In this section we compare the structure of the evolution at the second order to approaches of parton shower algorithms which aim to describe evolution at the same accuracy.The results here can also be used to assemble a full evolution using our results [9] and the double emission current [15].Putting β S,i = 0 for the evolution of the hard density operator we obtain Γ (2) Notice that both the double real, as well as the double virtual terms have a form which resemble a subtraction calculation in the sense that the "earlier" exchange or emission n − 1 constrained to be resolved by 1 − Ξ n,1 or 1 − Θ n−1,1 and the later one is pinned down to the evolution variable.Also Eq. 6.2 shows a similar structure, and in fact the first contribution removes the one-emission current after a virtual insertion from the oneloop/one-emission contribution and so is in one-to-one correspondence with directly using the result of [16], while the second term corresponds to a different subtraction for the correction of an emission which is resolved.It would thus be interesting to compare the structure of such subtractions with NLO-type subtraction scheme to O(α 2 S ) corrections of a splitting kernel as advocated in [30,31].Let us in particular emphasize how the resolution conspires with the subtraction in the cases above, e.g. for the real emission.In this case we can take Θ n,1 = 1 − Θn,1 θ(E n − µ S ) with some angular resolution Θn,1 which regulates collinear divergences.In order to guarantee that the double emission phase space is projected onto the singular region, we need to take Θ n,2 = 1− Θn,2 θ(E n−1 −µ S )θ(E n −µ S ).Inserting this into Eq.6.3 we find that i.e. each of the two terms is completely regulated by the angular and energy cutoffs, and subtractions of the iterated piece are applied in the ordered region only.
of the observable.In this case, the resolution criterion Θ must be accordingly synchronized with the observable in question.
where we have generically identified the collinear cutoff prescription with Θ λ (n n+1 ) and the emission's direction with n n+1 .We assume that the observable probes a veto on energies less than ρ in the angular out region, complementary to the in region.If we identify the energy resolution µ S with the observable resolution ρ then in fact the difference vanishes except if the emission becomes collinearly unresolved in the out region.In this case we expect that the change in between different collinear cutoffs will need to be compensated for by using the cutoff anomalous dimension Eq. 5.6 to obtain a stable result.

Aspects of a Hadronization Model
We will now focus on the evolution equation for S n .In fact it is interesting to make explicit that this object is expected to convert the hard partons plus n emissions into m hadrons eventually, and thus we write S n → S nm to highlight this fact and assume that observables will eventually be obtained at hadron level by integrating over the m final state momenta.Notice that the soft evolution is quite different in its form as compared to the evolution equation of A n .In particular, the cross-feed with different partonic multiplicities is now coming in from larger multiplicities, with the emission operators appearing in a form that they project down the n + 1 onto the n-parton colour space; also the anomalous dimensions are now appearing conjugate to the case of the operator A n , as otherwise they would not be able to be used as subtractions for the virtual contributions.Any hadronization picture consistent with amplitude level evolution must satisfy the evolution equation above; in what follows we will now mainly focus on initial conditions and how this perturbative evolution might tie in with existing models.To be complete, though, let us write the solution (in four dimensions) as where is the usual soft gluon evolution operator also entering the evolution of the operator A n , and we have fixed an initial condition We have also set µ S = µ R , thus including a running coupling which could be continued to a meaningful infrared model.Algorithmically, including such a model is not much different than the algorithms we have already presented in [6].In fact, one would now generate the evolution of the operator A n , at the amplitude level, down to an infrared scale µ S and in the leading order case considered here the evolution would essentially continue unaltered to even lower scales where it is terminated by an insertion of S NP nm .This simple pattern will be altered at higher orders; it should otherwise also not be considered an unnecessary complication since the above factorization truly shows that one can retain perturbative control of the soft scale µ S (meaning that our final cross section is independent of it to the order we have pinned the evolution down to) and then systematically hand over to the evolution from a high-energy 'end' of a hadronization model including effects such as colour reconnection, which is mediated by the evolution operators just as it is in the context of the evolution discussed in [18], and possibly related to further model improvements discussed in the literature [17,34].We should also note that depending on the Ansatz chosen for the non-perturbative initial condition such an algorithm might not only be very efficient, but the evolution needs to be inserted in order to arrive at a configuration we can associate with a hadronic final state.
As a start we will try and use such a picture to make contact with the cluster hadronization model [35][36][37][38].Let us therefore assume that the initial conditions for S nm will describe the transition from n partons into n − 1 primary clusters5 (assuming a q q + (n − 2)g final state), and so we can take where Cl n (w) projects onto the valid cluster configurations given by all colour flows which do not involve a 'singlet' gluon contribution, but including, and explicitly allow, states where the q q pair is colour connected and accompanied by a gluonic system forming a singlet on its own.These are the configurations which explicitly cannot originate from a leading-N , probabilistic parton shower.We then take Cl n (w) to be diagonal, so we can simply weight each cluster configuration by a distribution of the cluster decay products in phase space, say [τ |Cl n (w)|σ] = δ τ σ w σ δ cluster(σ) , where we take to describe the individual cluster decay weights.The emission operators will then map such a configuration onto one which is contained by merging clusters, and we have where we have abbreviated the Eikonal factor by ω ij , i.e.
where wn (w) links to w by a simple leading-N dipole transition, which in this context we in fact associate with the cluster fission process.The evolution operators will then mediate a colour reconnection dynamics as discussed in [17], however in this case the process will happen differently in the amplitude and the conjugate amplitude.In the colour flow picture, the 'singlet' gluon (or trace condition part) would become inert to this evolution just as it does in the hard evolution.This signals that we need to extend the non-perturbative initial condition in fact to include a process which describes the transition of the 'singlet' colour flow gluon into a single cluster.

Summary and Outlook
In this letter we have discussed how evolution equations in colour space, specifically those for soft gluons, originate from a general diagrammatic recursion, infrared subtraction, and the application of a renormalization program to re-arrange the individual contributions into a finite cross section.We have particularly paid attention to the role of the observable and the accuracy which the algorithm can deliver, and to the consequences of observables which are forced to resolve colour, either because we are considering hadronization corrections or factorize the partonic contributions in a more general way than needed for a jet cross section.This provided us with a new notion of the structure of power corrections which now appear encoded in an operator relation in colour space.They will only simplify to more standard approaches in case of jet cross sections and by enforcing unitarity in the sense that infrared subtractions should relate between real and virtual contributions in a one-to-one correspondence.We stress that, due to the nature of unstable particles in the case of electroweak evolution [22], this can not be achieved any longer since real and virtual corrections are not lined up anymore, and explicit projections on isospin states are performed by the measurement.The role of the measurement needs to be clarified as an integral part of the evolution.We will address this aspect in more detail in future work.The role of the soft, or measurement operator, is to encode how virtual and real contributions are re-arranged to cancel in the physical cross section.For QCD, this corresponds precisely to what we expect a hadronization model to do since in the presence of such dynamics infrared cancellations are obscured.On a more general ground, we can adopt the following analysis of the cross section, which we here present in a more sketchy way with details to be addressed later: is the general cross section we would strive to calculate, in terms of m experimentally observed particles from which we can calculate an observable u(φ m ), and the trace refers to a sum over degrees of freedom within the n constituent particles from which we could build up the m observed particles.Up until now we have collectively addressed as the measurement function.The above clearly must be a consequence of a factorization theorem which would start from general cross section and then exploits the properties of (classes of) hadronic observables.We therefore assume that our formalism will provide a generalization of cross sections with many identified hadrons or the track function formalism [39] to cases beyond purely collinear physics and in particular if correlations correlators among the hadrons are involved.Ultimately, we would start from a fully exclusive final state with S matrix element squared of the form Our current formalism has then assumed that the result will be dominated by those contributions which have all partonic lines put on-shell, and where each of the P i is given by the sum of constituents j p ij and the two factors are then convoluted by momentum integrations whereas the remaining integration then is over the partonic phase space Subject to the observable and leading regions the convolution measure can then be simplified further to yield the soft gluon algorithm we have been discussing now.It is worth noting that R leaves (d − 1) ( m i=1 ni − m) degrees of freedom to integrate over.Essentially the soft gluon algorithm can be constructed by transforming, for n i = ni , the pij into scaled versions close to the p ij and then to expand around scaling parameters near one (see [19,22] for more discussion on such mappings) which then do not appear in the conjugate amplitude anymore and can easily be integrated over.While have been making no assumptions about the precise form of U and have studied its mere presence, a careful analysis of the above then clearly reveals an operator definition.Also note that in this context, the soft factor can comprise operators with more particles overlapping with the physical, external states -this exactly corresponds to the ability of S n to absorb emissions which have been generated during the evolution of the hard process as the emission contribution in the evolution of S n corresponds to a transition from a larger to a smaller final state.
Beyond the soft gluon case we would need to allow the algorithm to populate momentum configurations and flavours possibly differently in between the amplitude and its conjugate so long as they serve to provide constituents consistent with the measured final states (in the sense of interpolating fields with non-vanishing overlaps).This might already apply in the case of hard collinear physics, where convolutions not only exist in the longitudinal momentum fractions but also regarding transverse momenta.The expansion around the partonic mass shells is then still possible, and a treatment of the algorithms along these lines will become crucial in the electroweak case.It will heavily rely on the fact that we maintain overall energy-momentum conservation in the amplitude and the conjugate, as advocated in [22].The main purpose of this analysis and outlook is, however, to show that we can touch ground with first-principle definitions, an avenue which we will explore further.Beyond this observation we note that correlations among hadrons with certain momentum balances might actually be able to probe whatever (colour) interference the hard process allows to, since this will be imprinted in the evolution much in the same way as it has been observed in the case of double parton distributions for those processes where colour interference will mix these objects [40].In fact, the general formula we envisage should be easily able to accommodate processes with incoming hadrons.

m
with k ≤ l loops and/or m ≤ n legs, but the colour structure of S (nc) l H (l) n

Figure 2 :
Figure 2: Contributions to the different emission and virtual correction operators.From top left to bottom right we show graphs contributing to D

R
m [dp] n [dp] {p} n |{P } m ) χβ ({p} n|{P } m ) G α (p n ) Ḡβ (p n)m-(truncated, on-shell)(8.3)where the sum over the field indices α, β does constitute Tr n , and the relevant Green's function is a tensor in all of these indices, subject to choosing a basis |α} of colour and spin (which we here have deliberately denoted with a curly bracket notation to distinguish them from the true quantum mechanical states): |m = |i |f is the product of initial and final states involved in the definition of the S matrix element, and χ α (p n ) are collectively denoting the external wave functions, i.e. in general contractions of interpolating fields with the physical external states, or (conjugates of) Bethe-Salpeter amplitudes for the states of interest.They are dictated by the interpolating fields ψ α we have chosen for the elementary or possibly composite external states, and one can also consider the optical theorem in order to arrive at a similar structure.R m is the product of the corresponding wave function renormalizations for the final state m, and the integrations over the offshell constituent momenta are constrained such that their sum equal the observed final state momenta P i .The truncation of the Green's functions is to be understood such that iterations of N m -PI-irreducible legs which combine into a composite state of N m constituents have been truncated.Momentum conservation of the observed particles then results in an integration over [dP ] m δ In the case of a prototypical non-global observable we can writeu(p 1 , ..., p n , p n+1 ) = (θ(ρ − p 0 n+1 ) + θ(p 0 n+1 − ρ)θ(n n+1 ∈ in))u(p 1 , ..., p n ) , (u(p 1 , ..., p n , p n+1 ) − u(p 1 , ..., p n