Azimuthal asymmetries in QCD hard scattering: infrared safe but divergent

We consider high-mass systems of two or more particles that are produced by QCD hard scattering in hadronic collisions. We examine the azimuthal correlations between the system and one of its particles. We point out that the perturbative QCD computation of such azimuthal correlations and asymmetries can lead to divergent results at fixed perturbative orders. The fixed-order divergences affect basic (and infrared safe) quantities such as the total cross section at fixed (and arbitrary) values of the azimuthal-correlation angle φ. Examples of processes with fixed-order divergences are heavy-quark pair production, associated production of vector bosons and jets, dijet and diboson production. A noticeable exception is the production of high-mass lepton pairs through the Drell-Yan mechanism of quark-antiquark annihilation. However, even in the Drell-Yan process, fixed-order divergences arise in the computation of QED radiative corrections. We specify general conditions that produce the divergences by discussing their physical origin in fixed-order computations. We show lowest-order illustrative results for cos(nφ) asymmetries (with n = 1, 2, 4, 6) in top-quark pair production and associated production of a vector boson and a jet at the LHC. The divergences are removed by a proper all-order resummation procedure of the perturbative contributions. Resummation leads to azimuthal asymmetries that are finite and computable. We present first quantitative results of such a resummed computation for the cos(2φ) asymmetry in top-quark pair production at the LHC.


Introduction
Angular distributions of final-state particles produced by high-energy collisions are known to be relevant observables for the understanding of the underlying dynamics of their production mechanism. In this paper we deal with azimuthal-angle distributions and related asymmetries.
In the case of collisions that are produced by spin polarized particles, azimuthal angles with respect to the spin direction can be specified and measured: the analysis of the ensuing azimuthal correlations is a much studied topic and the subject of intense ongoing investigations (see, e.g., ref. [1] and references therein).
Azimuthal correlations in spin unpolarized collisions are less studied. A relevant exception is the lepton azimuthal distribution [2] of high invariant mass lepton pairs that are produced in hadron collisions through the Drell-Yan (DY) mechanism of quark-antiquark annihilation. Angular distributions for the DY process are a well known topic that has been deeply studied at both the theoretical and experimental levels (see, e.g., refs. [3][4][5] and references therein).
In the case of spin unpolarized collisions, much attention to azimuthal asymmetries has been devoted in recent years within the context of a QCD framework based on the introduction of non-perturbative transverse-momentum dependent (TMD) distributions of quarks and gluons inside the colliding unpolarized hadrons. In particular, the TMD JHEP06(2017)017 distribution of linearly-polarized gluons [6] inside unpolarized hadrons plays a distinctive role since it leads to specific cos(2ϕ) and cos(4ϕ) modulations of the dependence on the azimuthal-correlation angle ϕ. TMD distribution studies of such modulations have been performed, for instance, for the production process of heavy quark-antiquark (QQ) pairs [7] and for the associated production of a virtual photon and a jet (γ * + 1 jet) [8]. Many more related studies can be found in the list of references of ref. [8].
In this paper we consider spin unpolarized collisions and, specifically, we focus our discussion on hadron-hadron collisions, although many features that we discuss are generalizable to (and, hence, valid for) lepton-hadron and lepton-lepton collisions. We consider the production of a system F of two or more particles with a high value of the total invariant mass, and we examine the azimuthal correlations between the momenta of the system and of one of its particles. Roughly speaking, the relevant azimuthal angle ϕ is related to the difference between the azimuthal angles of the transverse momenta q T of the system and p T of the particle. Considering high values of the invariant mass M of the system F (additional kinematical cuts on the momentum of the produced particle can be applied or required), we select the hard-scattering regime of the production process. In this regime, azimuthal correlations are infrared and collinear safe observables [9], and they can be studied and computed by applying the standard QCD framework of factorization in terms of perturbative partonic cross sections and parton distribution functions (PDFs) of the colliding hadrons. Starting from the results in refs. [10,11] and, especially, from those in ref. [12] and the analysis on the role of soft wide-angle radiation that is performed therein, we develop and present some general (process-independent) considerations on azimuthal correlations and asymmetries.
We point out with some generality that, in spite of the infrared and collinear safety of azimuthal-correlation observables, their computation at fixed perturbative orders leads to divergences. The divergences appear also in basic quantities such as the total cross section at fixed (and arbitrary) values of the azimuthal-correlation angle ϕ. By contrast, the azimuthally-averaged (azimuthally-integrated) total cross section is known to be finite and computable order-by-order in perturbation theory. The finiteness of the azimuthallyaveraged cross section also implies that the divergences are associated with non-trivial modulations of the ϕ dependence.
Examples of processes with fixed order (f.o.) divergences are heavy-quark pair (QQ) production, associated production of vector (V ) or Higgs bosons and jets (e.g., 'V + 1 jet' with V = Z, W ± , γ * ), dijet and diboson production. We support our theoretical discussion by performing the numerical calculation of cos(nϕ) modulations (for some values of n) at the lowest perturbative order for some specific processes. The numerical results are found to be consistent with the divergent behaviour if n = 2, 4, 6 in the case of top-quark pair (tt) production and if n = 1, 2, 4, 6 in the case of associated 'Z + 1 jet' production. A noticeable exception to the appearance of QCD divergences is the case of the DY process. However, even in the case of the DY process, we predict the presence of f.o. divergences in the perturbative computation of QED radiative corrections. Moreover, analogous divergences arise in the case of lepton-lepton collisions from the computation of radiative corrections in a pure QED context (e.g., in the QED inclusive-production process e + e − → µ + µ − + X, JHEP06(2017)017 or even in the simpler process γγ → µ + µ − + X). In summary, processes that at the leading order (LO) are produced with an exactly vanishing value, q T = 0, of the transverse momentum q T of the system F tend to develop sooner or later (in the computation of QCD and QED radiative corrections at subsequent perturbative orders) f.o. divergences.
A possible reaction to this state of affairs is to advocate non-perturbative stronginteractions dynamics and related non-perturbative QCD effects that can cure the f.o. divergences. This cure would spoil the common wisdom according to which non-perturbative contributions to infrared and collinear safe observables are (expected to be) power suppressed, namely, suppressed by some inverse power of the relevant hard-scattering scale (as set by the high invariant mass M of the system F ). Moreover, non-perturbative strong-interactions dynamics cannot be advocated in the case of f.o. divergences in a pure QED context.
We pursue a more conventional viewpoint within perturbation theory. We examine the origin of the singularities for azimuthal-correlation observables order by order in perturbation theory, and we proceed to resum the singular terms to all orders [11,12].
The singularities originate from the kinematical region where q T → 0. We identify two sources of singularities. One source [10,11] is initial-state collinear radiation in the case of systems F that can be produced with vanishing q T by gluon initiated partonic collisions. The other source [12] is soft wide-angle radiation in the case of systems F that contain colour-charged particles (or electrically-charged particles for QED radiative corrections). Collinear radiation is, per se, responsible for singularities in cos(nϕ) (and sin(nϕ)) harmonics with n = 2 and n = 4, only [11]. Soft radiation is responsible for singularities that, in principle, can affect harmonics with any (both even and odd ) values of n (n = 1, 2, 3, 4 . . . ).
At high perturbative orders collinear and soft radiation is mixed up, and the singularities are enhanced by logarithmic (ln(M/q T )) contributions. We discuss in general terms how the all-order resummation of the logarithmic contributions leads to (q T integrated) azimuthal asymmetries that are finite and computable. Such procedure is effective in both QCD and QED contexts. In the QCD case this does not imply that non-perturbative effects have a negligible quantitative role in the region of very-low values of q T , but these effects give small contributions to q T -integrated azimuthal-correlation observables. Within perturbative QCD, the most advanced theoretical treatment of singular azimuthal correlations that is available at present regards the process of heavy-quark pair production [12]. We use the process of top-quark pair production in pp collisions at the LHC as a prototype to present a first illustration of the resummation of singular azimuthal asymmetries at the quantitative level.
A correspondence can be established [10,11,13] between the TMD distribution of linearly-polarized gluons and singular azimuthal correlations that are due to initial-state collinear radiation in gluon initiated partonic subprocesses. Therefore, both the framework used in refs. [7,8] and our perturbative QCD treatment lead to corresponding cos(2ϕ) and cos(4ϕ) azimuthal modulations. In this respect, the TMD distribution of linearly-polarized gluons can be regarded as the extension of specific features of QCD collinear dynamics from perturbative to non-perturbative transverse-momentum scales. Much discussion of the JHEP06(2017)017 present paper is related to soft radiation effects that occur in processes with colour-charged final-state particles. Soft radiation produces azimuthal asymmetries (which are singular at f.o. and finite after resummation) for both quark (or antiquark) and gluon initiated subprocesses, and for cos(nϕ) harmonics with arbitrary even values (n = 2, 4, 6, . . . ) [12] and also odd values of n. These soft wide-angle radiation effects are unrelated to the TMD distribution of linearly-polarized gluons.
The paper is organized as follows. In section 2 we start our discussion on azimuthal asymmetries and we specify the conditions that lead to divergences in f.o. computations. In section 3 we discuss azimuthal asymmetries in two examples of hadron collider processes, i.e. the production of lepton pairs through the DY mechanism and the production of a tt pair, by contrasting the different behaviour of the corresponding azimuthal harmonics. In section 4 we start our analysis of the small-q T limit: in section 4.1 we focus on the smallq T behavior at f.o. and in section 4.2 we recall the transverse-momentum resummation procedure for the case of azimuthally-averaged cross sections. In section 5.1 we discuss the origin of singular azimuthal correlations and present illustrative lowest-order results for the Z+jet production process. In section 5.2 we outline the resummation procedure of singular terms in the case of azimuthally-correlated cross sections and we contrast the small-q T behavior expected for the resummed cross section with the known behavior of the azimuthally-averaged transverse-momentum cross section. In section 5.3 we focus on the n = 2 harmonic for tt production and we present first quantitative results of a resummed calculation at next-to-leading logarithmic accuracy. After the matching with the complete next-to-leading order (NLO) result, the resummed computation offers an effective 'lowestorder' prediction for the n = 2 harmonic. We also comment on the possible role of nonsingular terms. In section 6 we summarize our results.

Azimuthal correlations and asymmetries in fixed-order perturbation theory
Our discussion on azimuthal correlations has a high generality. To simplify the illustration of the key points, we consider the simplest class of processes, in which the produced highmass system in the final state is formed by only two 'particles' in generalized sense (pointlike particles and/or jets). We consider the inclusive hard-scattering hadroproduction process in which the collision of the two hadrons h 1 and h 2 with momenta P 1 and P 2 produces the triggered final state F , and X denotes the accompanying final-state radiation. The observed final state F is a generic system that is formed by two 'particles', f 3 and f 4 , with four momenta p 3 and p 4 , respectively. The two particles can be point-like particles or hadronic jets (j), which are reconstructed by a suitable (infrared and collinear safe) jet algorithm. As for the case of point-like particles, the most topical process is the production of a high-mass lepton pair through the DY mechanism of quark-antiquark annihilation. We consider many other cases such as, for instance, the production of a photon pair (γγ), JHEP06(2017)017 a pair of top quark and antiquark (tt) or a pair of vector bosons (V V ), in addition to dijet production (jj) and associated production processes such as vector boson plus jet (V j). The invariant masses of the two particles have a little role in the context of our discussion (and they do not affect any conceptual aspects of our discussion). The system F has total invariant mass M (M 2 = (p 3 + p 4 ) 2 ), transverse momentum q T and rapidity y (transverse momenta and rapidities are defined in the centre-of-mass frame of the colliding hadrons). We require that M is large (M Λ QCD , Λ QCD being the QCD scale), so that the process in eq. (2.1) can be treated within the customary perturbative QCD framework. We use √ s to denote the centre-of-mass energy of the colliding hadrons, which are treated in the massless approximation (s = (P 1 + P 2 ) 2 = 2P 1 · P 2 ).
The dynamics of the production process in eq. (2.1) can be described in terms of five kinematical variables: the total mass M , transverse momentum q T and rapidity y of the system F and two independent angular variables that specify the kinematics of the two particles f 3 and f 4 with respect to the total momentum q = p 3 +p 4 of F . These two angular variables are a polar-angle variable and an azimuthal-angle variable. To be definite and to avoid the use of 'exotic' variables, we refer to a widely-used set of angular variables and we use the polar angle θ and the azimuthal angle ϕ (of particle f 3 ) in the Collins-Soper (CS) rest frame 1 [2] of the system F .
The variable ϕ is the relevant variable for our discussion of azimuthal correlations. We remark that ϕ specifies the azimuth of one of the two particles in the system F with respect to the total momentum of the system. In particular, we also remark that we are not considering the relative azimuthal separation ∆φ = φ 3 − φ 4 (φ i = φ(p T i ), with i = 3, 4, is the azimuthal angle of the transverse-momentum vector p T i in the centre-of-mass frame of the colliding hadrons) between the two particles. However, we can anticipate (we postpone comments on this point) that our main findings are not specific of the CS frame, and they are equally valid for other azimuthal variables with respect to the system F (for instance, we can consider the azimuthal angle in a different rest frame of F or, simply, the azimuthal difference φ(p T i ) − φ(q T ), where φ(q T ) is the azimuthal angle of q T in the centre-of-mass frame of the colliding hadrons).
Using the kinematical variables in the CS frame, we can consider azimuthal distributions for the process in eq. (2.1). The most elementary azimuthal-dependent observable is the azimuthal cross section at fixed invariant mass, dσ and we can also consider less inclusive observables such as, for instance, the q T -dependent azimuthal cross section in eq. (2.3) and the multidifferential (five-fold) cross section in eq. (2.4):

JHEP06(2017)017
All these quantities are related (from the less inclusive to the more inclusive case) through integration of kinematical variables (for instance, the azimuthal cross section in eq. (2.2) is obtained by integrating eq. (2.3) over q 2 T ) and, in particular, the azimuthal integration of eq. (2.2) gives the total cross section (at fixed invariant mass) of the process: Obviously, we can also consider differential cross sections that are integrated over a certain range of values of the invariant mass. Since all the cross sections that we have just mentioned are infrared and collinear safe quantities [9], they can be computed perturbatively within the customary QCD factorization framework (see ref. [14] and references therein). The only non-perturbative (strictly speaking) input is the set of PDFs of the colliding hadrons. The PDFs are convoluted with corresponding partonic differential cross sections dσ that can be evaluated as a power series expansion in the QCD running coupling α S (M ).
This perturbative QCD framework is applicable at any finite (and arbitrary) fixed perturbative order. Despite this statement, the first main observation that we want to make is that the f.o. perturbative calculation of the azimuthal distributions can lead to divergent (and, hence, unphysical and useless) results. More specifically, the f.o. calculation of the azimuthal cross section of eq. (2.2) gives the following results: finite at any f.o. (DY production) divergent for any ϕ at some f.o. (tt, V j, jj, γγ, ZZ, W + W − , . . . , production) . (2.6) We mean that the perturbative computation of dσ/dM 2 dϕ for the DY process gives a finite result order-by-order in QCD perturbation theory, while in most of the other cases (some of them are listed in the right-hand side of eq. (2.6)) the computation gives a divergent (meaningless) result for any values of the azimuthal angle ϕ starting from some perturbative order. We note that the integration over ϕ of dσ/dM 2 dϕ gives the total cross section (see eq. (2.5)), which is known to be finite at any f.o.. Therefore, the divergent behaviour that is highlighted in eq. (2.6) originates from the azimuthal-correlation 2 component, dσ corr /dM 2 dϕ, of the azimuthal cross section: where the notation · · · av. denotes the azimuthal average and dσ/dM 2 is the total cross section in eq. (2.5).
To understand the origin of the divergent behaviour in eq. (2.6), we first comment about kinematics. If the system F has vanishing transverse momentum (q T = 0), the rest JHEP06(2017)017 frame of F is obtained by simply applying a longitudinal boost to the centre-of-mass frame of the colliding hadrons. Any additional rotation in the transverse plane of the collision leaves the system F at rest and makes the particle azimuthal angle ϕ ambiguously defined. Owing to the azimuthal symmetry of the collision process, if q T = 0 there is no preferred direction to define ϕ. In other words, considering the azimuthal angle ϕ (as defined in the CS frame or any other rest frame of F ) and performing the limit q T → 0, we have where φ 3 = φ(p T 3 ) and φ(q T ) are the azimuthal angles of the corresponding transversemomentum vectors in the centre-of-mass frame of the colliding hadrons. If q T = 0, φ(q T ) is not defined and, consequently, ϕ is not (unambiguously) defined. At the strictly formal level, to define the azimuthal correlation we have to exclude the phase space point at q T = 0. If we want to consider q T integrated azimuthal correlations (such as, for instance, the cross section in eq. (2.7)), we have to introduce a minimum value of q T (q T > q cut ) and eventually perform the limit q cut → 0. The divergences that are highlighted in eq. (2.6) do not appear by considering azimuthal correlations at fixed and finite values of q T (e.g., the azimuthal-dependent cross sections in eq. (2.3) and (2.4)). These divergences are related to the limit q cut → 0 after integration of the azimuthal correlations over q T . Our discussion on the limit q cut → 0 reconciles the appearance of divergences in eq. (2.6) with the perturbative criterion of infrared and collinear safety at the formal level: no divergence occurs provided q T = 0, namely, provided the azimuthal-correlation observable is well (unambiguously) defined. However, a single phase space point at q T = 0 is physically harmless. Therefore, the q T integrated azimuthal correlations (i.e., their limiting behaviour as q cut → 0) are finite and measurable quantities. Moreover, they are also finite and measurable if q cut is non-vanishing, or, equivalently if q T is not vanishing and fixed at an arbitrarily small value. The divergence in eq. (2.6) implies that the q T dependent azimuthal correlations become singular at small values of q T , and that this singular behaviour is not integrable over q T in the limit q T → 0. This unphysical behaviour of f.o. perturbative QCD at small values of q T and the divergence of the q T integrated azimuthal correlations certainly require some deeper understanding to have a QCD theory of physically measurable azimuthal correlations. As a possible shortcut, one can still use f.o. perturbative QCD but avoid the region of small values of q T . In this case one still needs some understanding of the phenomenon in order to assess the extent of the dangerous small-q T region where the f.o. predictions are 'unphysical' or, anyhow, not reliable at the quantitative level.
We note that the kinematical relation in eq. (2.8) implies that the azimuthal angle in the CS frame has no privileged role in the context of our discussion of azimuthal correlations in the small-q T region. The main features of the small-q T behaviour of azimuthal correlations that are discussed in this paper are basically unaffected by using azimuthal angles as defined in other rest frames of the system F , or, by using other related definitions of azimuthal angles. For instance, one can simply replace the CS frame angle ϕ with one of the relative azimuthal angles φ(p T i ) − φ(q T ) (i = 3, 4) in the centre-of-mass frame of the colliding hadrons. Alternatively, one can use the difference between the azimuthal angles (in the JHEP06(2017)017 centre-of-mass frame of the colliding hadrons) of the two transverse-momentum vectors p T 3 − p T 4 and q T . All these definitions of the relevant (for our purposes) azimuthal-angle variable turn out to be equivalent (because of eq. (2.8)) in the limit q T → 0 (or, at very small values of q cut ). In the following we continue to mainly refer ourselves to the azimuthal angle in the CS frame, although all the basic features that we discuss are unchanged by considering the other definitions of the azimuthal angle.
We also note that the discussion that we have presented so far about azimuthal correlations can be generalized in a straightforward way to consider the case in which the system F is formed by more than two particles. In the multiparticle case, we can simply examine azimuthal correlations that are defined by using the azimuthal angles of the various particles in F in a specified rest frame of F (such as the CS frame), by simply using the various relative azimuthal angles φ(p T i )−φ(q T ) in the centre-of-mass frame of the colliding hadrons, or by using some other related (i.e., equivalent in the limit q T → 0) azimuthal variables. Since the multiparticle cases do not present any additional conceptual issues with respect to the two-particle case, in the following (for the sake of technical simplicity) we continue to mostly consider the case of systems F that are formed by two particles.
The f.o. perturbative divergences of the azimuthal correlations originate from QCD radiative corrections due to inclusive emission in the final state of partons with low transverse momentum (soft and collinear partons). The origin is discussed in more detail in the following sections of the paper (see, in particular, section 5.1). Since the azimuthal correlations behave differently in different processes (as stated in eq. (2.6)), we anticipate here the conditions that produce the divergent behaviour.
Azimuthal correlations can have f.o. divergences if the final-state system F ({p 3 , p 4 }) can be produced by the partonic subprocess c 1 c 2 → F (see also eq. (4.3) and accompanying comments) where • at least one of the initial-state colliding partons c 1 and c 2 is a gluon; (2.9) • at least one of the final-state particles with momenta p 3 and p 4 carries colour charge. (2.10) The conditions in (2.9) and (2.10) follow from a generalization of the results in refs. [11,12]. The divergences arise from the computation of the QCD radiative correction to the partonic process c 1 c 2 → F (they do not arise in the computation of that partonic subprocess itself!). Specifically, the f.o. divergences originate from collinear-parton radiation [11] in the case of the condition (2.9) and from soft-parton radiation [12] in the case of the condition (2.10). We remark that one of the conditions in (2.9) and (2.10) is sufficient to produce the f.o. divergences. In particular, the condition (2.10) necessarily produces f.o. divergences, while the condition (2.9) produces divergences with some 'exception' (words of caution) in few specific cases (see below). Having discussed the source of f.o. divergences in general terms, we can comment on some specific processes.
The production of heavy-quark pairs such as, for instance, tt can occur through the partonic subprocesses of quark-antiquark annihilation (qq → tt) and gluon fusion (gg → tt). One of these subprocesses has initial-state gluons and, moreover, both t andt have QCD colour charge. Therefore, if F = {tt} both conditions (2.9) and (2.10) are fulfilled and the JHEP06(2017)017 corresponding azimuthal correlations have f.o. divergences (as stated in eq. (2.6)). The same reasoning and conclusions apply to the production of F = {V j} and F = {jj} (for instance, F = {V j} can be produced through qg → V j, where at the lowest-order j = q carries QCD colour charge). In the case of F = {γγ} production the final-state photons do not carry QCD colour charge. However, diphoton production at the next-to-next-to-leading order (NNLO) can occur through the subprocess gg → γγ (the interaction is mediated by a quark loop): therefore, due to the condition (2.9), the azimuthal correlations for diphoton production diverge starting from the N 3 LO computation. The cases F = {ZZ} and F = {W + W − } behave similarly to F = {γγ}, and they also lead to azimuthal correlations that have f.o. divergences.
We consider the DY process, where F = { } and the high-mass lepton pair originates from the decay of a vector boson V (V → ). The final-state leptons have no QCD colour charge and, therefore, the condition (2.10) is not fulfilled. The partonic subprocesses qg → V ( ) andqg → V ( ) are forbidden by colour conservation, and the partonic subprocess gg → V ( ) is forbidden by the spin 1 of the vector boson. The DY lepton pair can be produced through qq → V ( ), but this partonic subprocess has no initial-state gluons. Therefore, also the condition (2.9) is not fulfilled and the azimuthal correlation for the DY process have no f.o. divergences in the computation of QCD radiative corrections (as it is well known and recalled in eq. (2.6)). However, we remark that f.o. divergences do appear in the computation of QED (or, generally, electroweak) radiative corrections to azimuthal correlations for the DY process (this important point is discussed in more detail in section 3).
We can comment on the production of a system of colourless particles, such as F = {γγ} or F = {ZZ}, in the specific case (or, better, within the approximation) in which those particles originate from the decay of a spin-0 boson, such as the Standard Model (SM) Higgs boson H (e.g., H → γγ or H → ZZ). In this case the production mechanism gg → H (followed by the H decay) is allowed. However, due to the spin-0 nature of H, the H decay dynamically decouples from the H production mechanism: the angular distribution of the final-state particles (γγ or ZZ) with momenta p 3 and p 4 is dynamically flat (it has no dynamical dependence on the decay angles, since it can only depend on the Lorentz invariant 2p 3 · p 4 = M 2 , namely, on the invariant mass M of the produced pair of particles). As a consequence, although the condition (2.9) is fulfilled, there are no azimuthal correlations in this specific case and, hence, there are no accompanying f.o. divergences. The absence of azimuthal correlations follows from the requirement that the two final-state particles are due to the decay of a spin-0 boson. As a matter of principle at the conceptual level we note that, considering the final-state system F = {γγ} or F = {ZZ}, the various subprocesses that contribute to the production mechanism gg → F (subprocesses with and without an intermediate H, and corresponding interferences) are not physically distinguishable (this is, strictly speaking, correct although the applied kinematical cuts can sizeably affect the relative size of the various contributing subprocesses). As we have previously discussed, the production of such systems without the intermediate decay of a spin-0 boson leads to azimuthal correlations with f.o. perturbative divergences. Therefore, if the condition (2.9) is fulfilled, we can conclude that sooner or later (in the computation of subsequent perturbative orders) QCD radiative corrections produce non-vanishing azimuthal correlations and ensuing f.o. divergences.

JHEP06(2017)017
To the purpose of studying azimuthal correlations and presenting corresponding quantitative results, we find it convenient to introduce harmonic components of azimuthaldependent cross sections. In particular, we define the n-th harmonic: (2.11) where n is a positive integer (n = 1, 2, 3, . . . ). Note that, in view of our previous discussion on the origin of the divergences that are mentioned in eq. (2.6), we have introduced a minimum value q cut of q T (q T > q cut ). We usually consider values of q cut that are proportional to the invariant mass M of the system F (q cut = r cut M , with a fixed parameter r cut ), although also fixed values of q cut can be used. One can also consider n-th harmonics of more differential cross sections (e.g., differential with respect to y and cos θ as in eq. (2.4)). The n-th harmonic in eq. (2.11) is defined by using the weight function cos(nϕ); n-th harmonics with respect to the sine function can also be considered by simply replacing cos(nϕ) → sin(nϕ) in the integrand of eq. (2.11). In particular, the knowledge of the harmonics with respect to both cos(nϕ) and sin(nϕ) for all integer values of n (n = 1, 2, 3, . . . ) is equivalent to the complete knowledge of the azimuthal-correlation cross section in eq. (2.7). Note that the n-th harmonic in eq. (2.11) is not a positive definite quantity. Although the physical azimuthal cross section dσ/dM 2 dq 2 T dϕ is positive definite, the weight function cos(nϕ) (or sin(nϕ)) has no definite sign.
The n-th harmonic with weight cos(nϕ) (or sin(nϕ)) gives a direct measurement of the cos(nϕ) (or sin(nϕ)) asymmetry of the azimuthal-dependent cross section. The QCD computation of the n-th harmonic gives a finite result provided q cut is not vanishing. On the basis of eq. (2.6), in the limit q cut → 0 some harmonics (asymmetries) can be divergent if computed at some f.o. in QCD perturbation theory. We cannot draw general conclusions on harmonics with odd values of n (n = 1, 3, 5, . . . ), but from the results in refs. [11,12] we know that harmonics with even values of n (n = 2, 4, 6, . . . ) can have f.o. divergences. Specifically, if the condition (2.9) is fulfilled the harmonics with n = 2 and n = 4 have f.o. divergences [11], and if the condition (2.10) is fulfilled all the n-even (n = 2, 4, 6, 8 and so forth) cos(nϕ) asymmetries have f.o. divergences [12]. In section 3 we explicitly show quantitative results on f.o. computations of harmonics for the DY and tt production processes. In section 5.1 we also show f.o. results for V j production and we further comments on n-odd harmonics.
As mentioned in the Introduction and discussed in section 5.2, the f.o. divergences of the azimuthal asymmetries can be cured by a proper all-order perturbative resummation of QCD radiative corrections. The resummed computation leads to azimuthal correlations and asymmetries that are finite, as it is the case of the corresponding physical (and measurable) quantities.

Examples
The lepton angular distribution of the DY process is a much studied subject both experimentally and theoretically. Recent measurements of the lepton angular distribution for JHEP06(2017)017 Z production at LHC energies ( √ s = 8 TeV) are presented in refs. [3,4], together with corresponding QCD predictions from Monte Carlo event generators and f.o. calculations. A very recent phenomenological study of the DY lepton angular distribution is performed in ref. [5]. As for previous literature on the subject, we mainly refer the reader to the list of references in refs. [3][4][5].
The QCD structure of the DY lepton angular distribution is well known, and it is a consequence of the spin-1 nature of the vector bosons (V = Z, γ * , W ± ) involved in the DY mechanism. To be definite we consider the DY process at the Born level with respect to electroweak (EW) interactions. Within this Born level framework, the DY differential cross section is expressed in terms of a leptonic tensor (for the EW leptonic decay of the vector boson) and a hadronic tensor (for the hadroproduction of the vector boson). The QCD production dynamics is embodied in the hadronic tensor and the two tensors are coupled through the spin (helicity) correlations of the vector boson. Owing to quantum interferences, we are dealing with 9 helicity components of the five-fold differential cross section in eq. (2.4) (the hadronic tensor is a 3 × 3 helicity polarization matrix, due to the 3 polarization states of the vector boson), which, therefore, can be expressed [15] as a linear combination of 9 spherical harmonics or, better, harmonic polynomials (the leptonic tensor has a polynomial dependence of rank 2 on the lepton momenta) of the lepton angles θ and ϕ (we specifically use the CS frame). For the illustrative purposes of our discussion of azimuthal correlations, it is sufficient to consider the four-fold differential cross section dσ/(dM 2 dy dq 2 T dϕ), which is obtained by integrating eq. (2.4) over cos θ. We write where dσ/(dM 2 dy dq 2 T ) is the DY differential cross section integrated over ϕ, and the shorthand notation [dσ] I denotes differential cross section components that only depend on M, y, q T . The azimuthal dependence of eq. (3.1) is entirely given by the harmonics {cos ϕ, cos(2ϕ), sin ϕ, sin(2ϕ)} that are explicitly written in the right-hand side of eq. (3.1). The subscript I in [dσ] I is defined according to a customary notation in the literature [3,4,15,16], such that [dσ] I ≡ C I A I (M, y, q T ) × dσ/(dM 2 dy dq 2 T ), where the functions A I (M, y, q T ) (I = 0, 1, . . . , 7) are known as 'DY angular coefficients' and C I are normalization factors (specifically, we have C 2 = 1/4, C 5 = 1/2, C 3 = C 7 = 3π/16). A point that we would like to remark is that the QCD dependence of the azimuthal correlations of the DY cross section involves only four harmonics. The five-fold cross section in eq. (2.4) for the DY process is expressed in terms of 9 harmonic polynomials of θ and ϕ: however, their dependence on ϕ is given only in terms of the four harmonics that also appear in eq. (3.1).
At LO in QCD perturbation theory, the DY cross section is due to the partonic subprocess qq → V (→ ), which is of O(α 0 S ). At this order in α S , there is no azimuthal dependence. Azimuthal correlations start to appear at the NLO 3 through the tree-level partonic JHEP06(2017)017 processes at O(α S ) (qq → V + g and its crossing related channels, such as qg → V + q). At this order only the cross section components [dσ] 3 and [dσ] 2 in eq. (3.1) are not vanishing. The azimuthal harmonics sin ϕ and sin(2ϕ) receive non-vanishing contributions only starting from NNLO processes. More precisely, these non-vanishing contributions [17,18] are entirely due to one-loop absorptive (and time reversal odd) corrections to the partonic subprocess qq → V + g and its crossing related channels. Azimuthal correlations up to NNLO were first computed in ref. [15]. Azimuthal-correlation results at N 3 LO can be obtained, in principle, by exploiting recent progress on the computation of the NNLO corrections to 'V + 1 jet' production [19][20][21][22][23].
We recall that the cross section components [dσ] 2 and [dσ] 3 in eq. (3.1) are parity conserving, while [dσ] 5 and [dσ] 7 are parity violating. Moreover, [dσ] 3 vanishes in the limit of vanishing axial coupling of V to either quarks or leptons. Therefore, if V = γ * , only the cos(2ϕ) harmonic contributes in eq. (3.1), while in the case of the five-fold differential cross section of eq. (2.4) there is an additional non-vanishing azimuthal contribution that is proportional to sin(2θ) cos ϕ and it is due to [dσ] 1 .
All the quantitative results that are presented in this paper refer to pp collisions at the LHC energy √ s = 8 TeV. In the QCD calculations we use the set MSTW2008 [24] of PDFs at NLO.
To present some illustrative results for the DY process we consider on-shell Z production and its leptonic decay Z → e + e − in an electron-positron pair. Our QCD calculation is performed by using the numerical Monte Carlo program DYNNLO [25]. The EW parameters are specified in the G µ scheme and we use the following values: the Fermi constant is G F = 1.1663787 · 10 −5 GeV −2 , the mass of the Z boson is M Z = 91.1876 GeV and the mass of the W boson is M W = 80.385 GeV. We use equal values of the factorization (µ F ) and renormalization (µ R ) scales and we set µ F = µ R = M Z .
We specifically consider the numerical calculation of various cos(nϕ) asymmetries (see eq. (2.11)), and ϕ is the azimuthal angle of the electron in the CS frame. We evaluate the asymmetries at their non-trivial lowest order in f.o. perturbation theory and, therefore, we compute all the NLO tree-level partonic processes whose final state is 'Z(e + e − )+1 parton'. Our numerical results are presented in figure 1-left.
The n-th harmonics are integrated over the lepton polar angle and over the rapidity y and transverse-momentum q T of the e + e − pair, and they are computed as a function of r cut = q cut /M Z , where q cut is the minimum value of q T (q T > q cut ). The results in figure 1 are obtained for very small values of r cut in the range 5 · 10 −4 < r cut < 5 · 10 −3 . As we have already recalled, the azimuthal asymmetries for the DY process have no f.o. divergences. Indeed, the results for n = 1, 2, 4 that are presented in figure 1-left show that the corresponding cos(nϕ) asymmetries are basically independent of r cut for very small values of r cut and finite (the results in figure 1-left practically coincide with the numerical evaluation at r cut = 0). The result for the n = 4 asymmetry (blue line) in figure 1-left is consistent with a vanishing value, in agreement with the general expression in eq. (3.1) specific process. In the case of azimuthal correlations, the N k LO result can effectively represent a QCD prediction at some lower perturbative order. (the very small deviations that are observed for n = 4 in figure 1-left give an idea of the numerical uncertainties in our calculation of the various asymmetries). The result for the n = 1 asymmetry (black line) gives a non-vanishing value (it corresponds to the integral of the differential cross section component [dσ] 3 in eq. (3.1)). A non-vanishing value is obtained also for n = 2 (red line), and it corresponds to the computation of the cross section component [dσ] 2 in eq. (3.1). Note that the n = 2 result reported in figure 1-left is rescaled by a factor of 0.1. Therefore, by inspection of figure 1-left we can see that the cos(2ϕ) asymmetry is approximately a factor of 5 larger than the cos ϕ asymmetry.
In figure 1-left we also report the result for the n = 0 harmonic (dashed line) of 'Z + 1 parton' production. The n = 0 harmonic corresponds to the total (azimuthallyintegrated) cross section, and it receives contributions from both real and virtual emission subprocesses. Real and virtual terms are separately divergent and their divergences cancel in the total contribution. In figure 1-left we report the result of the n = 0 harmonic computed exactly as specified in eq. (2.11), namely, by applying a non-vanishing lower limit q cut = r cut M on q T . Therefore, our computation selects only the real-emission term due to the 'Z + 1 parton' subprocesses. As is well known (see also section 4.1), this realemission term diverges in the limit q cut → 0 and the divergent behaviour is proportional to ln 2 r cut : in figure 1-left we clearly see the increasing behaviour of the n = 0 harmonic as r cut → 0. In the actual computation of the total cross section, the n = 0 result in figure 1-left has to be combined with the NLO real-emission term at still smaller values of r cut and with the LO and NLO virtual terms, thus leading to a total finite result. For the sake of completeness, we report the value (including the numerical error of the Monte Carlo integration) of the NLO total cross section σ N LO DY that we obtain by using the same parameters as used in the results of figure 1-left: it is σ N LO DY = 1100 ± 1 pb. We note that σ N LO DY is roughly 100 times larger than the value of the cos(2ϕ) asymmetry in figure 1-left. The DY process is quite 'special' with respect to azimuthal correlations: the azimuthal correlations are finite order-by-order in QCD perturbation theory and their general QCD structure involves only 4 azimuthal harmonics (as shown in eq. (3.1)). For most of the other JHEP06(2017)017 hard-scattering processes (with few possible exceptions such as H → γγ, ZZ production, as discussed in section 2) azimuthal correlations behave differently: usually they have an azimuthal dependence that involves an infinite set of harmonics (all values of n) and in many cases f.o. QCD computations lead to divergences.
We consider the production of heavy-quark QQ pairs and we treat the heavy quark and antiquark as on-shell particles. Our discussion equally applies to any heavy quark, but we specifically consider top quarks since the on-shell treatment is more suitable in this case. In figure 1-right we present results for tt production that are obtained analogously to the DY results presented in figure 1-left.
Our QCD computation of tt production is performed by using the numerical Monte Carlo program of ref. [26], which includes QCD radiative corrections at NLO (it uses the NLO scattering amplitudes of the MCFM program [27]) and part of the NNLO contributions. We set µ F = µ R = m t and we use the value m t = 173.3 GeV for the mass m t of the top quark. We consider the azimuthal angle ϕ of the top quark in the CS frame, and in figure 1-right we present the numerical results of various cos(nϕ) asymmetries (see eq. (2.11)) after their integration over the invariant mass M of the tt pair. As in the case of the results in figure 1-left, the tt results are integrated over the polar angle of the top quark and over the rapidity y and transverse momentum q T (q T > q cut ) of the tt pair. The azimuthal asymmetries are evaluated (as a function of r cut = q cut /M ) at their non-trivial lowest order (i.e., O(α 3 S )) in f.o. perturbation theory and, therefore, we compute all the NLO tree-level partonic processes whose final state is 'tt + 1 parton'.
Before commenting on the results in figure 1-right, we recall some of the results in ref. [12]. In ref. [12] azimuthal correlations for tt production are computed in analytic form in the small-q T limit. At O(α 3 S ) the result of ref. [12] shows that the q T -dependent azimuthal-correlation cross section dσ corr /dq 2 T dϕ behaves proportionally to 1/q 2 T in the limit q T → 0 and, hence, it is not integrable over q T down to the region where q T = 0. The 1/q 2 T behaviour is proportional to a non-polynomial function of cos 2 ϕ that, therefore, leads to divergent cos(nϕ) harmonics for even values of n (n = 2, 4, 6, . . . ). The q T spectra of harmonics with odd values of n have instead a less singular q T behaviour at O(α 3 S ) and they are integrable over q T in the limit q T → 0.
The numerical results in figure 1-right are consistent with the convergent or divergent behaviour predicted in ref. [12]. The harmonic with n = 1 (black line) is not vanishing and basically independent of r cut for very small values of r cut . The sign of the n = 1 harmonic is negative (the n = 1 harmonic oft would be positive, analogously to the corresponding harmonic of the electron in the DY case of figure 1-left), and its absolute size (note that it is rescaled by a factor of 0.1 in figure 1-right) is roughly a factor of two smaller than the size of the n = 1 harmonic for Z(e + e − ) production (figure 1-left). The harmonics with n = 2 (red line), n = 4 (blue line) and n = 6 (magenta line) in figure 1-right have instead an increasing (in absolute value) size for small and decreasing values of r cut : this behaviour is consistent with a ln r cut dependence, as expected from the analytical results in ref. [12]. The results for n = 2, 4 and 6 in figure 1-right have no straightforward quantitative implications for physical azimuthal asymmetries since they refer to small values of r cut (and they eventually diverge in the limit r cut → 0). Nonetheless we observe that the absolute magnitude of the JHEP06(2017)017 n-even cos(nϕ) asymmetries decreases as n increases. As in the case of figure 1-left for the DY process, in figure 1-right we also present the O(α 3 S ) result of the n = 0 harmonic (dashed line) for the real-emission process 'tt + 1 parton'. Analogously to the DY process, the 'tt + 1 parton' contribution to the n = 0 harmonic diverges in the limit r cut → 0, and its dominant behaviour at small r cut is proportional to ln 2 r cut + O(ln r cut ). At small values of r cut the shape of the r cut -dependence of the n = 0 result is thus steeper than that of the results for the harmonics with n = 2, 4 and 6. After combining real and virtual contributions, the n = 0 harmonic gives the total cross section. The value (including the numerical error of the Monte Carlo integration) of the NLO total cross section that we find is σ N LO tt = 226.2 ± 0.1 pb. We note that σ N LO tt is roughly 200 times larger than the absolute value of the cos ϕ asymmetry (n = 1) in figure 1-right.
As we have anticipated in section 2 with a brief sentence, we expect (and predict) the appearance of f.o. perturbative divergences in the computation of QED radiative corrections to azimuthal correlations for the DY process. To explicitly explain this point we consider some analogies of the DY and tt production processes. We have illustrated the O(α 3 S ) divergences of the n-even cos(nϕ) asymmetries for tt production. Part of these divergences arise from the computation of the partonic subprocess qq → tt + g and, specifically, they arise from the kinematical region where the radiated final-state gluon g is soft and at wide angle with respect to the direction of the initial-state q andq. From this kinematical region the tt spectrum dσ/dq 2 T dϕ receives an azimuthal-correlation contribution that behaves as 1/q 2 T and that depends on the QCD colour charges of the final-state t andt. Specifically, this contribution is proportional to the Fourier transformation of the function D (1) in eq. (36) of ref. [12]. The first-order QED radiative corrections to the DY process involve analogous partonic processes, such as qq → Z(e + e − )+γ, with a soft and wide-angle photon γ that is radiated in the final state. These photon radiative corrections produce f.o. QED divergences to azimuthal asymmetries for the DY process. The DY analogue of the tt function D (1) in eq. (36) of ref. [12] is simply obtained by replacing the QCD colour charges of t andt with the QED electric charges of the DY final-state leptons (by simple inspection of eq. (36) in ref. [12], such replacement leads to a non-vanishing result).
More generally, we remark that the expression in eq. (3.1) for the DY process is valid only at the Born level with respect to the EW interactions (although the expression is valid at arbitrary orders in QCD perturbation theory). Including QED radiative corrections, we conclude that the azimuthal-correlation cross section dσ DY /dq 2 T dϕ for the DY process receives contributions from cos(nϕ) and sin(nϕ) harmonics with arbitrary values of n (not only n = 1, 2 as in eq. (3.1)). Moreover, the lowest-order perturbative QED computation of the cos(nϕ) contributions with n-even already leads to a singular (not integrable) behaviour in the limit q T → 0 and to ensuing divergent azimuthal asymmetries upon integration over q T . The divergences are switched on by QED radiative corrections and can receive additional contributions from powers of α S in the context of mixed QED×QCD radiative corrections. Obviously, similar divergences arise also in a pure QED context such as, for instance, in the QED computation of the process e + e − → µ + µ − + X. Eventually all these divergences originate from the QED analogue of the condition (2.10): the QED divergences arise if 'at least one the final-state particles with momenta p 3 and p 4 carries non-vanishing electric charge'.

JHEP06(2017)017
We note that lepton pairs with high invariant mass can also be produced throughout the partonic subprocess γγ → + − , in which the initial-state colliding photons arise from the photon PDF of the colliding hadrons. Strictly speaking, lepton pairs that are produced in this way are not physically distinguishable 4 from those that are produced by the DY mechanism of quark-antiquark annihilation. Since the ratio between the photon PDF and the quark (or antiquark) PDF is formally of O(α/α S ) (α is the fine structure constant), the subprocess γγ → + − can be regarded as a QED radiative correction in the context of the DY process. We remark that radiative corrections to γγ → + − also produce divergences in the computation of the lepton azimuthal correlations. These divergences (see section 5.1) originate from the abelian (photonic) analogue of the condition (2.9): the divergences arise if 'at least one of the initial-state colliding partons c 1 and c 2 is a photon'.
To cure the f.o. perturbative divergences of azimuthal correlations one may advocate non-perturbative strong-interactions dynamics and related non-perturbative QCD effects, which can be sizeable in the small-q T region. However, in the case of q T -integrated azimuthal correlations (see eq. (2.7)) the non-perturbative QCD dynamics should cancel divergent terms proportional to some powers of α S (M ), and this would imply that nonperturbative QCD effects scale logarithmically with M (i.e., these effects would not be suppressed by some power of Λ QCD /M in the hard-scattering regime M Λ QCD ), thus spoiling not only the finiteness but also the infrared safety of the azimuthal correlations. Moreover such non-perturbative cure of the divergent behaviour cannot be effective in the case of QED radiative corrections and, in particular, it cannot be conceptually at work in a pure QED context since QED is well-behaved in the infrared (small-q T ) region. In the following sections we show that the problem of f.o. divergences in azimuthal correlations has a satisfactory solution entirely within the context of perturbation theory. Namely, the resummation of perturbative corrections to all orders leads to (q T integrated) azimuthal asymmetries that are finite and computable. Such solution is effective in both QCD and QED contexts, although in the QCD case this does not imply that the non-perturbative effects have a negligible quantitative role in the small-q T region.

The small-q T region
The f.o. divergences of azimuthal correlations arise from the small-q T region and, eventually, from the q T behaviour at the phase space point where q T = 0. From a physical viewpoint, however, a single (isolated) phase space point is harmless and, in particular, non-vanishing azimuthal correlations (which are defined with respect to the direction of q T ) cannot be measured if q T ∼ 0. These physical requirements imply that the q T dependence of azimuthal-correlation observables has to be sufficiently smooth in the small-q T region and, in particular, in the limit q T → 0.
We specify more formally these smoothness requirements. Since physical measurements cannot be performed at a single phase space point, we split the q T integration region in the JHEP06(2017)017 two intervals where 0 ≤ q T ≤ q cut (the lower-q T bin) and q T ≥ q cut (the higher-q T region). We then consider a generic azimuthal-correlation observable (e.g., multidifferential cross sections as in eqs. (2.3) and (2.4), or n-th harmonics as in eq. (2.11)) and we denote by σ l (q cut ; ϕ) and σ h (q cut ; ϕ) the 'cross sections' that are obtained by q T integration of the observable over the lower-q T bin and the higher-q T region, respectively. The smoothness requirements are where σ tot (ϕ) denotes the 'total' cross section (i.e., the result obtained by integrating the observable over the entire region of q T ). These requirements specify azimuthal correlations that are physically well behaved. In particular, eq. (4.1) implies that non-vanishing azimuthal correlations cannot be physically observed if q T ∼ 0, and eq. (4.2) implies that the total azimuthal correlation is physically observable. Note, however, that eqs. (4.1) and (4.2) do not specify how dσ/dq 2 T dϕ exactly behaves in the limit q T → 0. Nonetheless, eqs. (4.1) and (4.2) imply that the sole phase space point at q T = 0 (where azimuthal-correlation angles are not defined) has no relevant physical role (the lower-q T bin has an analogous harmless role if q cut is sufficiently small).
The f.o. divergences of azimuthal-correlation observables are due to the fact that dσ/dq 2 T dϕ does not have a ('sufficiently') smooth dependence on q T at small values of q T if this quantity is computed order-by-order in perturbation theory. In particular, σ h (q cut ; ϕ) diverges in the limit q cut → 0 and, consequently, σ tot (ϕ) is not computable, whereas σ l (q cut ; ϕ) is definitely not computable (divergent) even if q cut has a finite value.
As we have observed in section 3, also the customary azimuthally-integrated (or azimuthally-averaged) cross section (i.e., the harmonic with n = 0 in eq. (2.11) or figure 1) does not fulfil the requirements in eqs. (4.1) and (4.2) (both σ l (q cut ; ϕ) av. and σ h (q cut ; ϕ) av. can separately be divergent in the limit q cut → 0). Therefore, the f.o. calculation of azimuthally-integrated cross section can have difficulties (as is well known) in describing the detailed q T shape in the small-q T region. In contrast to azimuthal-correlations observables, however, both σ l (q cut ; ϕ) av. and σ h (q cut ; ϕ) av. are computable for finite values of q cut . In particular, the total cross section σ tot (q cut ; ϕ) av. = σ l (q cut ; ϕ) av. + σ h (q cut ; ϕ) av. is always finite and computable (for arbitrary non-vanishing values of q cut ) within f.o. perturbation theory.
In the following sections we discuss the small-q T behaviour of cross sections in f.o. perturbation theory and after all-order resummation.

Perturbative expansion with azimuthal-correlation terms
Among all the hadroproduction processes of the type in eq. (2.1), we consider those that at the LO in perturbative QCD are produced by the following partonic subprocesses: JHEP06 (2017)017 where c 1 and c 2 (c i = q,q, g) are the initial-state massless colliding partons of the two hadrons h 1 and h 2 . In the case in which one or both particles (with momenta p 3 , p 4 ) in F is a jet, the notation in eq. (4.3) means that the jet is replaced by a corresponding QCD parton. Note that the LO process in eq. (4.3) is an 'elastic' production process, in the sense that F is not accompanied by any additional final-state radiation. 5 This specification is not trivial, since it excludes some processes from our ensuing considerations. Some processes are excluded because of quantum number conservation. For instance, if F includes two top quarks (not a tt pair) no LO process as in eq. (4.3) is permitted by flavour conservation. Some other processes are excluded because of their customary perturbative QCD treatment. For instance, if F contains a hadron, its QCD treatment requires the introduction of a corresponding fragmentation function and, consequently, F is necessarily produced with accompanying fragmentation products in the final state. Therefore, in the following we exclude the cases in which F includes one or two hadrons. 6 All the processes that are explicitly listed in eq. (2.6) (including the DY process) have LO partonic processes of the type in eq. (4.3).
The small-q T behaviour of azimuthal correlations for the processes that we have just specified is partly related to the behaviour of q T differential cross sections that are integrated over the entire azimuthal-angle region (or, equivalently, that are azimuthally averaged). The behaviour in the azimuthally-integrated (averaged) case is well known (starting from some of the 'classical' studies for the DY process [28][29][30][31][32][33][34][35]), and in our presentation we contrast it with the cases with divergent azimuthal correlations. To remark the differences between the two cases in general terms, we use a shorthand notation, which can be applied to final-state systems F with two or more particles (a more refined kinematical notation can be found in refs. [11,12]). A generic multidifferential cross section (e.g., eq. (2.3), eq. (2.4) or related observables) with dependence on q T and on the azimuthal-correlation angle is simply denoted by dσ/dM 2 d 2 q T , where q T is the transverse momentum of F . Additional kinematical variables that are possibly not integrated (such as rapidities and polar angles) are not explicitly denoted. The dependence on a generic (as discussed in section 2) azimuthal-correlation angle with respect to q T is denoted through the dependence on the directionq T = q T /q T of the two-dimensional vector q T . In practice such dependence will occur through functions of scalar quantities such as, for instance,q T ·p 3T . In particular, independent ofq T means absence of azimuthal correlations, and the azimuthally-integrated (averaged) cross section corresponds to the azimuthal integration (average) with respect to the directionq T of q T .
Owing to transverse-momentum conservation in the process of eq. (4.3), at the LO level F is produced with vanishing q T . Its corresponding LO q T cross section is proportional to As usual, the final-state collinear remnants of the colliding hadrons are not denoted in the partonic process. 6 Specifically, our ensuing discussion does not apply to azimuthal correlations of systems F that contain hadrons with momenta p3 or p4 (whereas it applies to infrared and collinear safe jets, which physically contain hadrons).

JHEP06(2017)017
and the proportionality factor is simply the LO total (q T integrated) cross section. In the presence of the LO sharp q T behaviour of eq. (4.4), NLO QCD radiative corrections are dynamically enhanced in the small-q T region. The dynamical enhancement is due to QCD radiation of low transverse-momentum partons (soft gluons or QCD partons that are collinear to the initial-state colliding partons), and it has the following general form: where the dots on the right-hand side denote contributions that are less enhanced ('nonsingular') in the limit q T → 0. The 'coefficients' a i (i = 0, 1, 2) and a corr in eq. (4.5) depend on the process and they are independent of q T (we mean they are independent of the magnitude q T of the vector q T ). The important point (see below) is that a corr does depend on the directionq T of q T , whereas a i (i = 0, 1, 2) do not depend on it. All these coefficients (a 1 , a 2 .a 3 , a corr ) can depend on the other kinematical variables (e.g., rapidities and polar angles). The QCD running coupling α S in eq. (4.5) and in all the subsequent formulae is evaluated at a scale of the order of M (e.g., we can simply assume α S = α S (M )).
The terms that we have explicitly written in the right-hand side of eq. (4.5) scale as 1/q 2 T (modulo logs) in the limit q T → 0 (i.e., they scale as 1/λ 2 under the replacement q T → λq T ). The non-singular terms (which are simply denoted by the dots) represent subdominant ('power-correction') contributions in the limit q T → 0. These are, for instance, terms of the type M/q T or, more generally, terms that are relatively suppressed by some powers (modulo logs) of q T /M . Independently of their specific form and of their azimuthal dependence, these non-singular terms have an integrable and smooth behaviour in the limit q T → 0. Because of these features, the non-singular terms and the ensuing azimuthal-correlation effects are well behaved in the small-q T region. Note, however, that the non-singular terms produce azimuthal effects whose actual size (and q T behaviour) depends on the specific definition of the azimuthal-correlation angle (see eq. (2.8) and related comments in section 2).
The expression in eq. (4.5) is the master formula that we use for our discussion of the small-q T behaviour of azimuthal correlations at NLO and higher orders. Considering the 'singular' terms (the NLO terms that scale as 1/q 2 T , modulo logs, in the limit q T → 0), the azimuthal dependence is entirely embodied in a corr (q T ). More precisely, the azimuthalcorrelation dependence has been separated (as in eq. (2.7)) and embodied in a corr (q T ) that, therefore, gives a vanishing result after azimuthal integration. Using the notation in eq. (2.7) we have a corr (q T ) av. = 0 . This azimuthal-correlation term is absent for the DY process (i.e., a corr (q T ) = 0 in this case), and in all the other processes it has no effect by considering q T dependent but azimuthally-integrated observables. The presence of such term produces a q T behaviour that definitely differs from the behaviour studied in the 'classical' literature [28][29][30][31][32][33][34][35] on the small-q T region and on QCD transverse-momentum resummation.

JHEP06(2017)017
The other ('classical') terms in eq. (4.5) are proportional to the coefficients a i (i = 0, 1, 2), and the symbol [ ] + denotes the ('singular') 'plus'-distribution, which is defined by its action onto any smooth function f (q T ) of q T . The definition is where g(q T ) = 1 q 2 T ln m (M 2 /q 2 T ) (m = 0, 1) in eq. (4.5), and µ 0 is a scale of the order of M (at the formal level, varying the value of µ 0 changes the plus-distributions and the coefficient a 0 , so that the right-hand side of eq. (4.5) is unchanged). The plus-distribution in eq. (4.7) is equivalently defined through a limit procedure as follows: (4.8) Note that the point at q T = 0 is at the border of the phase space, but, at the strictly formal level, it is inside the phase space (at variance with the case of azimuthal correlations, in which it is formally outside the phase space). This is essential to make sense of the LO q T differential cross section in eq. (4.4) and of the corresponding LO total cross section. This is also essential (see eqs. (4.7) and (4.8)) to transform the non-integrable q T behaviour of 1 into an integrable plus-distribution over the small-q T region. The plusdistribution [ g(q T ) ] + involves two terms: a term that is simply proportional to the function g(q T ) and a contact term (the term proportional to f (0) in eq. (4.7) or, equivalently, the term proportional to δ (2) (q T ) in eq. (4.8)). At the conceptual level, these two terms arise from a combination of real and virtual radiative corrections, which are separately infrared divergent. Real (soft and collinear) emission corrections to the LO process in eq. (4.3) produce the non-integrable terms 1 q 2 T ln m (M 2 /q 2 T ), while the corresponding virtual radiative corrections lead to the contact terms that eventually regularize the divergence at q T = 0.
The azimuthal-correlation term in eq. (4.5) is instead proportional to the divergent (not integrable) function 1/q 2 T , rather than to a corresponding plus-distribution. At the formal level one may be tempted to transform the azimuthal-correlation term into a plusdistribution through the replacement a corr (q T )/q 2 T → [ a corr (q T )/q 2 T ] + , but such a replacement would not be effective since namely, the plus-prescription is not able to regularize the non-integrable behaviour of azimuthal-correlation terms. The equality in eq. (4.9) is a consequence of the fact that the contact term (see eq. (4.8)) in the corresponding plus-prescription identically vanishes: indeed we have because the azimuthal integration of a corr (k T ) gives a vanishing result 7 (see eq. (4.6)). We note that a result analogous to eq. (4.9), namely [g corr (q T )] + = g corr (q T ) , (4.11) is valid for any azimuthal-correlation function g corr (q T ), i.e., any function such that g corr (q T ) av. = 0. Therefore, in general we can conclude that azimuthal-correlation terms that are not integrable in the limit q T → 0 cannot be regularized by the plusprescription. This general conclusion and, in particular, the equality in eq. (4.9) have a simple conceptual interpretation in the context of f.o. perturbation theory. Virtual radiative corrections to the process in eq. (4.3) have q T = 0 and, therefore, they cannot contribute to azimuthal-correlation terms. In particular, they cannot provide the real-emission azimuthal-correlation term with a (non-vanishing) contact term that transforms the nonintegrable real-emission contribution into an integrable plus-distribution. In summary, the azimuthal correlation of the q T cross section is a phenomenon that is necessarily produced by real emission (since azimuthal dependence requires q T = 0) and the f.o. divergences of the azimuthal asymmetries are the consequence of a complete mismatch with virtual emission, which is completely absent (i.e., it does not contribute) at the corresponding perturbative order. In the case of azimuthal-independent (plus-distribution) terms, instead, both real and virtual effects contribute at the NLO, although their relative contribution is highly unbalanced in the small-q T region. The r cut (q cut ) dependence of the results that are presented in figure 1 is fully consistent with the small-q T behaviour in eq. (4.5). The harmonics with n = 0 receive contributions only from the term that is proportional to a corr (q T ) in eq. (4.5), while a corr (q T ) does not contribute to the n = 0 harmonics (because of eq. (4.6)). We simply note that in the DY process we have a corr (q T ) = 0, whereas in the case of tt production a corr (q T ) is not vanishing (in particular, the n = 1 harmonic of a corr (q T ) vanishes, while a corr (q T ) has nonvanishing harmonics with n = 2, 4, 6). We also note that the results for the harmonics with n = 0 (in both the DY and tt production processes) follow from the q T shape of the plusdistribution terms in eq. (4.5): the q T integration of the plus-distributions over the region q T > q cut = r cut M produces the r cut dependence in figure 1, and this r cut dependence is exactly cancelled by that produced from the integration of the plus-distributions over the region 0 ≤ q T ≤ q cut .
We include a comment on an additional point related to the f.o. perturbative expansion of generic hard-scattering processes. In our discussion of eq. (4.5), we have considered processes that at the LO are initiated by the partonic subprocesses in eq. (4.3). We note that elastic production processes of the type in eq. (4.3) can also appear at some higher perturbative orders. For instance, as already mentioned in section 2, this is the case if F = {γγ}, {ZZ}, {W + W − }. These diboson production processes are initiated at LO by qq → F , and they also receive contributions from gg → F (the gluons are 7 The integral in eqs. (4.8) and (4.10) may be interpreted as the integral over the small-kT region of the loop momentum k T of the virtual NLO correction. After integration over the azimuthal angle of k T , only the azimuthal average gives a non-vanishing contribution to eq. (4.8), while the azimuthal-correlation effect in eq. (4.10) vanishes.

JHEP06(2017)017
coupled to γγ, ZZ, W + W − through a quark loop 8 ) starting from the NNLO. In these cases, the next-order radiative corrections to gg → F produce a small-q T behaviour that has exactly the same structure as that in the round-bracket term of eq. (4.5). Our following discussion throughout the paper equally applies to all the processes that can occur through subprocesses of the type in eq. (4.3) at some perturbative order (different types of colliding partons c 1 and c 2 can be involved at different perturbative orders).
Small-q T singular terms of the type in eq. (4.5) are present in the computation of dσ/dM 2 d 2 q T at each higher perturbative orders. The 1/q 2 T behaviour of both azimuthally-independent and azimuthally-correlated terms is enhanced by logarithmic factors, ln k (M 2 /q 2 T ), produced by multiple radiation of soft and collinear partons. There are at most two additional powers of ln(M 2 /q 2 T ) for each additional power of α S , so that the dominant singular terms have a double-logarithmic structure. The N k LO radiative corrections to dσ/dM 2 d 2 q T include singular terms that are proportional to α k , and this q T dependence is regularized by a plus-prescription only in the case of azimuthally-independent contributions. Owing to eq. (4.11), the plus-prescription is not effective for azimuthal-correlation terms.
The order-by-order singular behaviour that we have just discussed poses no problems for a basic quantity such as the total (q T and azimuthally integrated) cross section: the singular azimuthal-correlation terms cancel because of the azimuthal integration, and the plus-distribution terms are integrable over q T . As is well known, the plus-distribution terms (which remain after azimuthal integration) are not 'harmless' for differential cross sections. They lead to a sharp unphysical behaviour of the small-q T differential cross section at the NLO, and to large radiative-correction effects at any subsequent perturbative order, with an unstable f.o. expansion and poor predictivity for the detailed shape of the q T cross section in the small-q T region. The smooth physical behaviour of the q T cross section and the predictivity of perturbative QCD can be recovered by the all-order resummation (transverse-momentum resummation) of the logarithmically-enhanced plusdistribution terms.
In the case of azimuthally-sensitive observables, the disease of f.o. perturbation theory is more serious. The singular azimuthal-correlation terms are not integrable and, as discussed in section 2, basic quantities such as the total (q T integrated) cross section at fixed azimuthal angle (see eq. (2.6)) and the total (q T integrated) azimuthal asymmetries (see eq. (2.11)) can be divergent if they are evaluated in f.o. perturbation theory. In the presence of divergences, also the f.o. perturbative result for dσ corr /dM 2 d 2 q T at fixed (and finite) values of q T has an unphysical shape in the small-q T region, since the f.o. result follows a non-integrable (unphysical) behaviour.
In the following we briefly recall known results on transverse-momentum resummation for 'azimuthally-insensitive' observables. More precisely, we simply sketch some main points of transverse-momentum resummation that are relevant for our subsequent discussion on azimuthal correlations and asymmetries.

Transverse-momentum resummation with no divergent azimuthal correlations
In this subsection we limit ourselves to considering only singular contributions to dσ/dM 2 d 2 q T that are independent of the azimuthal-correlation angles. Order-by-order in perturbation theory, these contributions can be additively separated from non-singular terms (see, e.g., the 'dots' in the right-hand side of eq. (4.5)), and they can also be additively separated from singular azimuthal-correlation terms (see eqs. (2.7), (4.5), (4.11) and accompanying comments). Therefore, we are dealing with all the singular plus-distribution terms that appear in eq. (4.5) and in corresponding higher-order contributions. As a consequence of the additive separability, azimuthally-insensitive, azimuthally-integrated and azimuthally-averaged contributions have the same equivalent meaning in the context of the discussion in this subsection. In the case of processes whose azimuthal correlations are not singular (such as the DY process in the context of QCD radiative corrections), the singular terms that we are considering are the entire singular contributions to those processes.
Considering hard-scattering observables that are insensitive to divergent azimuthal correlations, the all-order QCD treatment of the singular contributions to the q T differential cross section in the small-q T region is conceptually well known [30,31,35]). In the cases in which the produced high-mass system F is formed by particles that carry no QCD colour charge (colourless system F ), for instance, in the cases of the DY process and of SM Higgs boson production, the treatment is fully developed at arbitrary logarithmic accuracy by using various methods and formalisms, such as direct QCD resummation (see ref. [36] and corresponding references therein), Soft Collinear Effective Theory (SCET) methods (see, e.g., refs. [37,38]), and transverse-momentum dependent (TMD) factorization (see, e.g., refs. [39][40][41][42]). In the case of processes whose final-state system F includes colour-charged particles (colourful system F ), the theoretical treatment is formally less advanced and only few specific cases have been treated beyond the leading-logarithmic (LL) level. For the specific process of heavy-quark pair production (F = {QQ}), transverse-momentum resummation has been explicitly developed [12,43,44] up to next-to-next-to-leading logarithmic (NNLL) accuracy. The process of dijet production (F = {jj}) has been explicitly studied [45,46] up to next-to-leading logarithmic (NLL) accuracy within the approximation of small values of the cone size R of the jets.
To the purposes of our subsequent discussion on azimuthal correlations, we limit ourselves to recalling some features of the transverse-momentum resummation formalism. Transverse-momentum resummation has to be carried out [30] in impact parameter (b) space to exactly implement the relevant constraint of transverse-momentum conservation for all-order multiparton radiation in the inclusive final state. The terms with plusdistributions (as in eq. (4.5)) at small values of q T become powers of logarithms, ln(bM ), in b space at large values of b (bM 1). These logarithmically-enhanced contributions can be resummed in b space. Then the q T cross section dσ/dM 2 d 2 q T is obtained by inverse where J 0 (x) is the 0th-order Bessel function, and the superscript 'res' denotes the resummed contribution to the q T differential cross section. Since we are dealing with the resummation of the azimuthally-insensitive (independent) plus-distribution terms, we have introduced the subscript 'az. in.', and we simply point out that these terms equally contribute to the azimuthally-averaged ('az. av.') q T cross section. The integrand Σ az.av. involves process-dependent and process-independent factors (see, e.g., refs. [11,12,35]). One of them is the Sudakov form factor S c (M, b) (c = q, g) of the colliding partons c 1 and c 2 of the partonic subprocess in eq. (4.3). The Sudakov form factor is universal (it is independent of F ) and it includes the entire effect of the resummation of the dominant double-logarithmic (DL) contributions (from radiation of partons that are both soft and collinear to the initial-state colliding partons) in b space.
Within the DL approximation, the Sudakov form factor has the exponential form 9 e −α S ln 2 (bM ) [30] that produces a very strong damping of the large-b region (bM 1) in the integrand on the right-hand side of eq. (4.12). The suppression of the large-b region is stronger than that produced by any power of 1/(bM ). The damping effect of the Sudakov form factor eventually leads to resummed perturbative predictions for the q T cross section that are physically well behaved (with a smooth q T dependence) in the small-q T region: the shape of the q T cross section is the result of a (computable) smearing of the LO δ-function behaviour in eq. (4.4). In particular, the qualitative behaviour of dσ (res) az.av. /dM 2 dq 2 T at very low values of q T can be examined by performing the limit q T → 0 of eq. (4.12). Using 10 J 0 (bq T ) = 1 + O(bq T ) in eq. (4.12), we simply have [30] dσ (res) az.av.

JHEP06(2017)017
also depend on non-perturbative 11 effects (see, e.g., ref. [47]) that affect the region of very large values of b (b ∼ > 1/Λ QCD ). We note that experimental results on q T cross sections are usually presented in terms of the differential cross section dσ/dq T , rather than dσ/dq 2 T . The two cross sections are directly related by just a factor of 2q T from kinematics (dσ/dq T = 2q T dσ/dq 2 T ). The azimuthally-averaged (or azimuthally-integrated) cross section dσ az.av. /dq T has a peak in the small-q T region and it vanishes linearly in q T as q T → 0: this behaviour is the consequence of the combined effect of the dynamical small-q T behavior in eq. (4.13) and of the kinematical suppression factor 2q T as q T → 0.
We have previously noticed (in sections 2 and 4.1) that in many processes the finalstate system F can be elastically produced (see eq. (4.3)) by subprocesses with different initial-state partonic channels. Our discussion about resummation equally applies to all these processes: the corresponding q T cross section has several components (one component for each contributing partonic channel) (see, e.g., refs. [11,12]) and each component has basically the same structure as in eq. (4.12).

Origin of divergent azimuthal correlations
In the context of QCD perturbation theory, singular azimuthal correlations in the small-q T region were observed in the theoretical studies of refs. [10][11][12].
Reference [10] deals with diphoton production (F = {γγ}) in hadron-hadron collisions. Considering the next-order radiative corrections (which are part of the complete N 3 LO corrections for diphoton production) to the gluon fusion subprocess gg → γγ, the authors of ref. [10] find that the cos(2ϕ) modulation of the azimuthal-dependent q T cross section dσ/dq 2 T dϕ behaves as 1/q 2 T at small values of q T . The production of a generic system F of two or more colourless particles (for instance, F = {γγ}, {ZZ}, {W + W − }) is considered in ref. [11]. Since the particles in F have no QCD colour charge, the elastic-production subprocesses (of the type in eq. (4.3)) that are permitted by colour conservation are qq → F and gg → F . The study of the small-q T region performed in ref. [11] shows that QCD radiative corrections to the gluon fusion channel, gg → F , lead to azimuthal correlations with a singular q T behaviour (the singular behaviour is instead absent for QCD corrections to the quark-antiquark annihilation channel, qq → F ). The singular behaviour is proportional to 1/q 2 T at the lowest order (i.e., the corresponding a corr (q T ) of eq. (4.5) is not vanishing) and it is enhanced by doublelogarithmic factors, ln 2 (M 2 /q 2 T ), at higher orders [11]. Not all the azimuthal harmonics have a singular behaviour in the limit q T → 0 (see, in particular, eqs. (76)-(86) and accompanying comments in ref. [11]): the singular harmonics are cos(2ϕ) (and sin(2ϕ)) starting from the lowest order, and cos(4ϕ) (and sin(4ϕ)) at higher orders (the sin(2ϕ) and sin(4ϕ) harmonics typically receive contribution from parity-violating effects).

JHEP06(2017)017
The study of ref. [12] considers heavy-quark pair production (F = {QQ}) in the small-q T region. The allowed elastic-production subprocesses (of the type in eq. (4.3)) are qq → QQ and gg → QQ, as in the case of the colourless systems F considered in ref. [11]. The important difference with respect to ref. [11] is that the produced heavy quarks (Q andQ) carry QCD colour charge, so that they are sources of additional QCD radiation and produce additional dynamical effects (both final-state effects and initial/finalstate quantum interference effects). Singular azimuthal correlations, analogous to those in ref. [11], are produced by QCD radiative corrections to the gluon fusion channel gg → QQ. Additional azimuthal correlations, with the same small-q T singular behaviour, are produced by the non-vanishing colour charges of Q andQ, and they are generated through QCD radiative corrections to both production channels qq → QQ and gg → QQ. The explicit lowest-order (which corresponds to a corr (q T ) in eq. (4.5)) and higher-order results of ref. [12] show that the small-q T singular behaviour of these additional azimuthal correlations affects the cos(nϕ) harmonics with arbitrary even values of n (n = 2, 4, 6, . . . ).
The analyses and results of refs. [11,12] involve some process-dependent features. However, the sources of singular azimuthal correlations that are identified therein have a process-independent (universal) dynamical origin. This allows us to draw some general conclusions for the entire class of processes discussed in sections 2 and 4.1 (some of these processes are mentioned in eq. (2.6)). There are certainly two sources of small-q T singular azimuthal correlations and related divergent azimuthal asymmetries for these processes. These two sources are i) initial-state collinear radiation from gluonic colliding states that produce the particles of the high-mass system F ; ii) soft wide-angle radiation (both final-state contributions and initial/final-state interference contributions) generated by colour-charged particles of the high-mass system F .
The process-independent origin of singular azimuthal correlations that we have just specified leads to the conditions (2.9) and (2.10) that we have anticipated in section 2. Specifically, the collinear origin at the point i) leads to the condition (2.9), and the soft origin at the point ii) leads to the condition (2.10). Some main consequences of the conditions (2.9) and (2.10) are already discussed in the subsequent part of section 2, and they are not repeated here. In the following we only present some additional comments. As for the initial-state collinear radiation at the point i), we refer the reader to eqs. (33)-(43) (and accompanying comments) in ref. [11] for a more detailed discussion of its process-independent dynamical origin. Collinear radiation from the initial-state colliding partons (c 1 and c 2 in eq. (4.3)) directly affects the kinematics (q T and rapidity) of the entire final-state system F , but it has no direct kinematical effects on the azimuthal dependence of the individual particles in F . The azimuthal dependence requires variations of the angular momentum of those particles. Therefore, collinear radiation can dynamically produce azimuthal correlations only through spin correlations induced by the collinear-emission process. Gluonic collinear splitting processes are indeed intrinsically spin JHEP06(2017)017 polarized. The produced linear polarization of the exchanged (t-channel) gluons leads to quantum interferences (between scattering amplitudes and their complex-conjugate amplitudes) that generate n = 2 (single helicity-flip contributions [10,11]) and n = 4 (double helicity-flip contributions [11]) azimuthal harmonics. The singular azimuthal correlations are due to the initial-state partonic splitting subprocess a → g + X (a = q,q, g is spin unpolarized and X denotes the accompanying collinear radiation) that dynamically generates a linearly-polarized colliding gluon. Singular azimuthal correlations are analogously produced by the (QED induced) collinear splitting subprocess a → γ + X, which generates a linearly-polarized colliding photon. In contrast, quark (antiquark) collinear splitting processes (a → q(q) + X) turn out to be spin unpolarized, as a consequence of helicity conservation in QCD radiation from a massless quark (antiquark). This difference between the collinear evolution of gluons (or photons) and quarks (antiquarks) is responsible for the specification of gluon (or photon) initial states at the point i) and in the condition (2.9). Obviously, QCD radiation from quarks can also produce spin correlations (and ensuing azimuthal asymmetries), as in the case of the DY process, but they are not singular in the collinear limit and, consequently, in the small-q T region.
The soft radiation at the point ii) refers to radiation of soft partons at wide angles with respect to the directions of the colliding partons and of the particles of the system F . If the system F contains only colourless particles, primary soft radiation only originates from the colliding partons c 1 and c 2 in eq. (4.3) (c 1 c 2 = qq or c 1 c 2 = gg, in this colourless case): this radiation is helicity conserving and it has a high degree of process independence [11]. If F contains colour-charged particles [12], they can radiate wide-angle gluons that kinematically produce azimuthal-correlation dependence. In these processes, we can have c 1 c 2 = qq, gg, qg and so forth, depending on the particles in F . Colour-charged particles in F and the colliding partons c 1 and c 2 simultaneously act as sources of primary QCD radiation, and the dynamical pattern of soft wide-angle radiation has a very high degree of process dependence. The process dependence originates by colour correlations with the underlying hard-scattering subprocesses. Moreover, within a single process, soft wide-angle radiation strongly depends on the kinematics of the particles in F because of initial/final-state interferences and related colour-coherence phenomena. All these features have a 'minimal' dependence on the actual partonic content of the colliding hadrons in the initial state.
We note that the collinear radiation and the soft radiation at the points i) and ii) have a different origin. However, starting from the second non-trivial perturbative order (F plus two or more additional partons in the final state), collinear radiation and soft radiation do not act as independent sources of singular azimuthal correlations. The two sources are indeed entangled by helicity and colour correlations, and they are analogously dynamically tangled with the helicity and colour structure of the underlying hard-scattering subprocesses (see, e.g., eqs. (22)- (25) and accompanying comments in ref. [12]).
Soft wide-angle radiation generally leads to divergent azimuthal harmonics with even values of n [12]. In principle, we have no reasons to justify the absence of divergent n-odd harmonics that can be produced by soft wide-angle radiation. In the case of QQ production we cannot exclude the appearance of divergent n-odd harmonics from subdominant JHEP06(2017)017 logarithmic contributions at high orders. We find it really difficult to justify the absence of divergent n-odd harmonics for processes with a more involved kinematical structure, such as those processes related to jet production (e.g., F = {V j} or F = {jj}). In the case of jets, the more involved kinematical structure arises, for instance, from jet selection cuts, and it certainly arises from the (infrared and collinear safe) procedure that is necessary to identify (define) hard-scattering jets.
To our knowledge, azimuthal correlations for jet production processes have not yet been studied throughout computations in QCD perturbation theory. In view of our comments on divergent harmonics with odd values of n, in the following we present an illustrative numerical investigation at the lowest perturbative order for the process of associated production of a vector boson V (V = Z, W ± , γ * ) and one jet j.
We consider the inclusive production process h 1 + h 2 → V j + X. The corresponding LO partonic processes are c 1 c 2 → V j (see eq. (4.3)), where c 1 c 2 = qq, qg,qg and j is a single final-state parton (with the corresponding flavour assignment, namely, g, q,q). The lowest-order contribution to azimuthal correlations is entirely due to all the NLO tree-level processes whose final state is 'V + 2 partons', where (after application of an infrared and collinear safe jet algorithm) one of the final-state partons is identified with the observed jet while the other parton contributes to the inclusive final-state radiation. As we have already discussed in section 2, both conditions (2.9) and (2.10) are fulfilled: therefore, at the lowest perturbative order we expect singular azimuthal correlations of both collinear (see point i) above) and soft (see point ii) above) origins. Note, however, that the collinear contribution can lead to singularities only for the cos(2ϕ) harmonic, since a singular cos(4ϕ) harmonic requires [11] a gg initial state at the LO and a final state with V j and at least two additional partons. The soft contribution is expected to be responsible for singular cos(nϕ) harmonics with all even values of n (n = 2, 4, 6, . . . ) [12] and, possibly, also with odd values of n.
We perform the lowest-order calculation of azimuthal harmonics by specifically considering the case of associated Zj production in pp collisions at the LHC ( √ s = 8 TeV). We use the numerical Monte Carlo program DYNNLO [25] with the same set-up and parameters as mentioned and used in section 3. We compute the azimuthal harmonics at fixed values of q T of the Zj system by integrating over all the other kinematical variables with the following selection cuts on the rapidity 12 η jet and transverse momentum p T jet of the jet j: |η jet | < 4.5 and p T jet > 20 GeV. We do not apply any lower limit on the invariant mass M of the Zj system since (as stated in section 3) we consider on-shell Z production and, consequently, M > M Z . The jets are defined by using the k ⊥ algorithm [48,49] or, equivalently, the anti-k ⊥ algorithm [50], since both algorithms give exactly the same results at this perturbative order. We use the value R = 0.4 of the jet radius R. We consider fully-inclusive V j production and, therefore, we do not select the observed jet j according to the hardness of its transverse momentum (the jet j of the Zj system is in turn one of the two partons of the final state 'V + 2 partons', independently of the value of its transverse JHEP06(2017)017 momentum). The conditions that we have just specified give an infrared and collinear safe definition of the cross section (and of related azimuthal correlations) for inclusive Zj production.
We define the azimuthal-correlation angle as ϕ = φ(p T jet )−φ(q T ), where φ(p T jet ) and φ(q T ) are the azimuthal angles of the transverse momenta p T jet and q T of the jet j and the system Zj in the centre-of-mass frame of the collision. As discussed and emphasized in section 2, this is one possible definition of the relevant azimuthal-correlation angle and the singular behaviour of azimuthal correlations is not affected by this choice. Note that we are not using the azimuthal-angle definition in the CS frame of the Zj system. Our choice of ϕ has a twofold purpose. Firstly, we take the chance of showing results for an alternative (with respect to the CS frame) definition of ϕ. Secondly, the definition of the CS rest frame of the Zj system requires the determination of the four-momentum of the jet (e.g., also the jet invariant mass matters) and such quantity is not always actually measured in experiments (anyhow, its determination is not necessary by using ϕ = φ(p T jet ) − φ(q T )).
We consider the q T spectrum dσ n /dq T of cos(nϕ) harmonics (see eq. (5.1)), and in figure 2 we present our lowest-order numerical results for Zj production. Note that in figure 2 we present the q T dependence of dσ n /dq T , whereas the azimuthal asymmetries in figure 1 are presented as functions of r cut = q cut /M and are obtained by integration of dσ n /dq T over q T with the lower limit q T > q cut . Note also that we are not showing results for dσ n /dq 2 T (see eqs. (4.5) and (5.1)), which is the q T differential cross section that we typically use throughout the text of this paper. In particular, since dσ n /dq T = 2q T dσ n /dq 2 T , the presence of singular harmonics (with n = 1, 2, . . . ) implies that dσ n /dq T ∝ 1/q T at the lowest perturbative order.
We comment on the results in figure 2. As in the case of figure 1, in figure 2 we also present the result for the harmonic with n = 0 (i.e., the azimuthally-integrated q T cross JHEP06(2017)017 section). Since we are considering finite values of q T (q T > 2 GeV), the plus-prescription (see eq. (4.5)) is not effective and the harmonic with n = 0 has a dominant behaviour that is proportional to 1/q T ln(M/q T ) at small values of q T : the corresponding numerical results in figure 2 are consistent with such behaviour. The other results in figure 2 regard the n-th harmonics with the specific values n = 1, 2, 4, 6. The numerical results for these azimuthalcorrelation harmonics have similar quantitative features: dσ n /dq T is small and it has a mild dependence on q T in the region where q T ∼ > 20 GeV; at smaller values of q T , the q T shape changes, and the absolute magnitude of dσ n /dq T increases rapidly as q T decreases. This small-q T numerical behaviour is consistent with the dependence dσ n /dq T ∝ 1/q T , and we interpret these results as a clear indication of singular azimuthal harmonics with n = 1, 2, 4, 6. As we have previously discussed, the harmonic with n = 2 has singularities of both collinear and soft origins. The initial-state collinear radiation at the point i) cannot give a singular contribution to the harmonic with n = 1, 4, 6: in the case of these harmonics the singularity is due to the soft wide-angle radiation at the point ii) (radiation at large angles with respect to the directions of the initial-state colliding partons and to the direction of the momentum of the jet j in the Zj system). If dσ n /dq T ∝ 1/q T at small values of q T , the size of the q T slope of dσ n /dq T roughly corresponds to the size of the n-th harmonic of the coefficient a corr (q T ) in eq. (4.5). From the comparison of the results with n = 1, 2, 4, 6 in the small-q T region of figure 2, we notice that the q T slope of dσ n /dq T with n = 1 is sizeable and this is a clear evidence of singular harmonics with odd values of n.

Resummation
The perturbative resummation of singular azimuthal-correlation contributions was considered in refs. [11,12]. Reference [11] presents an all-order treatment (i.e., the treatment is valid to arbitrary logarithmic accuracy) of azimuthal correlations, but it is limited to processes with colourless systems F (the individual particles in F have no QCD colour charge). Therefore, ref. [11] only deals with azimuthal correlations due to initial-state collinear radiation (see the point i) previously mentioned in section 5.1) for production processes that occur through the gluon fusion channel (gg → F + X). The all-order resummation for QQ production in ref. [12] is limited to an explicit treatment up to NNLL accuracy, but it deals with a specific example of the relevant case in which the produced final-state system F is colourful (the individual particles in F carry QCD colour charge). Therefore, the analysis and results of ref. [12] not only deal with collinear radiation, but they also deal with all the main conceptual (and technical) aspects related to azimuthal-correlation effects due to soft wide-angle radiation (see the point ii) previously mentioned in section 5.1). The all-order resummation structure in ref. [12] is quite different from that of azimuthally-insensitive q T cross sections (e.g., the case of the DY process [35]): the former includes resummation factors and helicity/colour correlations that originate from the collinear/soft contributions at the points i) and ii) in section 5.1. The resummation results of ref. [12] certainly have some process-dependent features (due to the specific production of F = {QQ}), but they are generalizable (although with specific differences and, possibly, complications at the technical level) to generic production processes of the type in section 4.1 with final-state colourful systems F (e.g., F = {V j}, {jj}). Therefore, we can consider the all-order re-JHEP06(2017)017 summation for QQ production in ref. [12] as the prototype of resummation for this entire class of processes and, consequently, we can present some general comments (by avoiding cumbersome technicalities) in the following.
The all-order resummation in ref. [12] is performed in b space, by treating azimuthallyindependent and azimuthal-correlations contributions on equal footing. Then the q T cross section is obtained by inverse Fourier transformation of the resummed result from b space to q T space. To facilitate the discussion of the comparison (analogies and differences) with the azimuthally-insensitive case (see section 4.2), it is convenient to consider the q T dependence of the n-th harmonic: The azimuthal asymmetry in eq. (2.11) is obtained by integration of eq. (5.1) over q T . Considering harmonics, the inverse Fourier transformation from b space to q T space can be recast in the form of a Bessel transformation. In particular, the n-th harmonic projects out the nth-order Bessel function J n (x) (see, e.g., section 6 in ref. [11]), and the final result of the resummation procedure of the singular azimuthal-correlation terms has the following (sketchy) form: dσ The resummation formula (5.2) allows us to discuss analogies and differences with respect to the corresponding formula (4.12) for the azimuthally-independent case (which coincides with the n = 0 harmonic). Obviously, we can also consider the n-th harmonic that is obtained through the replacement cos(nϕ) → sin(nϕ) of the weight function in eq. (5.1). Such a replacement does not change the Bessel transformation structure of eq. (5.2), although the factor Σ (res) n has different explicit expressions for cosine and sine harmonics. The factor Σ (res) n (M, b) in eq. (5.2) embodies the effect of the all-order resummation procedure in b space. In the following, we first comment on the perturbative expansion of eq. (5.2) and, then, we discuss the effect of perturbative resummation.
At the lowest perturbative order, the factor Σ (res) n (M, b) in the integrand of eq. (5.2) corresponds to the Fourier transformation of the singular contribution a corr (q T )/q 2 T in eq. (4.5). Therefore, Σ (res) n (M, b) is simply proportional to the n-th harmonic of the azimuthal function a corr (q T ) and its dependence on b is only due to the Bessel transformation of 1/q 2 T . This b dependence is that of the following integral over q T : We note that, having selected the n-th harmonic, the non-integrable singular function 1/q 2 T has a well defined Bessel transformation. We also note that this Bessel transformation is a constant function of b in the large-b region. 13 13 We can introduce an upper bound qT ∼ < M in the integral of eq. (5.3), as it would be appropriate to deal with the small-qT region (approximation). In the large-b limit (bM 1), the result of such integral is that in the right-hand side of eq. (5.3), modulo additive vanishing corrections (corrections that vanish as powers of 1/(bM )).

JHEP06(2017)017
At higher perturbative orders, singular azimuthal correlations produce a dependence of the q T cross section dσ/dM 2 d 2 q T on 1/q 2 T times powers of ln(M 2 /q 2 T ) (with at most two powers of ln(M 2 /q 2 T ) for each additional power of α S ). As discussed in section 4.1, such q T dependence is factorized with respect to the dependence on the azimuthal-correlation angle. Therefore, the b dependence of Σ (res) n (M, b) is due to the Bessel transformation of these singular functions of q T . For example, in the case of the 2nd harmonic the dominant (DL) next-order correction is proportional to the following q T space integral: . is the Euler number) is a numerical coefficient (of kinematical origin) that customarily appears in the context of Fourier (Bessel) transformations of logarithmic functions. In general, the typical singular behaviour of dσ n /dM 2 dq 2 T at small q T leads to contributions to Σ (res) n (M, b) that are proportional to the following basic integrals: The result of the basic integrals is are numerical coefficients that depend on n and m (see eq. (5.9)). Note that each of these integrals I We note that the result in eq. (5.6) for the integrals I [ n ] k can be obtained by using the generating function method that is used in appendix B of ref. [51]. The relevant generating function is with where Γ(z) is the Euler Gamma function. In particular, the coefficient d m in eq. (5.6) is the m-th derivative of d [ n ] (λ) with respect to λ at the point λ = 0: (5.9) We remark again that the singular (and non-integrable) azimuthal-correlation contributions in q T space at each f.o. have a corresponding nth-order Bessel transformation in b space (see eqs. (5.3)-(5.6)). Since the singular terms have an azimuthal dependence JHEP06(2017)017 that is factorized with respect to the q T dependence in q T space, our remark equally applies to the two-dimensional Fourier transformation from q T space to b space of the entire azimuthal-correlation component dσ corr /dM 2 d 2 q T of the cross section (the projection onto the n-th harmonic in q T space and the sum over an infinite set of harmonics in b space are procedures that can be applied with no further conceptual difficulties). The general f.o. structure of Σ (res) n (M, b) in the large-b region is formally similar to that of the corresponding function Σ (res) az.av. (M, b) for the azimuthally-averaged cross section. In both cases we are dealing with logarithmically-enhanced singular terms in q T space at small values of q T (q T M ) and with corresponding powers of ln(bM ) and constant terms in b space at large values of b (bM 1). In both cases the all-order resummation procedure in b space deals with the resummation of ln(bM ) terms, and the formally dominant contributions at each f.o. are of DL type (two additional powers of ln(bM ) for each additional power of α S ). In the case of azimuthal correlations, the ln(bM ) and constant terms originate from the Fourier (Bessel) transformation of non-integrable functions of q T (see eq. (5.5)). This origin can be contrasted at the formal level with that of the analogous contributions to the azimuthally-averaged case, in which the ln(bM ) (and constant) terms arise from the Fourier (Bessel) transformation of terms that are plus-distributions (and δ-function contributions) of q T . The basic integrals that occur in the azimuthally-averaged case are obtained from the right-hand side of eq. (5.5) throughout the replacement J n → J 0 and the replacement of the q T functions with the corresponding plus-distributions. The explicit computation of these integrals for the azimuthally-averaged case is presented in the appendix B of ref. [51] (see, in particular, eqs. (129) and (141) of the arXiv version of ref. [51] or, equivalently, eqs. (B.18) and (B.30) of the corresponding published version).
We now discuss the all-order behaviour of the resummed result in eq. (5.2). The b-space resummed cross section Σ (res) n involves several factors (see section 6 in ref. [11] and section 2 in ref. [12]) and, in particular, the explicit dependence on the n-th harmonic is definitely related to the specific multiparticle system F that is produced in the hard-scattering process. Although the resummed functions Σ (res) az.av. and Σ (res) n in eqs. (4.12) and (5.2) are different, the analysis in refs. [11,12] shows that the dominant DL behaviour of these two b-space cross sections is controlled by the same universal (process-independent) factor, namely, the Sudakov form factor S c (M, b) that we have already recalled in section 4.2. The presence and universality of this common factor, S c (M, b), in the resummed result of both cross sections in eqs. (4.12) and (5.2) follow from the physical origin of the DL contributions to both cross sections: the DL terms originate from radiated partons that are soft and collinear to the initial-state colliding partons c 1 and c 2 of the corresponding main production channels (see eq. (4.3)) in the small-q T region.
As in the case of azimuthally-independent contributions (see our comments in section 4.2), the Sudakov form factor produces a strong suppression of the large-b region (bM 1) in the integrand on the right-hand side of eq. (5.2). This damping effect leads to resummed perturbative predictions for the q T dependence of the n-th azimuthal harmonic dσ n /dM 2 dq 2 T that have a smooth and, especially, integrable behaviour in the small-q T region. The integrable behaviour can be checked by performing the limit q T → 0 of eq. (5.2). Owing to the Sudakov-type suppression of the large-b integration region, we can proceed JHEP06(2017)017 (as in section 4.2) to use the small-q T approximation J n (bq T ) = O((bq T ) n ) in the right-hand side of eq. (5.2), and we obtain [11] dσ (res) n 10) The small-q T behaviour in eq. (5.10) is integrable at q T = 0 for every harmonic (n = 1, 2, 3, . . . ). This is a highly non-trivial result of the all-order resummation procedure of terms that individually behave in a very singular way as 1/q 2 T (modulo powers of ln(M 2 /q 2 T )). In the case of the azimuthally-averaged cross section, as the result of the resummation procedure, the fixed-order singular behaviour is turned into the constant behaviour of eq. (4.13). In the case of the n-th azimuthal harmonic the effect of resummation is even more substantial, since the 1/q 2 T singular behaviour is turned into the suppressed power-like behaviour, q n T , of eq. (5.10). We thus expect that the shape of the resummed q T spectrum, dσ n /dM 2 dq 2 T , for the n-th azimuthal harmonic can be substantially different from the shape of the q T spectrum for the azimuthally-averaged cross section. Independently of the detailed q T shape, the all-order perturbative resummation has a 'dramatic' effect on the total production rate of the azimuthal asymmetries. The total (q T integrated) azimuthal asymmetries turn out to be finite (and computable) after resummation, whereas they can be divergent at fixed perturbative orders. In contrast, the total (q T integrated) azimuthally-averaged cross section is finite order-by-order in perturbation theory, and it is basically not affected by resummation effects in the small-q T region.
As we have discussed in section 4.1, the f.o. divergences of azimuthal asymmetries are produced by primary real emission of one final-state parton, which is absolutely not compensated by virtual radiation at the corresponding perturbative order. At higher perturbative orders, real emission (which is strictly necessary to produce azimuthal asymmetries) is accompanied by ensuing virtual and real radiative effects. Throughout all-order resummation, virtual radiative effects turn out to be dominant in the small-q T region, and they produce the Sudakov form factor suppression that makes the azimuthal asymmetries finite (and with a smooth q T shape at small values of q T ).
The small-q T behaviour in eqs. (4.13) and (5.10) is the result of a combined effect of QCD dynamics and kinematics. We conclude this subsection with a brief discussion of this statement.
We note that the q T dependence in eqs. (4.13) and (5.10) is eventually driven by the corresponding small-q T behaviour of the Bessel functions J n (bq T ) in the right-hand side of eqs. (4.12) and (5.2). The presence of the Bessel function J n in eqs. (4.12) and (5.2) has a kinematical origin: the Bessel functions simply arise from the two-dimensional Fourier transformation from b space to q T space after projection onto the n-th harmonic component. This kinematical origin, however, does not imply that the small-q T behaviour in eqs. (4.13) and (5.10) is purely produced by (all-order) kinematical effects. Such small-q T behaviour has indeed a dynamical origin, since it is due to the strong Sudakov-type suppression of the large-b integration region (see eqs. (4.12) and (5.2)) that is embodied in the resummed components Σ (res) az.av. and Σ (res) n . To make clear the role of this strong suppression, we can consider an explicit counterexample as given by a b space function Σ n (M, b) that is is power-like suppressed by the factor (1/bM ) 2λ in the large-b region 14 (we assume that λ is positive, and we may also assume that λ is arbitrarily large if we want to obtain a suppression effect that is quantitatively strong). The Bessel transformation (see the righthand side of eq. (5.2)) of such a power-like function leads 15 to a q T differential cross section dσ n /dM 2 dq 2 T that is equally power-like suppressed as (q 2 T ) λ−1 in the small-q T region, with no power-like dependence on the order n of the Bessel function J n . The Sudakov-type suppression in b space acts differently on the azimuthally-averaged cross section and on the azimuthal asymmetries in q T space. In the case of the azimuthally-averaged cross section dσ (res) az.av. /dM 2 dq 2 T , resummation leads to a smearing of the LO δ-function behaviour (see eq. (4.4)), and the resummed q T distribution broadens at small values of q T . In the case of the azimuthal asymmetry dσ (res) n /dM 2 dq 2 T , resummation leads to a suppression of the 1/q 2 T behaviour at the lowest order, and the resummed q T distribution is power-like suppressed at small values of q T . In both cases, all-order kinematical effects have a nonnegligible role. Transverse-momentum conservation, which is 'exactly' (provided q T M ) implemented through b space resummation, is eventually responsible for the small-q T behaviour in eqs. (4.12) and (5.2), which is not a DL exponentially suppressed q T behaviour. This feature is (partly) a consequence of the two-dimensional nature of the kinematical conservation law of the transverse-momentum vector q T , which is the relevant kinematical variable for the small-q T behaviour of the azimuthally-averaged cross section and of the azimuthal asymmetries. In contrast, we can recall the effect of resummation on observables, such as event shapes, that are controlled by kinematical variables (e.g., energies and invariant masses) with one-dimensional conservation laws: near the exclusive boundary of the phase space, these observables have a DL exponential suppression (see, e.g., ref. [52]), a suppression effect that is stronger than that produced by any fixed power of the kinematical distance from the phase space boundary.
As we have discussed in section 3, QED radiative corrections can also produce singular azimuthal asymmetries (for instance, in the DY process qq → Z(e + e − ) + X). The QCD resummation procedure that we have discussed in this subsection can equally be applied to QED radiative corrections, and it conceptually leads to the same results and resummation effects that we have just discussed for the QCD case. The same reasoning and conclusions are basically valid in a pure QED context, such as, for instance, in the QED computation of e + e − → µ + µ − + X (although the resummed b-space expression and the QED Sudakov form factor acquire an explicit dependence on the mass of the colliding leptons in the initial state).

Azimuthal asymmetries in tt production
To illustrate resummed results on the q T dependence of azimuthal asymmetries at the quantitative level, we choose the specific process of tt production. Our choice has several reasons, besides the fact that tt production is certainly a process of high phenomenological relevance. One reason is that for tt production we have explicit theoretical control [12] of 14 This is exactly the large-b behaviour that is produced by the small-qT non-singular terms (i.e., the terms that are simply denoted by the dots on the right-hand side of eq. (4.5)). 15 The reader can easily check this result, for instance, by inspection of the Bessel integral in eq. (5.7).

JHEP06(2017)017
singular azimuthal correlations and their QCD resummation. Moreover, due to the nonvanishing QCD colour charge of the produced t andt, in this process we are dealing with singular azimuthal-correlation effects of both collinear and soft origin (see the points i) and ii) in section 5.1). These effects appear in the QCD radiative corrections to both production channels of quark-antiquark annihilation (qq → tt + X) and gluon fusion (gg → tt + X).
Finally, singular azimuthal correlations first arise at the NLO level with respect to the computation of the tt total cross section and, therefore, they are expected to produce effects that have a relatively-large size.
Perturbative resummation of singular azimuthal correlations is under theoretical control [11] also in the case of production processes of colourless systems F (e.g., F = {γγ}, {Zγ}, {ZZ}, {W + W − }). However, in these cases we are sensitive only to singular contributions of collinear origin (see the point i) in section 5.1) and not to those of soft origin (see the point ii) in section 5.1). These singular contributions of collinear origin arise only from QCD radiative corrections to the sole gluon fusion production channel gg → F + X.
The colourless system F is typically produced at the LO throughout the quark-antiquark production channel (qq → F ). Owing to the absence of direct coupling of the gluons to the colourless particles in F , the production channel gg → F is suppressed by some powers of α S with respect to the channel qq → F : typically, gg → F enters at the NNLO in the computation of the total cross section, and the singular azimuthal asymmetries due to gg → F + X first appear at the N 3 LO. Therefore, there is a formal mismatch (suppression) of two powers 16 of α S of the expected relative size of singular azimuthal-correlation effects (which turn out to be finite after QCD resummation) in going from tt production to the production of a colourless system F . Owing to all the reasons that we have just discussed, we think it is more 'interesting' to focus our illustrative quantitative analysis on the case of tt production.
We consider tt production in pp collisions with the centre-of-mass energy √ s = 8 TeV at the LHC. In section 3 (see figure 1-right) we have considered the computation of cos(nϕ) harmonics, where ϕ is the azimuthal angle of the top quark in the CS frame (the differential cross section is integrated over the polar angle of the top quark and over the rapidity y and invariant mass M of the tt pair). The n-th harmonics in figure 1-right are integrated over the transverse momentum q T of the tt pair (with the lower limit q T > q cut ). In this subsection we present numerical results on the q T dependence of the n-th harmonic, and we explicitly consider dσ n /dq T rather than dσ n /dq 2 T . Since dσ n /dq T = 2q T dσ n /dq 2 T , the presence of singular harmonics implies that dσ n /dq T ∝ 1/q T (modulo logs) as q T → 0 in f.o. perturbative calculations.
We specifically consider the cos(2ϕ) harmonic. We present results at the lowest perturbative order (figure 3-left) and including all-order resummation (figure 3-right). 16 This formal conclusion (based on counting the powers of αS) has a caveat, since the channel gg → F can receive a quantitative enhancement with respect to the channel qq → F from the possibly large luminosity of the gluon PDF. In particular, the gluon PDF can be larger than the antiquark PDF of the proton. Roughly speaking, the gluon PDF enhancement can quantitatively compensate the effect of a single power of αS (see, e.g., the results in ref. [53] for the diphoton, F = {γγ}, production process).

JHEP06(2017)017
(a) (b) Figure 3. Transverse-momentum spectrum dσ n /dq T for the n = 2 harmonic in tt production. In the left panel the exact NLO result (blue) is compared to the 1/q 2 T singular component (black), which is split in its collinear (black dashed) and soft (blue dashed) parts. The singular components are plotted with opposite overall sign. The finite result (NLO-sing) that is obtained by subtracting the singular component from the exact NLO result is also shown (red). In the right panel the NLO result (blue) is compared to the purely resummed result (red) and to the complete matched result (black), computed as described in the text.
The set-up and parameters of the lowest-order calculation are the same as those described and used in section 3. Specifically, we compute the NLO contribution to tt production that is due to all tree-level partonic processes whose final state is 'tt + 1 parton'. The result of the q T dependence of the cos(2ϕ) asymmetry at NLO is shown in figure 3-left: the cross section increases by decreasing q T and, in the small-q T region, the increasing behaviour is consistent with the 1/q T behaviour predicted in ref. [12]. The NLO result in the first q T -bin (which includes the region where q T = 0) is not shown in figure 3-left since it would be divergent to +∞. The singular 1/q T behaviour predicted in ref. [12] in analytic form can be evaluated numerically, and the corresponding result is also presented in figure 3-left. Actually, in figure 3-left we plot the NLO singular result with the opposite overall sign (we do that for better visibility, since this avoids the presentation of close histograms). We remark that the NLO singular result in figure 3-left exactly behaves as 1/q T (i.e., it corresponds to the 2nd harmonic of the term a corr (q T )/q 2 T in eq. (4.5)), with no additional q T dependence due to NLO contributions that are non-singular in the limit q T → 0 (i.e., the terms denoted by dots in the right-hand side of eq. (4.5)). In figure 3left we also present the numerical result (it is denoted by the label 'NLO − sing') that is obtained by subtracting the 1/q T singular part from the complete NLO result. Therefore, 'NLO − sing' is expected (predicted) [12] to be finite in the limit q T → 0 and, moreover, from it we can numerically extract the behaviour of the dominant non-singular contribution to the NLO asymmetry in the small-q T region. By direct inspection of figure 3-left we see that 'NLO − sing' is indeed finite: in the limit q T → 0 the finite result appears to be numerically consistent with a non-vanishing (actually, negative) constant value. This

JHEP06(2017)017
implies that, in the case of the cos(2ϕ) asymmetry, the dominant non-singular behaviour of dσ n /dq T is proportional to a constant. 17 As we have discussed, singular azimuthal correlations can have collinear and soft origins (see the points i) and ii) in section 5.1). In the case of tt production, the collinear and soft contributions were separately evaluated in ref. [12]. In figure 3-left we report the numerical results of the splitting of the complete NLO singular contribution in its collinear and soft parts. We note that the collinear and soft parts have the same relative sign, so that there is no numerical compensation between them. Both the collinear and soft parts behave as 1/q T : removing one of these two parts from the singular contribution would lead to a divergent result for the combination 'NLO − sing' in the limit q T → 0. The overall size of the soft part is smaller than that of the collinear part, however, the soft part is definitely not negligible: the ratio between the soft and collinear parts is approximately 0.3. Moreover, we remark that the relative size of collinear and soft contributions depends on the type of colliding hadrons and on the centre-of-mass energy of the collision. 18 We also comment on the results of figure 3-left in the large-q T region. We note that 'NLO − sing' is negative, and this is what we can generally expect at sufficiently large values of q T . The NLO calculation is expected to be physically well behaved at large q T and, hence, the NLO result should be suppressed stronger than 1/q T . On the contrary, according to our definition, the small-q T singular approximation behaves exactly as 1/q T over the entire q T range. As a consequence, at large values of q T the singular approximation overshoots the NLO result, and the result of 'NLO − sing' turns out to be negative.
We now move to consider resummation for the q T dependence of the cos(2ϕ) asymmetry. Our quantitative results are presented in figure 3-right. Before discussing the results in figure 3-right, we briefly comment on the actual content of our resummed calculation. Details on the calculation will be presented in a forthcoming paper, in which we also consider resummation results for the customary azimuthally-integrated q T cross section in tt production.
We use the resummed theoretical results of ref. [12] for tt production by applying the implementation formalism of ref. [51]. This formalism has so far been applied to transversemomentum resummation of azimuthally-integrated (or azimuthally-insensitive) cross sec- 17 In the case of the DY process, QCD radiative corrections do not produce singular harmonics in the limit qT → 0. At the lowest perturbative order the small-qT behaviour of the cos(2ϕ) harmonic is dσDY,n=2/dqT ∝ qT ln(M/qT ) [54] (i.e., we have [dσ]2 ∝ ln(M/qT ), where [dσ]2 is the qT cross section in eq. (3.1)), and the cos(2ϕ) asymmetry vanishes in the limit qT → 0. We note that the dominant non-singular term for the cos(2ϕ) asymmetry in tt production is enhanced by a power of 1/qT (modulo logs) with respect to the DY process. 18 This remark follows from the following point. The soft contribution affects both the gluon fusion (gg → tt + X) and qq annihilation (qq → tt + X) channels, whereas the collinear contribution only occurs through the gluon fusion channel. These two production channels contribute to the NLO result with a relative size that depends on the gg and qq PDF luminosities. In the case of tt production in pp collisions at √ s = 8 TeV, the qq luminosity is much smaller than the gg luminosity. However, by decreasing the relative size of the gg luminosity, the relative effect of the soft (vs. collinear) contribution would increase. For instance, the soft/collinear ratio increases in going from pp collisions to pp collisions at the same or smaller (e.g., Tevatron) centre-of-mass energy. Similar PDF related considerations about the relative effect of collinear and soft contributions are valid for production processes of other high-mass systems F .

JHEP06(2017)017
tions for many production processes (see, e.g., ref. [55] and related references therein). The complete result dσ for the azimuthal-correlation component of the q T differential cross section is obtained by using a matching procedure that combines the resummation of the singular contributions (dσ (res) ) with the evaluation of non-singular (finite) contributions (dσ (fin) ) at a given f.o. in perturbation theory. We write dσ = dσ (res) + dσ (fin) , and we work at NLL+NLO accuracy for the resummed (NLL) and f.o. (NLO) contributions. Therefore, the f.o. part dσ (fin) that contributes to the complete matched (res + fin) result dσ in figure 3right is exactly the 'NLO − sing' contribution in figure 3-left. The resummed part dσ (res) is first computed in b space, where we are dealing with the resummed factor Σ (res) (M, b). The resummed factor Σ (res) (M, b) is proportional to the contribution from the singular part of the NLO result (see figure 3-left) and, very roughly speaking, it is multiplied by a form factor exp{G(M, b)} that is evaluated in resummed form up to NLL accuracy. We mean that in the large-b region (bM (M, b), where G LL (M, b) and G N LL (M, b) embody the resummation ( ∞ k=1 ) of the LL terms α k S ln k+1 (bM ) and the NLL terms α k S ln k (bM ), respectively. We remark an important difference between G LL and G N LL . The LL contribution 19 G LL (M, b) is essentially process independent and it basically depends on the initial-state colliding partons of the partonic subprocesses 20 gg → tt + X and qq → tt + X. By contrast, G N LL (M, b) includes an important amount of processdependent information. In particular, it includes the effect of soft wide-angle radiation for the tt production process. This effect is controlled by a colour-dependent anomalous dimension, namely, the soft anomalous dimension Γ (1) t in ref. [12] (see eqs. (16), (17) and (33) in ref. [12]), and we resum the effect in complete exponentiated form throughout the diagonalization of Γ (1) t in colour space. At the technical level, we note that we compute Σ (res) (M, b) in two-dimensional b space and we obtain the q T cross section by numerically performing the Fourier transformation to q T space. The computation of the cos(nϕ) harmonic is also performed by numerically integrating over ϕ. Therefore, we do not perform the projection onto the n-th harmonic in analytic form (see eq. (5.2)), and our implementation can be directly applied to the numerical evaluation of n-th harmonics with different values of n.
Our quantitative resummed results on the q T dependence, dσ/dq T , of the cos(2ϕ) harmonic are presented in figure 3-right. We show the purely resummed result (dσ (res) ) and the complete final result after the inclusion of the matching term (dσ (fin) ) at NLO. In figure 3-right we also report the NLO result (see figure 3-left) for direct comparison with the resummed results. The error bars in the complete matched result denote the numerical errors of our calculation (for simplicity of presentation the numerical error bars are not reported in the other histograms of figure 3). By direct inspection, we observe that the purely resummed result quantitatively differs from the complete matched result, and we comment about that below. Since we are considering the results after integration over the 19 If we restrict the evaluation of G LL to DL accuracy (i.e., only k = 1), exp{G LL } is exactly the DL approximation of the Sudakov form factor that we have mentioned in sections 4.2 and 5.2. 20 In the gluon fusion channel G LL (M, b) is the same as the corresponding contribution to Higgs boson production [51,56]. Analogously, in the qq annihilation channel G LL (M, b) is the same as the corresponding contribution to the DY process [55].

JHEP06(2017)017
invariant mass M of the tt pair, the characteristic hard scale of the q T cross section is of the order of 2m t (m t is the mass of the top quark), which is the minimum value of M , and the q T range that is shown in figure 3-right includes the regions of small and intermediate values of q T . The high-q T region, q T ∼ > 2m t , is not shown in figure 3-right. We first comment on the results in the large-q T region of figure 3-right. We observe that, starting from intermediate values of q T (e.g., q T ∼ > 150 GeV), the complete resummed result tends to nicely agree with the NLO result. Note that we are not using a direct (though, possibly smooth) switching procedure from the purely resummed result to the NLO result at some intermediate value of q T . The agreement between the complete result and the NLO result at large q T is thus a non-trivial consequence of the matching procedure. In the absence of the contribution from the matching term (dσ (fin) ), the purely resummed result sizeably overshoots the NLO result at intermediate and large values of q T (for instance, at q T ∼ 350 GeV the purely resummed result is about 8 times larger than the NLO result). This feature is due to the fact that the size of dσ (res) is driven by the small-q T singular part of the NLO calculation, and the singular part is not a good quantitative approximation of the NLO calculation at intermediate and large values of q T (we have already discussed that and the corresponding overshooting effect in our comments about figure 3-left). In the complete (matched) result, the effect of this approximation is compensated by adding dσ (fin) (i.e., the 'NLO − sing' result of figure 3-left), which sizeably decreases dσ (res) at large q T . Once dσ (res) and dσ (fin) are combined, the difference between the resummed and NLO results is much reduced at intermediate values of q T , and the higher-order (i.e., beyond NLO) resummation terms that are included in our resummed calculation do not have a large residual effect in that q T region.
We now comment on the results in the small-q T region of figure 3-right. In this region the resummed results unavoidably differ from the NLO result: the NLO result is much larger at small values of q T , and it diverges to +∞ if q T → 0. According to our previous discussion (see eq. (5.10)), the purely resummed result behaves as dσ/dq T ∼ q 3 T in the limit q T → 0, and the numerical results in figure 3-right are consistent with this expectation. The maximum value of the resummed cross section is in the region (q T bin) where 30 GeV< q T <60 GeV. Using our implementation of the theoretical results of ref. [12], we have performed a corresponding resummed calculation of the customary (azimuthallyintegrated) q T cross section dσ/dq T of the tt pair and we find [57] that its q T shape is quite different from that of the cos(2ϕ) harmonic in figure 3-right (as expected from the general comments below eq. (5.10)): dσ az.av. /dq T has a maximum at q T ∼ 15 GeV and a much softer q T spectrum in the small-q T region. The qualitative features of the complete result for the 2nd harmonic are similar to those of the corresponding purely resummed result. However, the inclusion of the non-singular NLO contribution (dσ (fin) ) has a relevant quantitative effect, not only at intermediate values of q T (as previously discussed) but also in the small-q T region (for instance, in the q T bin where 30 GeV< q T <60 GeV, the purely resummed result turns out to be reduced by about 40%). We have also checked the quantitative effect of the NLL resummation terms that are due to the inclusion of the soft anomalous dimension Γ (1) t . We find that also Γ (1) t has a substantial effect in the small-q T region: dσ/dq T would be quantitatively smaller and with a harder q T shape by removing the effect of Γ

JHEP06(2017)017
By integrating dσ/dq T for the cos(2ϕ) harmonic over the entire kinematical region of q T , we find the total value σ tt,n=2 = 3.0 pb. In comparison, the value of the tt total cross section (σ N LO tt = 226 pb) is about 70 times larger than σ tt,n=2 . We remark that our quantitative results for σ tt,n=2 and also for the q T dependence of the cos(2ϕ) harmonic have to be regarded as 'effective' lowest-order predictions within (resummed) perturbative QCD. We do not attempt to quantify the theoretical accuracy of these predictions. As is well known, customary procedures (e.g., studies of dependence on factorization and renormalization scales) that can be applied to estimate theoretical uncertainties of lowest-order predictions typically fail in reproducing the quantitative effect of higher-order contributions (this is especially true in the case of QCD calculations for hard-scattering processes in hadron collisions at high energies). A first quantitative estimate of the theoretical uncertainty of the lowest-order predictions requires a computation at the subsequent perturbative order. Once the results at two subsequent orders are known, their relative difference can be used for a quantitative estimate of the theoretical uncertainty. At the present lowest-order level, we limit ourselves to adding few comments. The bulk of the contribution to the cos(2ϕ) harmonic originates from the small and intermediate regions of q T ; we find that the high-q T region (say, q T ∼ > 2m t ) gives only a few percent contribution to the total value of σ tt,n=2 . The NLL resummed contributions due to the soft anomalous dimension Γ (1) t are not negligible: removing the effect due to Γ (1) t in our calculation, we find that the total value of σ tt,n=2 decreases by about 30%. The total value of σ tt,n=2 would be much larger by considering the purely resummed term (dσ (res) ) and no matching with the complete NLO calculation: removing the matching contribution (dσ (fin) ), the value of σ tt,n=2 would increase by roughly a factor of two. This highlights the importance of non-singular NLO effects at small and intermediate values of q T .
In our resummed calculation we numerically perform the Fourier transformation from b-space to q T -space, and we obtain the q T cross section dσ/dq 2 T dϕ. We can then numerically evaluate different harmonics. We have carried out an exploratory numerical computation of the cos(4ϕ) asymmetry for tt production, and we find results with the expected features, analogously to the cos(2ϕ) asymmetry. Nonetheless, the size of the cos(4ϕ) asymmetry is very small and we do not present corresponding resummed results. As shown in figure 1-right, the n = 1 azimuthal asymmetry for tt production is not divergent if computed at its lowest perturbative order, and the corresponding q T spectrum is not singular. In view of our discussion in section 5.1, the next-order QCD correction to the n = 1 harmonic can produce a singular q T spectrum, so that the corresponding QCD prediction can require a resummed calculation at a subdominant logarithmic accuracy with respect to that considered in ref. [12].
We add some general comments on the contribution of non-singular terms (e.g., the terms denoted by dots in the right-hand side of eq. (4.5)) to azimuthal correlations at low values of q T . In the case of the azimuthally-averaged cross section dσ az.av. /dq 2 T , the perturbative resummation of singular terms produces a constant behaviour in the small-q T region (see eq. (4.13)) and, typically (e.g., in absence of additional kinematical cuts), the non-singular terms also have a constant behaviour at low values of q T (order-by-order in JHEP06(2017)017 the perturbative expansion the ratio between non-singular and singular terms if formally of O(q 2 T /M 2 ), modulo logarithms, in the limit q T → 0). As a consequence, resummed and f.o. non-singular contributions to dσ az.av. /dq 2 T behave similarly at low values of q T , although the non-singular contributions are perturbatively (and quantitatively) suppressed. In the case of the cross section dσ n /dq 2 T for the n-th harmonic, the resummation of the singular terms produces an enhanced power suppression (∝ q n T ) of dynamical origin (see eq. (5.10)) in the small-q T region and, therefore, non-singular terms are eventually more relevant at low values of q T if they are treated at f.o. in perturbation theory. We note that this reasoning neglects the fact that non-singular terms can also have logarithmic enhancement (powers of ln(M/q T )) order-by-order in perturbation theory and, consequently, an appropriate resummation treatment of non-singular terms can be required. For instance, QCD resummation of non-singular terms has been investigated (see ref. [58] and references therein) for the DY process. Note, however, that non-singular corrections to singular azimuthal correlations are quite different from azimuthal-correlation effects in the DY process 21 and they are expected to have a high degree of process dependence. We also recall that the non-singular terms produce azimuthal-correlation effects whose actual size and q T behaviour depend on the specific definition of the azimuthal-correlation angles (see eq. (2.8) and related comments in section 2). In particular, the q T dependence of non-singular terms can be redistributed among different n-th harmonics by considering azimuthal correlations that refer to different specifications of the azimuthal angle.
As for non-singular terms, at present we do not have a detailed theoretical understanding that goes beyond their treatment at f.o. . A better understanding of non-singular terms is certainly relevant to have an accurate quantitative control on the detailed q T shape of azimuthal correlations in the region of very-low values of q T . We note, however, that nonsingular terms have a relatively-mild effect on q T integrated azimuthal correlations. For instance, since non-singular terms are integrable in the limit q T → 0, their perturbative logarithmic enhancement still gives contributions to the total (q T integrated) harmonic that are parametrically controlled by powers 22 of α S . Moreover, after resummation of the singular terms (and, basically, as a consequence of the behaviour in eq. (5.10)), the bulk of the contribution to dσ n /dq 2 T is located at relatively-large values of q T (the bulk of the contribution to dσ az.av. /dq 2 T is typically located at smaller values of q T ) and, also in this region the logarithmic enhancement of non-singular terms is expected to produce relatively-mild effects (approximately, also these effects are parametrically controlled by powers of α S ). In summary, we think that an improved (beyond f.o. perturbation theory) treatment of nonsingular terms is relevant especially to determine the detailed (point-by point) q T shape of dσ n /dq 2 T at low values of q T , while even the total effect of non-singular terms in the low-q T region (e.g., in the lower q T -bin, q T < 30 GeV, of figure 3) can be relatively well under control within a f.o. treatment. 21 For instance, in a previous footnote we have already remarked that in the small-qT region the 2nd harmonic cross section dσn=2/dq 2 T has non-singular terms of O(M/qT ) for tt production and of O((qT /M ) 0 ) = O(1) for the DY process. 22 This effect is analogous to the corresponding effect on the total (qT integrated) azimuthally-averaged cross section: its value is reliably computable at f.o. independently of the presence of integrable logarithmically-enhanced, and even singular (plus-distributions), contributions in the small-qT region.

JHEP06(2017)017
We note that, at very low values of q T , also non-perturbative QCD effects can be relevant. In the case of the resummed contribution to the q T cross section, a customary (though much simplified) procedure to model non-perturbative effects consists in supplementing the Sudakov form factor with a non-perturbative form factor that mostly contributes at large values of b in the b-space resummation formula (see eq. (4.12)). A typical parametrization of the non-perturbative form factor is e −g b 2 (g being a parameter whose size, g ∼ O(1 GeV 2 ), is determined by effects at non-perturbative scales), and it produces a strong damping effect a very large values of b. In the resummation formula these non-perturbative effects act in a b-region that is already strongly suppressed by the Sudakov form factor and, consequently (after Bessel transformation from b space to q T space), they mainly affect only the region of low values of q T [30,47]. Such a model of non-perturbative effects can also be applied to the resummed contribution to the n-th harmonic (see eq. (5.2)), and the non-perturbative parameter can also depend on the harmonic (e.g., g → g n ). We may try to apply a related model also to non-singular terms. We consider (part of) non-singular terms of the n-th harmonic, we transform them to b-space and then we supplement them with a non-perturbative form factor e −gn b 2 . The non-singular terms in b-space are powersuppressed at large values of b, so that the non-perturbative form factor still acts on a b-region that is dynamically suppressed (correspondingly, the non-perturbative form factor mainly affects the region of low values of q T ). Nonetheless, at very large values of b, the exponential suppression that is produced by the non-perturbative form factor is so strong that, in the limit q T → 0, dσ n /dq 2 T behaves 23 as q n T , and the non-singular terms (those that are treated in this way) eventually produce the same q T behaviour as the resummed contribution to dσ n /dq 2 T . We cannot argue that such a model is physically justified (such a non-perturbative treatment of non-singular terms is certainly not applicable to azimuthal correlations in a pure QED context). We have mentioned it only to notice that in the region of low values of q T non-singular terms can be affected by both perturbative logarithmic enhancement and non-perturbative contributions. A study of non-perturbative effects in azimuthal-correlation cross sections is definitely beyond the scope of this paper.

Summary
In this paper we have considered high-mass systems formed by two or more particles that are produced in hadron collisions, and we have discussed the azimuthal correlations between the system and one of its particles. We refer to particles in an extended sense, namely, pointlike particles and QCD jets. Our main findings can be summarized as follows.
Despite the infrared-safe nature of the azimuthal distribution, the f.o. QCD computation of this observable can be divergent starting from some perturbative order. This conclusion holds for a large class of processes, including tt, V j, jj, γγ, ZZ, W + W − production, while the corresponding QCD calculation for the DY process is finite to any (arbitrary) perturbative order. We have supported this conclusion by performing QCD calculations at the first non-trivial order for the DY, Z+jet and tt production processes.