Towards all-order factorization of QED amplitudes at next-to-leading power

We generalise the factorization of abelian gauge theory amplitudes to next-to-leading power (NLP) in a soft scale expansion, following a recent generalisation for Yukawa theory. From an all-order power counting analysis of leading and next-to-leading regions, we infer the factorized structure for both a parametrically small and zero fermion mass. This requires the introduction of new universal jet functions, for non-radiative and single-radiative QED amplitudes, which we compute at one-loop order. We show that our factorization formula reproduces the relevant regions in one- and two-loop scattering amplitudes, appropriately addressing endpoint divergences. It provides a description of virtual collinear modes and accounts for non-trivial hard-collinear interplay present beyond the one-loop level, making this a first step towards a complete all-order factorization framework for gauge-theory amplitudes at NLP.


Introduction
Deepening our understanding of gauge theory scattering amplitudes in the limit where radiation is soft has important phenomenological benefits as well as significant intrinsic value. For n-particle scattering processes in QED with the emission of an additional soft photon with momentum k (as in figure 1), the scattering amplitude M n+1 can be expressed as a power expansion in the energy E = k 0 of the photon, where the leading power (LP) term has scaling M LP n+1 ∼ 1/E, and the next-to-leading power (NLP) contribution is of order M NLP n+1 ∼ E 0 . Crucially, the coefficients in this expansion can be expressed in terms of simpler objects, which relate the radiative amplitude to the non-radiative or elastic amplitude M n . Such relations go under the name of factorization or soft theorems. Their physical interpretation rests upon the long wavelength of soft radiation not being able to resolve the hard scattering. However, at each subsequent power in (1.1) more is revealed. In this paper we investigate aspects of factorization at NLP, and focus in particular on the objects that can appear at higher orders in perturbation theory. The simpler objects that enter in the coefficients of (1.1) describe soft and collinear dynamics and are universal, i.e. they do not depend on the particular scattering process. For instance, it is well known that the LP term in eq. (1.1) takes the universal form [1,2] where p µ i and q i denote the momentum and electric charge (in units of the elementary charge e) of the i-th hard particle, and ε µ (k) is the polarisation vector of the soft photon. The soft function S (0) n describes a set of eikonal interactions between the external particles and the emitted soft photon; in other words, soft radiation at LP is sensitive only to the direction and charge of the emitting particle. (The expression in eq. (1.2) is at lowest order in the coupling e, as indicated by the superscript (0), and receives loop corrections.) For multiple soft photons the function S n can be calculated as the vacuum expectation value of a set of Wilson lines, one for each hard emitting particle, expanded to the appropriate order in the coupling.
The factorization in eq. (1.2) is not only of theoretical interest, but also relevant for phenomenology. The 1/E singularity in the soft limit enhances soft radiation in scattering processes. Measurements that are sensitive to soft radiation involve a small scale, and the corresponding cross section contains large logarithms of the ratio of this small scale and the scale of the hard scattering. Such large logarithms potentially spoil the convergence of the expansion in the coupling e, a problem that can be addressed by resummation. The development of theorems such as the one in eq. (1.2) led to the proofs of factorization [3][4][5], but also constitutes the first step towards resummation, as it allows one to decompose a multi-scale amplitude (or cross section) into the product of simpler single-scale functions.
The resummation of large logarithms in QCD originating from the LP term in the equivalent expansion of eq. (1.1) has been an active research topic for many years. Resummation of soft gluon radiation at LP has been systematically applied to most processes of interest at lepton and hadron colliders. In the seminal papers [6,7] soft gluon resummation in Drell-Yan and DIS was achieved by means of diagrammatic techniques, and in [8,9] it was shown that soft radiation can be described in terms of Wilson lines, whose exponentiation properties are at the basis of resummation. Later, by means of similar diagrammatic techniques, resummation was extended to more processes, including those with coloured particles in the final state, see e.g. [10][11][12][13][14][15]. Soft gluon resummation by means of renormalisation-group techniques was first studied in [16] and in a different method in [17,18], and this has more recently also been accomplished using effective field theory techniques, see e.g. [19][20][21][22][23][24][25][26].
By contrast, the factorization and resummation of the NLP contribution in eq. (1.1) is still under much investigation. This NLP term exhibits a more involved structure: emitted (next-to-)soft radiation becomes sensitive to the spin of the hard particles, and starts to reveal details of the internal structure of the hard interaction. The first insight into the structure of M NLP n+1 was already achieved a long time ago in papers by Low, Burnett and Kroll [27,28], who realised that, in the case of massive emitting particles, the structure of the NLP term is dictated by gauge invariance, by means of Ward identities. This early formulation, now known as the "LBK" theorem, was proven [29] to hold only in the region k 0 m 2 /Q, with Q the centre of mass energy. For m 2 /Q < k 0 < m, the LBK theorem must be extended to account for NLP contributions arising from soft photons emitted from loops in which the exchanged virtual particles have momenta collinear to the external particles (i.e. having a small virtuality, while retaining momentum components which are large compared to the soft radiation), which can be taken into account by a radiative jet [29].
The factorization proposed in [29] was later confirmed by some of us in [30], and extended to non-abelian theories in [31], unveiling many features of soft radiation at NLP. For instance, NLP radiation at next-to-leading order (NLO) in perturbation theory has a universal structure, where the radiative matrix element squared can be expressed as a reweighing of a kinematics-shifted non-radiative matrix element squared. This reweighing factor is universal, in the sense that it only depends on the (colour) charge of the particles participating in the scattering [32,33]. These developments enabled the soft gluon resummation at NLP at leading logarithmic (LL) accuracy for Drell-Yan [34] (which has also been achieved using effective field theory methods [35], discussed more below), proving earlier conjectures [36][37][38]. Other methods, based on evolution equations, have been used in [39][40][41][42]. A phenomenological study of NLP resummation at LL accuracy for prompt photon production was carried out in [43].
While the radiative jet function significantly extends the applicability of the factorization beyond the original LBK theorem, there are additional types of radiative jets beyond one loop. These describe for example a soft emission from multiple highly energetic particles in the same collinear direction emanating from the hard scattering, which has been studied for Yukawa theory in [44]. The extension of this analysis to QED, classifying all types of radiative jets, is necessary to describe soft radiation in QED at NLP to all orders in perturbation theory, and is addressed in this paper.
Factorization theorems for amplitudes or cross sections involving soft emissions have also been studied using an effective field theory approach, specifically within Soft-Collinear Effective Theory (SCET) [45][46][47][48][49]. SCET describes the soft and collinear limits of QCD as separate degrees of freedom, each with their own Lagrangian. The elastic amplitude M n in the factorization in (1.2) is encoded in terms of effective n-jet operators and their corresponding short-distance coefficients, which capture the contribution from hard loops. At LP, soft emissions from hard particles can be described by Wilson lines, as in the diagrammatic picture. Beyond LP, these soft emissions follow from time-ordered non-local operators made out of soft and collinear fields, with insertions of the power-suppressed soft and collinear Lagrangian. Within SCET it is possible to define matrix elements which are equivalent to the radiative jets of the diagrammatic approach [50][51][52]. Several investigations have been conducted within SCET, including but not limited to soft gluon corrections, such as studies of the anomalous dimension of power-suppressed operators [53][54][55], the basis of power-suppressed hard-scattering operators for several processes [56][57][58], the application to subtractions [59][60][61][62][63][64][65] and resummation of NLP LLs in a variety of processes [35,[66][67][68][69][70][71][72].
SCET provides a systematic approach, as each operator and Lagrangian term have by construction a definite power counting. When all operators are included that are consistent with symmetries up to the desired power, this completeness ensures that the resulting factorization is valid to all orders in the coupling constant. Consequently, factorization can, at least formally, be extended beyond NLP, by simply adding more power-suppressed operators. On the other hand, the diagrammatic approach is often more direct compared to the full effective field theory treatment, and may also offer a way (as we will see in this paper) to address so-called endpoint singularities in convolution integrals. These convolutions between ingredients in the factorization are only well defined in dimensional regularisation, thus posing a challenge for SCET, where one first renormalises each ingredient in the factorization theorem to derive the renormalization group equations needed for resummation, causing these convolution integrals to become divergent. (However see [69] for recent progress in addressing this issue.) One may hope that the diagrammatic approach will provide an easier path to resummation, as resummation exploits exponentiation properties of soft radiation and the replica trick [73,74], which can be carried out within dimensional regularisation.
In order to set the stage for the development of a factorization theorem for the radiative amplitude M NLP n+1 in eq. (1.1), let us recall in more detail where the original LBK theorem fails. Within the latter treatment, one separates the radiative amplitude in two contributions: one in which the radiation is emitted from the external legs, plus another in which the radiation is emitted from a particle within the hard scattering kernel. Schematically this is written as where the two terms are represented in figure 2. The amplitude M int n+1 can be obtained by means of the Ward identity from M ext n+1 . Concerning the amplitude involving an emission from the external legs, consider as an example the emission from an outgoing fermion i: in this case M ext n+1 takes the form where M n represents the elastic amplitude (stripped off the spinorū(p i )). Within the LBK theorem, one expands the amplitude in the soft momentum k. As discussed, this expansion gives correct results only in the regime p i · k/Q m 2 /Q. But for parametrically small masses m 2 /Q < p i · k/Q < m, naively expanding the elastic amplitude in the soft momentum k misses a contribution in which the soft photon is emitted from internal particles collinear to the external leg, which is thus included neither in M ext n+1 , nor in M int n+1 . From a practical perspective, it is known that a correct expansion of M n+1 can be obtained by means of the method of regions [75,76]: one splits a given loop integral into the sum of several integrals, in which the loop momentum is assumed to take different scalings, related to the scaling of the momenta of the external particles. Within each region one is allowed to expand the integrand in the small scales in that region, and then the sum over all regions is expected to provide the correct expansion of the original integral. In particular, using this method it is possible to check that the missing contribution in eq. (1.4) indeed arises from the region where the loop momentum has scaling collinear to the external particles [77,78].
From the point of view of factorization, eq. (1.4) tells us that in order to understand the factorization properties of the radiative amplitude M n+1 we also need to obtain the correct factorization structure of the elastic amplitude M n , in presence of a small offshellness p i · k ∼ m 2 Q 2 . At leading power the factorization structure is known, see for instance [79,80], and it takes the following schematic form: In this equation the jet and soft functions J j and S describe long-distance collinear and soft virtual radiation in M n . These functions are universal, i.e. they depend only on the colour and spin quantum numbers of the external states, and determine also the structure of collinear and soft singularities of the elastic amplitude. Note that one must divide each jet by its eikonal counterpart J i , to avoid the double counting of soft and collinear divergences. Given this premise, our first task is thus to determine the analogue of eq. (1.5) at NLP. In the absence of soft radiation, the elastic amplitude would only depend on hard scales. Following [44], we will consider the external fermions to have a parametrically small mass m, providing us with a variable for the power expansion. We derive the power counting, which we then generalise to the case of massless particles, and obtain an allorder NLP factorization formula. In either case new, universal jet functions are required at the NLP level, consisting of multiple collinear particles along the same direction, probing the hard scattering process. We restrict ourselves to the first non-trivial jet function, which we calculate for a parametrically small fermion mass up to NLP. Subsequently, we perform checks at one-and two-loop level that validate the obtained jet function and the corresponding part of the factorization formula. We carry out a similar analysis for single-radiative amplitudes in the massless fermion scenario. We stress that our study is exploratory in nature, a full characterisation of the radiative jet functions is left to future work. Having in hand the factorization for both massless and massive fermions would allow the application of our results to a larger class of scattering processes of interest at the LHC, including the production of heavy coloured particles, such as top quarks, or scalar quarks and gluinos in supersymmetry.
The outline of the paper is as follows: in section 2 we carry out the power counting analysis that underpins the (non-radiative) all-order NLP factorization formula presented in that section for two scenarios: one with a parametrically small fermion mass m and one for massless fermions. We focus on the massive case in section 3, computing the first non-trivial jet function and performing checks at the one-and two-loop level. In section 4, we consider single-radiative amplitudes in the massless fermion scenario, again performing one-and two-loop checks. We conclude in section 5, while certain technical aspects are relegated to appendices.

From power counting to factorization
The factorization of an n-particle scattering amplitude with emission of a soft photon M n+1 crucially depends on the factorization properties of the corresponding elastic amplitude M n . Therefore we start our analysis by extending eq. (1.5) to NLP. Specifically, we set out to obtain a classification of the jet-like structures, consisting of virtual radiation collinear to any of the n external hard particles, contributing at subleading power. Phrased differently, we wish to derive which jet functions contribute up to NLP in a parametrically small scale, corresponding to a fermion mass or a soft external momentum.
In the following we will distinguish two fermion mass scenarios. One of which is the truly massless theory (m = 0), the standard approximation in high energy calculations, where it is well understood that (virtual) collinear effects beyond the LBK theorem play an important role. In the other scenario, we consider fermion masses to be non-zero but parametrically small, and in fact comparable to the scale associated to soft emissions: we assume m ∼ λQ, such that p i · k ∼ m 2 ∼ λ 2 Q 2 Q 2 , where k is a soft momentum. This more intricate small-mass approximation could be of phenomenological importance if soft gluons are emitted from particles with an intermediate-size mass, as mass effects may be comparable in size to the aforementioned collinear effects. Resummation of resulting NLP threshold logarithms (beyond LL accuracy) would, in that case, require a proper understanding of massive radiative jet functions. In addition, this second scenario may prove useful for the resummation of logarithmic mass terms, log(m/Q), even in non-radiative processes, where the small fermion mass m and the hard scale Q are the only scales in the problem.
We derive our results by power counting the pinch surfaces, that underlie the collinear (and soft) contributions we wish to describe in terms of jet (and soft) functions, for a general QED scattering amplitude. This was done recently for Yukawa theory in [44] for the same two mass scenarios. The pinch surfaces are the solutions of the Landau equations [81] and are represented by reduced diagrams in the Coleman-Norton picture [82]. In these diagrams, all off-shell lines are shrunk to a point, while the on-shell lines are kept and may be organised according to the nature of the singularity they embody, be it soft or collinear. This results in the general reduced diagram of figure 3, in which one distinguishes a soft "blob" containing all lines carrying solely soft momentum, n jets J i comprised of lines with momenta collinear to the respective external parton and lastly, a hard blob H collecting all contracted, off-shell lines. This picture seems unaltered by the presence of parametrically small fermion masses, because the limit of small m yields the same singular pinch surfaces as the massless theory. In support of this claim, we analysed the QED massive form factor using the method of regions [75,76,83], finding that soft and collinear modes are sufficient to correctly reproduce the singularity structure in this limit. This analysis is presented in appendix A.
To carry out the power counting we use light-cone coordinates. For each external momentum p i , we introduce two light-like vectors n i andn i , defined by These vectors are normalised such that n i ·n i = 1, and by definition n 2 i =n 2 i = 0. In any of these coordinates, a generic vector v decomposes as where v + = v·n i , v − = v·n i . A scalar product of two vectors then reads . . . Figure 3: The reduced diagram for a general, vector boson induced QED process with n well-separated hard particles in the final state. Ellipses denote the presence of an arbitrary number of photon/(anti-)fermion lines.
Of course, this needs further specification in which of the N collinear directions the decomposition is carried out. Adopting the notation k µ = (k + , k ⊥ , k − ), we associate the following scaling to lines that are soft or collinear to the i-th external leg Soft: k µ ∼ Q λ 2 , λ 2 , λ 2 , Collinear: in terms of the light-cone coordinates corresponding to p i . The scaling of the normal coordinates parametrises the contribution of soft and collinear lines around the singular surface, which is reached for λ → 0. Away from this limit, power counting in λ thus amounts to the ordering of finite contributions of different size and proves to be a valuable technique. We focus on virtual corrections to a hard scattering configuration, for which all invariants s ij = (p i + p j ) 2 ∼ Q 2 involving external momenta are large compared to the energy of the radiated soft photon in M n+1 . Requiring the soft momentum to be of order λ 2 rather than λ guarantees that the photon is soft with respect to all particles in the elastic amplitude. Whenever we refer to a NLP quantity in this paper, we mean that it is suppressed by up to two powers in λ with respect to the leading power contribution. This nomenclature originates from strictly massless (m = 0) QED where power corrections arise only through scales associated to soft emissions p i · k ∼ λ 2 Q 2 . In case of parametrically small masses (m ∼ λQ) power suppressed terms at O(λ) do occur, but we apply the same definition nonetheless.
Using the momentum scaling in eq. (2.4), we start by deriving the superficial degree of divergence of a particular reduced diagram G contained in figure 3, which is simply the λ-scaling of this diagram, G ∼ λ γ G . Specifically, we wish to determine how γ G depends on the structure of G. We will see that, in practice, γ G can be expressed as function of the number of fermion and photon connections between the hard, soft and collinear subgraphs and, in presence of fermion mass, on the internal structure of the soft subgraph. Such a formula tells us, at any perturbative order, which pinch surfaces contribute up to NLP and guides us in setting up a consistent and complete NLP factorization framework for QED. 1

Power counting rules for individual components
In order to derive an expression for γ G it is convenient to set up a catalogue of the degree of divergence of the individual components first. For massive (m ∼ λQ) and massless (m = 0) QED, these rules vary slightly and we derive them explicitly here.
Given eq. (2.4), the propagator for a collinear, massive fermion scales as A massless collinear fermion obeys the same rule as the mass term is subleading in the numerator (O(λ) versus O(1)) and of equal size in the denominator (both O λ 2 ). For soft fermion lines a difference does arise; for non-zero mass while for a massless fermion one finds instead Since we aim at determining the order at which each configuration start contributing, we will only keep track of the most singular contribution to γ G , and discard the subleading terms in eq. (2.5) and (2.6). The singular structure of eq. (2.6) is uncommon because the denominator is not strictly on shell, since p 2 ∼ λ 4 while m 2 ∼ λ 2 . In fact, this singularity is entirely determined by the fermion mass. Intuitively, because of their mass, soft fermions are integrated out, an aspect that would be worth investigating from an effective theory perspective. This momentum configuration, which contributes to the singular structure of scattering amplitudes despite being off shell, bears similarity to Glauber gluons, scaling as (λ 2 , λ, λ 2 ). Our power counting shows that these momentum configurations could affect scattering amplitudes only beyond NLP.
In gauge theories the rules for vector boson vertices depend on the choice of gauge. For power counting purposes, the axial gauge is particularly convenient since non-physical degrees of freedom do not propagate. The latter is a direct consequence of the form of the photon propagator, which reads 8) 1 We stress that the power counting analysis performed here is analogous to the one for Yukawa theory, which underlies the results of ref. [44]. Nonetheless, we deem it instructive to show this derivation in detail. We summarise additional results for Yukawa theory in appendix D.
with the choice of the reference vector r µ fixing the gauge. Eq. (2.8) satisfies while contracting with the propagating momentum results in which no longer has a pole in k 2 . Together, eqs. (2.9) and (2.10) show that scalar and longitudinal polarisations do not propagate in the chosen gauge. This choice of gauge makes it particularly convenient to derive the suppression effects associated to vertices involving gauge bosons [84]. These are effective rules, in the sense that they are not evident from the vertex factor as obtained from the QED Lagrangian, but follow from an interplay with the adjacent lines. To make this concrete, consider the expression for the emission of a photon from a collinear fermion line with momentum p µ , which is proportional to First, we point out that the first two terms are always power suppressed: the first one is per definition of order λ 2 , while the second term is of order λ even if the photon emission is collinear, as the dominant component vanishes due to (γ − ) 2 = 0 (for a soft photon the second term is manifestly of order λ 2 ). If the photon is soft, p µ / p in the third term of eq. (2.11) dominates, being of order λ 0 . In that case, no suppression is caused by the vertex. However, if the photon is collinear to the fermion lines extending from the vertex, we can write p µ = p + k + k µ + O(λ). From eq. (2.10) we then conclude that there is no dominant contribution to on-shell scattering amplitudes from the third term in eq. (2.11). Hence, in axial gauge, a suppression of λ is associated to each emission of a collinear photon from a collinear fermion line. With a different choice of gauge, the presence of longitudinal polarisations would erase this suppression effect of vertices, and individual diagrams would exhibit a harder scaling. In physical observables such polarisations cancel due to Ward identities, and the extra λ suppression would become evident when summing over a gaugeinvariant set of diagrams. We summarise the rules for QED vertices in table 1. The scaling of a photon line is determined by the common factor 1 k 2 in eq. (2.8), and is therefore ∼ 1 λ 2 and ∼ 1 λ 4 for respectively collinear and soft particles. A further suppression of the degree of divergence results from integration over loop momenta, where the measure d + d − d 2 ⊥ provides a suppression of, respectively, λ 4 and λ 8 for collinear and soft loops. These results are, together with the rules for propagators, presented in table 2.

Constructing the overall degree of divergence
We started the derivation of a formula for γ G by obtaining the power counting rules for the basic constituents of any diagram. Here we use these results to obtain power counting formulae for the soft (γ S ) and collinear (γ J i ) sub-diagrams independently. Subsequently, we consider the effect of connections between all sub-diagrams (γ S↔H ,γ J i ↔H and Table 1: Power counting rules for QED vertices, depending on the soft or collinear nature of the field. These rules apply to massive and massless fermions alike. Table 2: Power counting rules for loop integrals and propagators for photons and fermions. If no rule is specified for m ∼ λQ, the scaling is identical to m = 0. as well as the connections to the external particles (γ ext J i ). The degree of divergence of a reduced diagram G with n-jets will thus be given by We begin with γ J i and consider a blob of collinear lines, without any external attachments. According to the rules of table 2 the associated degree of divergence is where I =Ĩ f +Ĩ γ denotes the total number of fermion and photon lines internal to the isolated blob, L the number of loops and V the number of vertices. We use Euler's identity and note that diagrams without external legs (i.e. vacuum bubbles) have three internal lines per pair of vertices: I = 3 2 V . As a result, the degree of divergence of a collinear sub-diagram is independent of its internal structure: For the soft sub-diagram one needs to distinguish between the different mass cases. We start with  Applying Euler's identity in eq. (2.14) and exploiting the fixed ratio of the number of fermion and photon lines to the number of vertices in a QED vacuum bubbleĨ f = 2Ĩ γ = V , we obtain Next, we must account for the contribution to the overall degree of divergence arising from the connecting lines between hard, soft and jet sub-diagrams of the general reduced diagram in figure 3. Besides the explicit powers of λ associated to lines themselves, they affect the power counting of the disconnected sub-diagrams by splitting internal propagators and adding vertices to both sub-diagrams. In figure 4 we show these effects on a generic sub-diagram A resulting from either a photon or fermion connection to a sub-diagram B, depending on the internal line that is probed. For fermion connections, a fermion antifermion pair is inserted to conserve charge in both sub-diagrams. 2 The effect per fermion is simply half that of the combined fermion anti-fermion insertion. An additional suppression effect arises from the loops that are formed in this process.
Consider the connection between a jet and the hard sub-diagram first. A connecting (collinear) photon line adds also a collinear fermion line and an all-collinear QED vertex to the collinear blob, as shown in figure 4a. According to the rules listed in table 2 and 1, such a connection enhances the degree of divergence by −2 − 2 + 1 = −3. Each connecting fermion line gives the same effect, as found by using the aforementioned procedure: a fermion anti-fermion pair adds in total four collinear lines and two all-collinear vertices, such that the enhancement of the degree of divergence due to a single fermion line is f connecting lines give rise to N (i) − 1 collinear loops. Summing up, we find The reduced diagrams considered here are amputated, meaning that there is no propagator, and thus no power counting, associated to the external leg itself. Therefore, connecting the jet to an external fermion leg gives a further enhancement of the degree of divergence of where only the vertices and additional collinear propagators due to the (pairwise) fermion insertion are counted. Similarly, connecting the jet to an external photon gives one additional collinear fermion line an an all-collinear vertex, such that also in this case In contrast to γ J i ↔H and γ ext J i , the degree of divergence associated to the connection between the soft and hard sub-diagram is not suppressed by vertices since all lines are soft. Therefore, we only need to count the m γ + m f soft connections themselves, as well as the additional lines created in the soft blob by these insertions (one soft fermion per photon insertion; one soft photon and an additional soft fermion for a pairwise fermion insertion). 3 Including the loop suppression, we find Finally, we consider the n f connections between the soft sub-diagram and the jets, which affect both the sub-diagrams involved. 4 Also, these connections will form an additional loop by closing a path through H, J i and S, giving a total of n The number of soft photon and fermion lines connecting to the hard sub-diagram, denoted by mγ and m f , should not be confused with the fermion mass m. 4 The main difference in power counting compared to Yukawa theory arises from this interaction. In Yukawa theory, each scalar emission from a collinear fermion line is suppressed by a factor of λ, such that power counting rules for all-collinear and all-soft vertices are identical in QED and Yukawa theory. However, vertices for soft-collinear interactions are suppressed by λ in Yukawa theory, but are not suppressed in QED.
Combining ingredients according to eq. (2.12) gives We emphasise that the number of internal fermion linesĨ f in the soft sub-diagram denotes the number of lines in the isolated blob, before connections to the hard and jet functions have been accounted for. It is more intuitive to express this in terms of the total number of internal fermion lines in the amputated soft function, I f , for which we disregard the actual fermion connections to other blobs, but retain the effect that the connections have on the soft blob itself. Either a single photon attachment or a pairwise (anti-)fermion insertion adds a fermion line to the soft sub-diagram, as indicated in figure 4, giving the relation Inserting eq. (2.24) in eq. (2.23) gives which are the final expressions for the overall degree of divergence for a reduced diagram G with n-jets. The massless result in eq. (2.25a) is the analogue for N -jet production in QED of the power-counting formulae first derived in [84,85] for cut vacuum polarisation diagrams and wide-angle scattering amplitudes in a broader class of theories. The massive result in eq. (2.25b) is the equivalent of the equation obtained in [44] for Yukawa theory, which we also re-derived. We present this and other results for Yukawa theory in appendix D.

NLP factorization of QED amplitudes
Equipped with eq. (2.25), we can determine which reduced diagrams G contribute up to NLP in λ. For the class of diagrams considered in the previous section, which have an arbitrary number of purely virtual corrections, we see that γ G ≥ 0, independent of the number of hard particles in the final state. The γ G = 0 diagrams contain at most logarithmic singularities, while the γ G > 0 are finite and give a vanishing contribution in the λ → 0 limit. For small but non-zero values of λ, the γ G = 0 diagrams form LP contributions, with the γ G > 0 diagrams acting as power corrections. Eventually, we wish to develop a factorization formalism that allows one to resum NLP threshold logarithms associated to soft final-state radiation, which requires us to study the factorization of radiative amplitudes. Dressing the non-radiated graphs with a single, soft emission will enhance the degree of divergence by −2, by the splitting of a soft/collinear fermion line. 5 So for these radiative amplitudes, a LP contribution will be O(λ −2 ) instead, with the NLP corrections of O(λ 0 ). 6 Therefore, we will list all purely virtual reduced diagrams G characterised by γ G ≤ 2. Since we study the abelian theory, we restrict our analysis to (anti-)fermions in the final state, although the power counting formulae of eq. (2.25) describe processes involving hard final-state photons as well. As a minimal example we study the amplitude for γ → ff , but stress that the jet functions that appear there cover the general case of n (anti-)fermions.
The leading power configuration at γ G = 0 is obtained for N      at NLP only in case m = 0, while for m ∼ λQ eq. (2.25b) yields γ G = 5. 7 In either mass scenario, the soft blob may be connected to the jets by an arbitrary number of photons.
Since the reduced diagrams of figure 5 and figure 6 encode all relevant soft and collinear configurations up to NLP, we may immediately cast them into entries in the factorization formula. Starting at leading power, figure 5a yields the well-known factorization formula where the tensor product ⊗ denotes a contraction of spinor indices. The hatted vectors contain the dominant momentum component onlŷ where the light-cone vector n µ i is defined in eq. (2.1). The jet function has the operator definition while the soft function S is given by a product of Wilson lines, For simplicity, we assume that the potential overlap between the soft and collinear regions has already been accounted for in a redefinition of the jet functions.
Following the reasoning of [44] for Yukawa theory, we assume that a similar factorization picture holds at next-to-leading power, with each class of reduced diagrams described by a different jet function. As far as the hard-collinear sector is concerned, this means that the leading power formula in eq. (2.26) is supplemented with four types of contributions, To improve readability, we suppress the arguments of the factorization ingredients and introduce the indices i, j, labeling the collinear sectors. We will clarify this notation further momentarily. The first line describes the effect of figure 5b and starts contributing at order λ. This implies that at order λ 2 we may expect a dependence of the hard function on the perpendicular momentum component of the collinear photon emerging from it, which can be re-expressed in terms of the H i (f ∂γ) function, as will be shown shortly. The second line describes the classes of diagrams (a) and (b) in figure 6, while diagram (c) corresponds to the third line. These contributions, as well as the f ∂γ-term, are strictly O(λ 2 ), which implies that the soft function appearing in those terms is given by the leading-power definition of eq. (2.30). While for massless fermions the same reasoning applies to the f γ-term, in the massive case the soft function could in principle receive O(λ) corrections. Since we focus on hard-collinear factorization, we do not explore this possibility in detail. For the same reason, we will not supplement our factorization formula with terms corresponding to reduced diagrams (d) − (f) with additional connections to the soft function. We leave the identification and investigation of the corresponding terms for future work. Eq. (2.31) is formally identical to the counterpart for massive Yukawa theory [44], as the collinear sectors of the two theories exhibit the same scaling modulo the replacement of scalars with photons.
We now clarify the shorthand index notation. In the simplest non-trivial example of the f γ-jet and hard functions, we define The arguments of the jet function denote the momentum flowing through the fermion and photon leg, respectively, while in the hard function the index i also specifies which of the n hard momenta has been shifted in presence of the additional collinear emission. In analogy with eq. (2.27),ˆ µ i = + i n µ i denotes the large component of the momentum flowing in the photon leg. In principle, in the spirit of the LP factorization, one would like to replace p i withp i in the argument of the hard function, thus neglecting the small components in the external momenta too. This can be done in the massless theory, where the jet functions start contributing at O(λ 2 ). However, the massive theory allows for odd powers in the λ expansion, so that an overall O(λ 2 ) term can also originate from an order λ correction from both the hard and a jet function. This effect forces us to retain some subleading components in the argument of eq. (2.32), as will be made clear in the explicit calculation in section 3.2.1. In contrast to eq. (2.26), the ⊗-product in eq. (2.31) involves, besides spinor index contractions, convolutions over the leading momentum components and additional Lorentz contractions over spacetime indices carried by the photon leg. Explicitly, for the first line of eq. (2.31) The other terms in eq. (2.31) involve a straightforward generalisation of the notation in eq. (2.32). In presence of more than two legs (as for f γγ), the corresponding hard function acquires an additional argument, and the p i are shifted accordingly. As is clear from eq. (2.33), the hard functions depend only on the large momentum componentˆ (and not on the full ). However, since we want the NLP formula to be accurate at O(λ 2 ), we cannot set µ =ˆ µ at the level of amplitudes, but we need to keep also its transverse component µ ⊥ . This can be rephrased as a Taylor expansion in the transverse momentum around zero, thus identifying the two terms with respectively the f γ-and f ∂γ-contributions in eq. (2.31), where by definition the ρ ⊥ in the second term is absorbed in J i (f ∂γ) . In eq. (2.34), we generically denoted with H the part of the amplitude that is not explicitly described by the soft and collinear functions. In the traditional factorization approach, this would be obtained via a subtraction algorithm, while in the effective field theory it corresponds to a Wilson coefficient obtained from matching to full QED. Both approaches would require matrix element definitions of the jet functions in eq. (2.31), as well as of the NLP soft function.
Gauge invariance of each separate ingredient would then be manifest. This systematic analysis requires further investigation of the interplay between jet functions and (generalised) Wilson lines, which we leave for future work. In absence of an operator definition, we will extract in sec. 3.1 the f γ-and f ∂γ-jet functions from a generic matrix element, assuming the validity of the picture above. This necessitates in turn a diagrammatic definition of the hard function, of which we will give explicit examples in section 3 and 4. As a consistency check on this setup, we will show that eq. (2.31) with these functions reproduces the (hard-)collinear region of one-loop (two-loop) diagrams.
Our approach provides important insight into the NLP behaviour of gauge theories. Firstly, we can explicitly test the categorisation of factorization ingredients given by the power-counting in eq. (2.23). Secondly, we shed light on some dynamical subtleties that are not fully accounted for by simpler factorization theorems as [29,30]. Finally, the explicit calculations presented here show how to deal with the endpoint contributions that result from a factorization structure consisting of convolutions rather than direct products.

Hard-collinear factorization for massive fermions
We now turn to the study of some of the ingredients entering the factorization picture, in the regime where the fermion mass is parametrically small. For definiteness we focus on the f γ-and f ∂γ-contributions to our NLP factorization formula, corresponding to N γ = 1 and N f = 1 in eq. (2.23). We calculate these jet functions at one-loop order in section 3.1 and validate them as well as the factorization structure through one-and two-loop calculations in section 3.2. The f γ-term is particularly relevant, according to our power counting formula, it already contributes at O(λ), for a parametrically small fermion mass. However, the f ∂γ-term is further suppressed by a power of the transverse momentum component. Thus, to appreciate the interplay between the two functions, we will need to carry out the calculation to O(λ 2 ). This level of accuracy, and the fact that we consider QED, constitutes an important generalisation of the analogous functions presented in [44] where Yukawa theory was studied. We stress that at O(λ 2 ) accuracy one also needs the fffand f γγ-jet functions, which contribute from two-loop onwards. We leave the calculation of these ingredients for future work. At this order other interesting aspects such as endpoint contributions come into play.
Although in this section we do not consider external soft radiation, our analysis of the non-radiative factorization ingredients is an important step towards generalising soft theorems to gauge theories at NLP, in the case of parametrically small masses. In addition, this scenario carries intrinsic interest to collider phenomenology, since precise measurements of cross sections may benefit from classifying and possibly resumming logarithms of small fermion masses at NLP. Examples are charm mass effects in B decays [86], initial-state mass effects in heavy-quark induced processes [87][88][89], bottom mass effects in Higgs production and decay [69,90], and tt production at a future linear collider, where the top mass could serve as a soft scale. Understanding the NLP factorization structure is a necessary intermediate step towards resumming such mass effects at this level of accuracy.

The massive f γ-jet
In the following we carry out an explicit derivation of the one-loop expressions for two of the jet functions that enter the NLP factorization formula for massive QED, as presented in eq. (2.31). The detailed calculation of these quantities sheds light on the subtleties involved beyond leading power. Moreover the functions we extract are process-independent, and could therefore be used in other QED calculations.
To keep expressions compact, we choose a reference frame such that the momentum p µ of the external particle that defines the jet has no perpendicular component p µ = (p + , 0, p − ), with p − p + . Unit vectors in the collinear and anti-collinear direction are then denoted by n µ = (1, 0, 0) andn µ = (0, 0, 1). From now on, we also set the electric charge q = −1. The power counting in eq. (2.23) was derived in axial gauge, which we also use for the jet functions. We will verify that, in a gauge which only allows for physical polarisations, the predicted power counting works on a diagram-by-diagram basis. For simplicity we select light-cone gauge, setting r 2 = 0 in eq. (2.8). Furthermore, we choose the reference vector in the anti-collinear direction, r µ = r −nµ . Note that r − then cancels in the photon propagator in eq. (2.8), leaving We extract the jet functions from the diagram in figure 7. Working in D = 4 − 2 dimensions, the corresponding amplitude is given by with the numerator factor and a generic n-loop hard function H (n) (f γ)ν . We first rearrange its transverse-momentum dependence by Taylor expanding in ρ ⊥ , as in eq. (2.34). 8 Retaining terms up to O(λ) 8 We recall that throughout this paper ρ ⊥ is the D-dimensional perpendicular component of ρ , as defined by means of a Sudakov decomposition: ρ = ·n n ρ + ·nn ρ + ρ ⊥ . Similarly we define η νρ ⊥ = η νρ −n νnρ −n ν n ρ .
we trade the initial hard function for two objects that depend only on the fraction of the large component of the loop momentum x = + /p + . Comparing eq. (3.2) with the first line of eq. (2.31), allows us to extract the jet functions, In eq. (3.5), we switched from the dominant loop momentum component + to the momentum fraction x, which determines the convolution between hard and jet functions. The x-integration range is a priori (−∞, +∞), but is in fact restricted to (0, 1) by noting that the integral over − vanishes if the two poles lie on the same side of the integration contour. In eq. (3.6) the denominators have homogeneous λ-scaling, but the numerator still needs expanding. As expected from the power counting rule for an all-collinear vertex, we find, in axial gauge, that N is O(λ); therefore, the f γ-jet starts at the same order, while the additional term ⊥ causes the f ∂γ-jet to begin at O(λ 2 ). Performing the expansion leaves us with three independent numerator structures, The first two lead to straightforward integrals, and follow from closing the integration contour at infinity in the − complex plane, evaluating the residue of the integrand at the pole − = − 2 ⊥ /(2 x p + ) − iη, and solving in turn the resulting integral over transverse momentum. The third one is more subtle, since the integrand does not vanish fast enough at the boundary to apply Jordan's lemma. Instead, we can isolate the troublesome term, introduce a Schwinger parameter, and integrate the minus component to a Dirac delta, This endpoint contribution at x = 1 corresponds to the limit where the photon leg carries all the momentum along the +-direction and the fermion line becomes soft.
Results for the integrals relevant to computing eq. (3.6) are collected in eq. (B.1) in appendix B. Having carried them out, we conclude These are the one-loop expressions for the f γ-and f ∂γ-jet functions in QED, as derived in light-cone gauge. As expected, the f γ-jet function starts at order λ ∼ m/Q, while the f ∂γ-jet has pure λ 2 scaling. This is due to the additional factor ρ ⊥ , which arises in the expansion of H (f γ) in eq. (3.4), which it absorbs in the definition of eq. (3.5). As a result only the structure α ⊥ β ⊥ in eq. (3.7) survives in the numerator. Since m is the only small scale in these functions, the mass expansion coincides with the power expansion.

Testing NLP factorization with the method of regions
Equipped with the result of eq. (3.9), we will now test the factorization formula eq. (2.31) in a process with two final-state jet directions, at both one-and two-loop order. Specifically, we wish to see whether this formula reproduces the (hard-)collinear limit of full, unfactorized amplitudes which at face value should be described by the f γ-and f ∂γ-jet functions. We will isolate the part of the amplitude that we want to compare with, using the method of regions [75,76,83,91]. This is a well-tested tool for expanding (loop) amplitudes in kinematic limits where the various scales entering the amplitude are largely separated in magnitude. It is particularly useful for the dissection of loop integrals, by defining regions where the virtual modes have momenta of a certain size as compared to a particular scale in the problem. In this case we use the small ratio of scales λ = m Q to select momentum regions where a virtual photon is hard, soft or collinear to either of the highly energetic particles in the final state. Once the regions have been defined, one may expand the integrand for each region in λ (up to an arbitrary order), which simplifies its structure. The integration is still carried out over the full momentum space, which allows for easy evaluation, but one must be careful not to overcount contributions that appear in multiple regions. Most of the time this causes no issue, as each region has a specific associated energy scale (in our case, (m 2 /µ 2 ) − for a collinear region and (2p + 1 p − 2 /µ 2 ) − for the hard region), which inhibits any cross-talk between such regions. Finally, by summing over all relevant regions, one obtains the result of the full integral up to the chosen order in λ.
In presence of just two jets in the final state, we choose a frame in which the jets are back to back. Given the light-cone decomposition of p µ 1 , we identify the direction collinear to p µ 2 as the anti-collinear direction. This yields the following regions Collinear : k µ ∼ Q 1, λ, λ 2 , Anti-collinear : k µ ∼ Q λ 2 , λ, 1 . In the following sections we refrain from considering all regions, but use this tool to extract only the contribution from the collinear (hard-collinear) region of the full amplitude, which is relevant for our one-loop (two-loop) test.

One-loop test
For one-loop accuracy we calculate the collinear region of figure 8, which should be described by contracting the one-loop f γ-and f ∂γ-jet functions with corresponding tree-level hard functions. We will carry out the regions calculation in axial gauge, as we did for the f γ-and f ∂γ-jet functions.
The full amplitude for the diagram in figure 8 reads The expansion of the integrand in the collinear region is obtained by rescaling the momentum components of both the (collinear) loop momentum and the (anti-)collinear external momenta, according to eq. (3.10). We further exploit our freedom of frame choice to set the perpendicular momentum components of the external momenta to zero. In practice, it is convenient to project onto a single set of light-like vectors in the plus-and minus-direction, using The denominator in eq. (3.11) is expanded as where all propagator denominators have now a homogenous λ-scaling. The numerator in eq. (3.11) is suppressed by one power of λ, allowing us to drop every term but the leading one from the denominator expansion in eq. (3.13), including the explicitly shown O(λ 2 ) term. By discarding higher power corrections in the numerator too, we readily calculate the collinear region up to NLP from this expression using standard techniques: we perform the Dirac algebra using the Mathematica package FeynCalc [92,93], Feynman parametrise the homogeneous denominators, shift the loop momentum and remove odd integrands, evaluate the momentum integrals through standard tensor integrals and finally integrate the Feynman parameters in a convenient order. We find (3.14) The subscript C indicates that this result is expanded in the collinear region. As expected, in axial gauge the diagram obeys the power counting, strictly contributing only at NLP. This expression has a single pole in , which receives both UV and IR contributions. The UV term regulates divergences that would be subtracted by one-loop renormalisation; the remainder has a collinear (rather than soft) origin.
The vector nature of the electromagnetic current and the Sudakov decomposition we employ in eq. (3.12) limit the possible Dirac structures that can appear in the result to γ α , n α , andn α . In particular, we observe that γ α occurs only in even powers of the mass expansion, while n α andn α multiply odd powers of the mass. In this specific case, the structuren α is absent due to a cancellation which, as our two-loop check will make clear, is accidental.
Turning to the factorization approach, we note that the collinear photon in figure 8 is emitted from an anti-collinear external line, such that the propagator before the emission carries a hard momentum. Consequently, this diagram should be described by the convolution of the one-loop f γ-and f ∂γ-jets with the respective (tree-level) hard functions, as claimed above. We expect the following factorization structure up to O(λ 2 ) The hard functions are extracted from figure 9, and its f γ-and f ∂γ-parts separated according to the prescription of eq. We emphasise that as we are interested in the first two orders in λ, we cannot ignore O(λ) terms in the numerator of eq. (3.16), as they will combine with the leading term in the f γ-jet (eq. (3.9)). In particular, we cannot drop the mass term. However, we can do so in eq. (3.17), since the f ∂γ-jet is proportional to two powers of the mass (thus O(λ 2 )). After some Dirac algebra, we find The integral over the energy fraction x is easily performed, yielding indeed the result in eq. (3.14). This provides a first check of our jet functions. We observe that the singular structure of the collinear region is entirely reproduced by the jet functions of eq. (3.18), while the convolution with the respective hard functions does not generate any additional pole in . We stress that the endpoint contribution, described by the Dirac delta function in eq. (3.18), is essential to obtain the correct result.

Two-loop test
We now proceed with a more strenuous test, based on the same method. The goal is to validate the factorization of a fermion-anti-fermion-production amplitude into f γ-and f ∂γ-jet functions if the hard function is loop-induced. A minimal diagram suited for this task is given in figure 10. It consists of an off-shell, one-loop vertex correction, described by the hard loop momentum 2 , which gets probed by a collinear fermion-photon pair forming the 1 -loop on the upper leg. We recall that focusing on one particular diagram is justified in axial gauge, where the power counting holds on a diagram-by-diagram basis and the factorization picture is derived. 9 Naturally, a complete evaluation of such a process would require us to determine the full hard function, necessitating the calculation of additional diagrams. We now proceed with the region expansion. The full two-loop amplitude reads  where for brevity we omitted the Feynman prescription iη in each of the square brackets in the denominators. The numerator structure reads To carry out the integrals in eq. (3.20) we use the same techniques as the one-loop example. The main difference is the presence of two-loop integrals, but due to the regions expansion the added complexity is limited. In the presence of masses and axial-gauge propagators, numerator structures proliferate, which makes the calculation computationally more intensive. However, as in the one-loop case, the numerator (3.21) scales as λ, which allows us to neglect O(λ 2 ) terms from the denominator expansion. In fact, only the denominator in eq. (3.20) mixing the two loop momenta generates O(λ) terms, through the expansion which is a consequence of our frame choice, p 1⊥ = p 2⊥ = 0. For the 1-loop-hard 1-loopcollinear (HC) region we thus obtain where Γ 1,2 denote the following combinations of Euler gamma functions Figure 11: Diagrammatic representation of the one-loop matrix element H (1) It is instructive to compare this result with its one-loop equivalent in eq. (3.14). Despite the more involved expressions for the coefficients, the basic structure is similar, but now all three different spin structures γ α , n α andn α contribute, with the latter arising only at order m 2 . There are other important differences, though. First, due to the more involved dynamical structure, the result features the two independent Γ-combinations in eq. (3.24). Second, since now a hard and a collinear loop are present at the same time, both scale ratios (m 2 /µ 2 ) − and (−2 p + 1 p − 2 /µ 2 ) − show up in the prefactor. However, setting for convenience the renormalization scale equal to the hard scale will remove the second factor. Upon expansion in this will result in logarithms of m 2 /(2p + 1 p − 2 ) at NLP, as for the oneloop case. These are small-mass logarithms that ideally would be resummed by a complete factorization framework.
We now continue with the calculation of the corresponding hard functions, and check that the convolution with the jet functions in eq. (3.9) reproduces our region calculation. Similar to the one-loop example, we can extract the hard functions by Taylor expanding the hard matrix element represented in figure 11, according to the prescription of eq. (3.4). The unexpanded amplitude reads from which we separate the f γ-and f ∂γ-term, Here the common denominator and the numerator structures are Note that the derivative in the transverse component defining the f ∂γ-term as in eq. (3.4) can act either on the spin structure in the numerator, or on the 1⊥ -dependent denominator of eq. (3.25), generating the two structures displayed in curly brackets. In the f ∂γnumerator structure we already dropped O(λ) terms, since we know that this structure enters the factorization formula in a convolution with a jet function that is already O(λ 2 ). We can now solve the integral with standard techniques. The presence of many different spin structures at this stage renders the intermediate expressions for the hard functions rather cumbersome, therefore we will not show them here. Taking the convolutions with the jet functions in eq. (3.9) yields the partial results shown in eq. (B.3) in the appendix. As one can readily verify, their sum correctly reproduces the hard-collinear region result obtained in eq. (3.23).
We will now examine the pole structure of eq. (3.23) in light of the equivalent factorization result. Interestingly, we note the presence of triple poles. One overall inverse power of is due to the single pole in the jet functions in eq. (3.9), while the remaining factor of 1/ 2 has two distinct origins. First, the hard functions contain explicit double poles since they describe both hard and soft physics. In section 4.3 we will extensively comment on this effect when examining the hard function in eq. (4.19). Second, the hard functions contain 1 1 x terms that produce an additional pole upon convolution with the respective jet functions. These endpoint singularities arise in the limit where the dominant momentum component of the collinear photon vanishes. Their origin is thus different from that of the (finite) endpoint contributions captured by δ(1 − x) in eq. (3.9), which describe the soft quark limit. Endpoint singularities appear in factorization studies using SCET, too [52,55,68,69], and seem inevitable at NLP. Since the expressions involved are unrenormalised quantities expressed in D = 4 − 2 dimensions, the endpoint singularities are easily regulated.

Hard-collinear factorization for massless fermions
In this section, we focus on the scenario of negligible fermion masses, m = 0, which is the standard approximation in high-energy collisions for light quarks. In the previous section we have seen that the fermion mass entered the non-radiative f γ-jet function through the overall scale factor m 2 /µ 2 − and a second-order polynomial in m. Removing this scale from the problem will thus have a serious impact on the ingredients in the factorization framework. Virtual loop corrections to the f -jet, as well as all loop-induced, genuine NLP jet functions (like the f γ-jet), are rendered scaleless and do not contribute. However as we are ultimately interested in threshold effects associated to soft final-state radiation, we are required to compute radiative jet functions. For such functions a new scale arises, set by the dot product of the (external) momenta of the emitting fermion and the soft photon. In massless radiative jet functions this small scale takes the place of the mass as the collinear scale.
As in the previous section, we will focus on the f γ-and f ∂γ-jet functions. The radiative functions will be obtained from the non-radiative counterparts by inserting a soft photon on any of the collinear fermion lines. Checking the factorization properties of gaugeinvariant sets of diagrams will allow us to make a convenient choice of gauge. While axial gauge proved to be practical for power counting the pinch surfaces that underlie the factorization ingredients, Feynman gauge is more suited for complex calculations. Therefore, we will extract the radiative jet functions here using the latter, and apply this gauge choice consistently in the calculation of the hard functions. The presence of longitudinally polarised photons in Feynman gauge will modify the power counting for individual diagrams. As a consequence, each diagram calculated in this section may contain spurious LP terms, which must cancel upon summing over a gauge-invariant set of diagrams.
The radiative jet functions we extract are process-independent quantities, which describe collinear physics regardless of the underlying hard scattering event. Similar to the massive case, we validate the expressions obtained by convolving these jets with appropriate hard functions by means of a one-and two-loop method of regions calculation. These non-trivial checks show that the all-order factorization formula for elastic amplitudes in eq. (2.31) provides a good starting point for the factorization (and potentially resummation) of threshold effects due to soft final state radiation.
In particular, our calculations show that the NLP factorization formula presented in [29,30] does not suffice for one-loop accuracy at the matrix element level, which has also recently been noted in a SCET context [52], although they do work at the cross section level for the cases studied there. Beyond one-loop such factorization formulae do not capture the intricate hard-collinear interplay for matrix elements that our current approach does account for.

The radiative, massless f γ-jet
At the lowest order in perturbation theory the radiative f γ-and f ∂γ-jet receive contributions from the two diagrams in figure 12. Both functions are defined in the same manner as in the massive fermion case and their evaluation relies on similar techniques. As before, we drop terms beyond NLP. In this case, we apply this constraint to the more involved denominator structure too, expanding denominators whose scaling is inhomogeneous in λ, as we did in the method of regions calculation of the amplitude. For example, after rescaling the momentum components by the appropriate powers of λ, the denominator of the innermost propagator is expanded as assuming collinear and soft scaling for and k respectively and abbreviating the homogeneous denominator by Figure 12: Contributions to the radiative f γ-and f ∂γ-jet.
Note that, having set m = 0, the external momentum can be chosen to be strictly in the +-direction p = (p + , 0, 0). Using this approach, we require six one-loop integrals to evaluate the contributions to the f γ-and f ∂γ-jet, which are listed in appendix C.1. We find the following result for these jet functions We point out that eq. . Note that eq. (4.3) does not contain the δ(1−x) term which appeared in the non-radiative jet function for the massive fermion case (eq. (3.9)) which was associated to the soft quark limit. In principle, one might expect a similar contribution here, but the numerator supplements the standard integrals of eqs. (C.1c) and (C.3c) with sufficient powers of (1−x) to suppress such a term. The radiative f -jet is known to have a Ward identity [29] relating it to its non-radiative counterpart (order by order in perturbation theory) via where q = −1 for the jets considered here. Similarly, we expect In particular, since J (f γ)ν and J (f ∂γ)νρ consist solely of scaleless integrals (for m = 0) and thus vanish in dimensional regularisation, we should find By contracting eq. (4.3) and eq. (4.4) with k µ , one finds that eq. (4.7) is indeed satisfied, which serves as a first check on these jet functions. Figure 13: Diagrams contributing to the collinear region of the 1R1V dijet production amplitude.

NLP factorization of the collinear sector at the one-loop level
Following the same approach as in 3.2, we wish to test the factorization structure of radiative amplitudes in the collinear sector, by means of a comparison to a method-of-regions computation of the single-real single-virtual (1R1V) correction to a dijet production process. A similar factorization/regions analysis has been carried out in refs. [30,77] for Drell-Yan production, at the same loop order but at the cross-section level instead. In [30] a NLP factorization formula for radiative amplitudes was derived from the standard LP factorization picture of purely virtual amplitudes, while here we start from the generalised NLP factorization formula of eq. (2.31). The difference between these two NLP approaches at the one-and two-loop level will be highlighted in the remainder of section 4. The diagrams that contribute to the collinear sector at the one-loop order are shown in figure 13. In diagrams (a) and (b) the collinear loop attaches only to the upper leg, meaning that there is only one fermion connection between the hard interaction and the part of the diagram containing the collinear dynamics. Therefore, these diagrams are predicted to factorize in terms of the one-loop radiative f -jet (see figure 14), the Born-level hard scattering amplitude (with amputated legs) and a trivial jet function for the opposite leg: M and H (0)α (f ) = −ieγ α . The one-loop radiative f -jet is readily computed by standard techniques and the result reads Note that this is strictly a NLP quantity, while from the non-radiative power counting formula (eq. (2.25)) one may have expected a contribution at LP. This power suppression is a radiative effect that only starts at one-loop order and is therefore not captured by the general power counting formula (the tree-level result does have a LP contribution). It is Figure 14: One-loop contributions to the radiative f -jet.
however fully consistent with [29], in which the radiative f -jet is defined to account for the NLP effects induced by soft emissions from collinear loops. As before, we may evaluate the contributions to the collinear sector of diagrams (a) and (b) in figure 13, using again the method of regions. Keeping terms up to NLP this approach yields For massless fermions we may choose p µ 1 = (p + 1 , 0, 0) and p µ 2 = (0, 0, p − 2 ), such that the standard (massless) Mandelstam variable t = (p 1 − k) 2 = −2p + k − . From eqs. (4.8) and (4.9) on the one hand and eq. (4.10) on the other, we see immediately that the regions result coincides with the factorization result which, given the trivial factorization structure of these diagrams, is perhaps not surprising.
For diagrams (c) and (d) in figure 13 we expect a factorization analogous to eq.
where s = (p 1 + p 2 ) 2 = 2p + 1 p − 2 . Since the f γ-and f ∂γ-jet function are strictly NLP quantities, we have discarded NLP corrections to both hard functions, as they would affect the full amplitude only at NNLP. Substituting eqs. (4.12) and (4.13) together with eqs. (4.3) and (4.4) into eq. (4.11) and simplifying the Dirac structure, we find p ℓ ℓ n k µ Figure 15: Additional contributions to the radiative f -jet, denoted by J (1)µ (f ) , in a simplified NLP factorization framework. Longitudinally polarised collinear photons that probe the hard scattering are described by the Wilson line interaction.
Since these missing terms vanish upon contraction with the conjugate amplitude, the simplified factorization approach did suffice in the 1R1V cross-section calculation presented in [30].

Hard-collinear interplay at the two-loop level
We now move to (single) radiative amplitudes at two-loop order (denoted as 1R2V) and carry out a similar test. At this loop order there is a more involved interplay between the hard and collinear sector, as the dominant component of the collinear momentum of the virtual photon may interfere with the hard loop. This hard-loop effect is not power suppressed, and has to be properly accounted for in the factorization picture in order to reproduce the exact NLP amplitude.
Our main effort here will be to explore this subtle interplay and therefore we (again) compare to a hard-collinear region with a method of regions calculation. The relevant diagrams for that purpose are shown in figure 16, where the collinear momentum is denoted by 1 and the hard momentum by 2 . We identify these diagrams through the following considerations.
First, the soft photon must originate from the collinearly enhanced region, rather than from the hard loop. Otherwise this would be described by a different term in the factorization formula, as stated by the Low-Burnett-Kroll theorem [27,28]: a soft finalstate emission from the hard scattering is described by a derivative with respect to either one of the external hard momenta, acting on the non-radiative hard scattering amplitude. 11 Second, the ordering of the virtual photon attachments is crucial. This is best seen from a Coleman-Norton analysis, in which hard, off-shell lines are shrunk to a point. In fact, it is strictly the attachment on the upper leg that matters, since the fermion propagators on the lower leg are shrunk to a point irrespective of the ordering. A propagator that is not part of the hard loop but which carries both an anti-collinear external momentum as well as a collinear loop momentum, obeys a hard scaling too. This implies that we can treat the planar-topology diagrams (c) and (d) in figure 16 as well as the crossed-topology diagrams (g) and (h) on equal footing. To see what happens if one inverts the order of attachments on the upper leg, let us consider diagram (c) as an example. In that case the outer loop Figure 16: Diagrams contributing to the hard-collinear region of the 1R2V dijet production amplitude, with 1 ( 2 ) denoting the collinear (hard) loop momentum. We use the shorthand notation 2 = 2 − 1 .
would be hard and therefore shrunk to the tree-level hard scattering vertex, as shown in figure 17. The supposedly collinear photon line would now form a tadpole-like attachment to the hard scattering vertex. However, this configuration cannot describe an on-shell line since it does not coincide with any classical trajectory [82], and does not contribute to the scattering amplitude. The collinear photon must thus attach to the upper leg outside of the hard loop, in order for the diagram to develop a hard-collinear region. Lastly, we note that these diagrams naturally contain a doubly collinear region too, which would be described by a higher-order radiative f γ-and corresponding f ∂γ-jet, as well as the radiative f γγ-jet, all contracted with a tree-level hard function. In a complete description of the doubly collinear region at this loop order, one would even expect contributions from the radiative fff -jet, which would be an interesting analysis by itself. This region does not overlap with the hard-collinear region we explore here, and thus we leave it to future work.
Analogously to the one-loop order, we foresee a pair-wise factorization of the diagrams in figure 16, by collecting those graphs differing only by the position of the radiated photon. For diagrams (a) and (b), we have with H (1)α (f ) (p 1 , p 2 ) the one-loop form factor. The hard function combined with the trivial jet function on the anti-collinear leg reads Figure 17: Coleman-Norton picture that arises from attaching the hard photon to the right of the collinear photon on the upper leg. The tadpole-like configuration of the supposedly collinear photon does not coincide with a classical trajectory. Hence this ordering of attachments does not contribute to the scattering amplitude.
Eq. (4.19) contains explicit double poles, while the unrenormalised hard function may only contain single poles of a UV nature; this double pole is thus of IR origin. In the method of regions, the appearance of IR poles in the hard region is a common phenomenon if the soft region is scaleless. (For the diagrams defined in figure 16 this is indeed the case, as is easily verified by assigning 2 a soft scaling according to eq. (3.10) and expanding denominators in λ.) Scaleless integrals are set to zero, which typically follows from a cancellation of IR and UV poles. Isolating this UV pole in the soft region and absorbing it in the hard region would cancel the double pole there, thus moving the double pole associated to soft physics from the hard to the soft region. We do not address this mixing of the hard and soft physics, as it affects the method-of-regions calculation and the hard function in the exact same way, while the collinear sectors, which are the focus of this study, do factorize entirely from the rest.
Turning to the remaining diagrams in figure 16, (c) to (h), we expect these to factorize according to with the one-loop f γ-and f ∂γ-hard functions extracted from figure 18. The calculation of these functions is deferred to appendix C.2 for conciseness. Upon evaluation of eqs. (4.18) and (4.20) we find Figure 18: Contributions to the one-loop hard functions H (f γ) (1) and H (f ∂γ) (1) .
We have combined the diagrams (c) to (h) rather than giving results per diagram pair, and have denoted combinations of gamma functions by the former coinciding with the second combination in eq. (3.24). The results of eq. (4.21) and eq. (4.22) are verified by calculating the hard-collinear region of the diagrams in figure 16. Given that the calculation is set up in a similar way as for the massive case, we will not provide further details for the sake of brevity. In particular, we find agreement between the factorization and regions results per diagram pair (a)+(b), (c)+(d), (e)+(f ) and (g)+(h). For the first three pairs we verified the exact agreement to all orders in , while for the last pair we compared series expansions in instead. This is due to the crossed topology of diagrams (g) and (h), which complicates the regions calculation by entangling Feynman parameters, yielding hypergeometric functions of the form 3 F 2 (a 1 , a 2 , a 3 ; b 1 , b 2 ; 1) upon integration. These multiply the second gamma function combination Γ 2 and are expanded up to and including finite terms (O 0 ) using HypExp [95,96]. By expanding the (exact) coefficients of the Γ 2 combination in the factorization result up to the same order, we verified their consistency.
The Γ 2 combination is in fact the signature of the mixing between hard and collinear loop momenta: starting at two loops, it originates from terms in the f γ-and f ∂γ-hard functions that carry an additional factor of x − , as seen in eq. (C.8). The appearance of those terms is, in turn, tied to an effective shift in the scale of the hard function; while the loop integration in H (f ) knows only the single scale 2 p + 1 p − 2 = s, the H (f γ) and H (f ∂γ) functions are sensitive to the dominant component of the collinear photon through 2 + 1 p − 2 = x s, giving additional x dependence. 12 Indeed, no Γ 2 combination is present for the f -jet factorization of diagrams (a)+(b) in eq. (4.21). As we see here, this effect is naturally captured by the NLP factorization formula of eq. (2.31). A simplified NLP factorization as in eq. (4.16), strictly in terms of f -jets, cannot do so: the complete factorization of the collinear and hard sector is an over-simplification of the intricate dynamics at play here.
Lastly, we point out that eqs. (4.21) and (4.22) contain at most 1/ 3 poles, while at two-loop order a maximal soft-collinear overlap would generate 1/ 4 poles. These leading singularities are captured by the soft function. 13 This means that for NLP threshold resummation purposes the collinear sector is needed starting at NLL accuracy. This has been noted before in the calculation of the collinear region of the 1R1V and 2R1V correction to Drell-Yan production in [77] and [78] respectively. Indeed, [34] showed that NLP threshold logarithms in Drell-Yan and single Higgs production are resummed at LL accuracy through an exponential next-to-soft function.

Conclusions
In this paper we have formulated a next-to-leading power factorization formula for njet production processes in QED, based on a power counting analysis for both zero and parametrically small fermion masses. We have thus begun the generalisation of the original Yukawa theory analysis of [44] to gauge theories. A factorization of degrees of freedom at NLP, following arguments such as in [10], could be an important step towards resummation of NLP (threshold) logarithms beyond leading-logarithmic accuracy.
We focused on the interaction that contributes at the first sub-leading power in λ, and computed the f γ-and f ∂γ-jets, which are universal quantities, up to order λ 2 . We first considered massive fermions, for which we calculated the non-radiative jet functions. To have a direct correspondence between the (next-to-)leading regions and their power counting we used axial gauge. We were able to test the factorization formula by comparing the convolution of jet functions and hard parts with a regions calculation of the two-jet amplitude to one loop. Subsequently we successfully tested these parts of the predicted factorization formula at the two-loop level by comparing the combined result against the hard-collinear region of a two-loop diagram. In particular, we pointed out the existence of two classes of endpoint contributions, one of which singular, that we deal with within dimensional regularisation.
For massless fermions we ensured the presence of a small scale, following [44], by adding an extra soft photon emission, so that the invariant of this photon momentum with a jet direction provides an analogue to the squared small fermion mass. We thus presented 12 The exact form of this x-dependence varies by diagram, as it is dictated by the denominators in the loop, and thus made explicit upon integration over the Feynman parameters that combine them. 13 Upon a correct assignment of poles, they would appear in the hard-hard region instead, by the same mechanism discussed before.
results for the radiative f γ-and f ∂γ-jet instead and tested the factorization of a radiative amplitude in a similar way as before. Here we found that the present approach reproduces all collinear contributions in the one-loop matrix element, contrary to the simplified NLP factorization of [29,30]. In addition, we noted a subtle interplay between hard and collinear modes, which is correctly accounted for in our factorization formula. We conjecture that our factorization formula is sufficiently general to factorize QED amplitudes up to NLP at arbitrary loop orders, and may thereby pave the way for the development of a similar factorization for QCD amplitudes.
We focused in our analysis on testing a specific part of the factorization framework. That is, the one that accounts for non-trivial hard-collinear interplay at the two-loop level, via novel jet functions which describe double hard-collinear interactions. In addition, similar tests should be carried out for the other ingredients (such as the triple hard-collinear interactions and the soft sector). Moreover, the jet functions used in this analysis are extracted from a generic jet-like scattering amplitude, rather than being derived from an operator definition. Such definitions would make the gauge invariance of the separate factorization ingredients manifest, thereby formalizing the NLP factorization framework for QED. These further steps, together with the extension to QCD, are part of ongoing work.
Acknowledgments EL acknowledges support from the Dutch NWO-I program 156, "Higgs as Probe and Portal". JSD and WW are supported by the D-ITP consortium, a program of NWO funded by the Dutch Ministry of Education, Culture and Science (OCW). LV acknowledges support from the Fellini -Fellowship for Innovation at INFN, funded by the European Unions Horizon 2020 research programme under the Marie Sk lodowska-Curie Cofund Action, grant agreement no. 754496. WW and LZ are supported by the ERC grant ERC-STG-2015-67732. This is also based upon work from COST Action CA16201 PARTICLEFACE supported by COST (European Cooperation in Science and Technology).

A Momentum regions for parametrically small fermion mass (m ∼ λQ)
The power-counting analysis in section 2 assumes the scaling k µ ∼ Q λ 2 , λ 2 , λ 2 , k µ ∼ Q 1, λ, λ 2 for soft and collinear momenta respectively. This follows from an analysis of the infrared structure of the scattering amplitude, which allows one to associate the pinch surfaces to momenta configurations that are soft and collinear. Here we complement this analysis by performing an expansion of the amplitude in momentum regions. This method provides an alternative approach for singling out the momentum configurations which are relevant for a given amplitude, in the presence of parametrically different scales, constituting a useful check for our assumptions in section 2. It also gives us the opportunity to briefly discuss differences between the two approaches.
To illustrate this second method, we focus on the scalar integral associated to the 1-loop diagram in figure 8. For fermions with mass m > 0 the integral reads e γ E Γ(1 + ) 2 m 2 μ 2 m 2 2 F 1 1, 1 + ,

B.1 Integrals for jet functions
For the calculation of the jet functions in the massive fermion case, we require with a common factor Note the different signs for the iη prescriptions in the denominators. As a consequence, the poles lie on opposite sides of the − integration contour if and only if 0 < x < 1, which restricts the convolution domain in x to that range.

B.2 Partial two-loop results
In the following we show expressions for the convolution of the jet and hard functions, that serve as a two-loop check of the result obtained in section 3.2.2. Here we list the f γ-and f ∂γ-terms separately: For diagram (b) in figure 12 a set of slightly more involved integrals is needed: with the f ∂s-function we derived. As for QED, we remark that a full treatment of the collinear sector at this order would require including f ss-and fff -jets, which however start contributing at two-loop order. Similar to section 3.2.2, we validated the factorization formula (D.3) using the method of regions. To this end, we expanded the two-loop diagram analogous to figure 10 in the hard-collinear region, where now photons are replaced by scalars, and verified that such a region is reproduced by the convolution between the jet presented in eq. (D.3) and the hard functions. We thus provided a check of the formalism of [44] beyond one loop and beyond O(λ).