Resummation of Boson-Jet Correlation at Hadron Colliders

We perform a precise calculation of the transverse momentum ($\vec{q}_T$) distribution of the boson+jet system in boson production events. The boson can be either a photon, $W$, $Z$ or Higgs boson with mass $m_V$, and $\vec{q}_T$ is the sum of the transverse momenta of the boson and the leading jet with magnitude $q_T=|\vec q_T|$. Using renormalization group techniques and soft-collinear effective theory, we resum logarithms $\log(Q/q_T)$ and $\log R$ at next-to-leading logarithmic accuracy including the non-global logarithms, where $Q$ and $R$ are respectively the hard scattering energy and the radius of the jet. Specifically, we investigate two scenarios of $p^J_T \lesssim m_V$ or $p^J_T \gtrsim m_V$ in $Z$+jet events, and we examine the $q_T$ distributions with different jet radii and study the effect of non-global logarithms. In the end we compare our theoretical calculations with Monte Carlo simulations and data from the LHC.


Introduction
The production of photons, W 's, Z's and Higgs bosons are important processes which allow us to test the Standard Model and extract its fundamental parameters. With precise calculations of the cross sections, they also give opportunities to search for physics beyond the Standard Model through deviations from Standard Model predictions. At hadron colliders, initial state radiation caused by the strong interaction contributes and necessarily affects the boson distributions. Moreover, energetic jets can be produced in association with the boson production. Therefore the understanding of the boson distribution is always complicated by the presence of hadronic activities in the events, which are governed by quantum chromodynamics (QCD). Using fixed-order calculations in perturbative QCD one can systematically improve the description of hadronic radiation. However, in certain regimes the fixed-order perturbative expansion diverges so that an all-order resummation is necessary for the validity of theoretical predictions. This happens when characteristic energy scales relevant in the process become hierarchical so that large logarithms of scale ratios can spoil the validity of fixed-order calculations. The transverse momentum p T distribution of the lepton pair in the Drell-Yan process is a classic example which requires the resummation of log(M V /p T ) in the regime of p T m V , where M V is the vector boson mass. This can be achieved to all-orders using the standard formalism by Collins, Soper and Sterman (CSS) [1]. Alternatively, the logarithms can also be resummed using renormalization group (RG) techniques and soft-collinear effective theory (SCET) [2][3][4][5] (see [6] for a review) as discussed in [7][8][9][10][11][12][13][14][15]. More recently, the resummation has been performed at next-to-next-to-next-to-leading logarithmic accuracy [16][17][18][19].
In this paper, we study the situation in which the boson has significant transverse momentum recoiling against hadronic activities consisting of jets in the final states. Specifically, we consider the q T distribution of the boson and the leading jet system where q T is the vector sum of the transverse momenta of the two objects and q T = | q T |, as illustrated in figure 1. In the 2 → 2 scattering, boson+jet back-to-back limit the value of q T is zero, although the boson and the jet can have large transverse momenta. In the small q T regime, the soft and collinear emissions induce large logarithms of log(Q/q T ) which need to be resummed, where Q represents the hard scattering energy. The situation is similar in di-jet production where q T is defined as the vector sum of the transverse momenta of the two leading jets, and the resummation of log(Q/q T ) at next-to-leading logarithmic (NLL) accuracy without non-global logarithms (NGLs) [20,21] was carried out using the CSS formalism in [22,23]. Similarly, the log(Q/q T ) resummation at NLL level was also performed for photon(γ)+jet [24], Z+jet [25] and top quark+jet production [26,27]. More recently, the NLL resummation in γ+jet production was also carried out using SCET [28].
The purpose of this paper is to derive an all-order expression in SCET for the systematic resummation of log(Q/q T ) in boson+jet production at small q T , including the resummation of the associated jet radius logarithms log R as well as NGLs 1 . The precise understanding of this observable in proton-proton collisions then forms the baseline of such hard probes in nucleus-nucleus collisions where a hot and dense QCD medium called the Figure 1. Boson+jet production in hadron collisions. Here p V and p J are the momenta of the color singlet boson and the jet, and R is the jet radius. By definition q T = p J T + p V T . The modes relevant for the observable q T include the soft modes with momentum p s , and the collinear modes along the two beam directions (n 1 and n 2 ) and the jet direction (n J ). Small-angle soft modes are taken as an independent degree of freedom from those emitted from the jet at wide angle, and its momentum is denoted as p t . The n 1 -collinear and n 2 -collinear modes and soft modes all have a transverse momentum ∼ q T , while the n J -collinear modes carry most of the jet momentum. quark-gluon plasma (QGP) is produced. Through interactions with the medium, jets in the event can be significantly modified while the color-singlet boson remains intact that can serve as a robust reference of the hard scattering process. This makes boson+jet production a useful channel for studying the properties of QGP though the relation between transverse momentum broadening and energy loss of jets in high-energy nuclear collisions [45], which requires a proper resummation of large logarithms [24,46,47]. The kinematic information of the boson+jet system has been explored quite extensively [48][49][50][51][52][53][54]. For example, the q T , the boson-jet momentum imbalance X JV ≡ p J T /p V T , and the azimuthal angle decorrelation |∆φ JV |: the azimuthal angle between the jet and the boson as measured along the beam direction, have been experimentally studied in Z+jet [55][56][57][58][59] and γ+jet [60] events at the LHC.
The rest of the paper is organized as follows. In section 2, we analyze all the relevant degrees of freedom which contribute to q T . We give a detailed derivation of our factorized expression (2.27) using a two-step matching procedure in SCET. In section 3, we discuss the renormalization of all the bare functions entering (2.27) and give an all-order resummation formula in (3.13). We explain the relation between our resummation formula with those in [24,25,28]. The anomalous dimensions relevant for the NLL resummation are also given in this section. In section 4 we analyze the Sudakov double logarithms, while in section 5.2 we perform the resummation of log(Q/q T ) at NLL accuracy for Z+jet production, including log R and NGL resummation. In section 6 we summarize and discuss some intriguing issues for future studies. In appendices A and B, we list the tree-level amplitudes of partonic V +jet production and the anomalous dimensions used in this calculation. In appendix C we give the LO singular terms for q T distribution.

Factorized Expression for q T in Boson+Jet Production
We derive a factorized expression of the differential cross section dσ/d 2 q T for the process where X stands for all the produced particles in the event except for the boson and the particles of the leading jet. As defined previously, the observable q T is the sum of the transverse momenta of the boson and the leading jet. The boson can be either the W , Z, γ or the Higgs boson. We will first focus on the q T region where Λ QCD q T Q and Q is a hard scattering energy scale depending on the leading-jet transverse momentum p J T (and, for a massive boson, the boson mass m V ). We will discuss the factorization of the cross section in SCET and resum large logarithms using RG techniques.

Degrees of Freedom
For q T Q, the dominant contributions to q T come from soft particles in all out-of-jet directions or collinear particles along the beam directions, with transverse momenta of the order q T . As illustrated in figure 1, such radiation can be either soft with p s ∼ q T , or collinear to the two beam directions n 1 and n 2 with p n 1 T ∼ q T or p n 2 T ∼ q T . In this paper all the calculations are carried out in the small R limit. In this case, the small-angle soft mode along the jet direction can be singled out as an independent degree of freedom [32,34,61,62]. Such soft radiation is sensitive to the jet direction and the jet boundary and will be referred to as the coft mode in the following discussions, cf. [32,34]. While wideangle, soft radiation is only sensitive to the total color charge of the jet, the coft mode can resolve any possible collinear constituents of the jet. If a coft radiation is emitted outside the jet it will contribute to the observable q T . We need to consider the coft mode in order to account for multiple out-of-jet radiation and resum the potentially large logarithms of log R.
The above kinematic analysis shows that the relevant SCET degrees of freedom for the calculation of q T in this process include the following modes as illustrated in figure 1 2 , where λ = q T /Q is the power counting parameter. The auxiliary light-like vectorsn i satisfy n i ·n i = 2 for i = 1, 2 and J and we choosen 1 = n 2 andn 2 = n 1 . As we shall discuss in section 4, Q is the hard scattering scale in the process and it may be parametrically different from p J T if the boson is massive. Here all the momenta p µ = (n i · p,n i · p, p n i ⊥ ) are expressed using light-cone coordinates with light-like vectors n i andn i , as denoted by the subscripts n ini in (2.2). We denote the transverse momenta perpendicular to the two beam directions n 1 and n 2 by the subscript T , and the transverse momenta perpendicular to n J by the subscript ⊥ .

Derivation of the Factorized Expression
We derive the factorized expression using SCET with all the degrees of freedom in (2.2). The derivation is carried out in a two-step procedure similar to the one in [32,34].

Matching of QCD onto an intermediate SCET
The intermediate SCET (more specifically, SCET II [3]), includes n 1 -, n 2 -and n J -collinear fields and one soft gluon field. In this step, the hard mode is integrated out and encoded in the Wilson coefficients C. The local collinear gauge invariance demands that collinear fields along different directions do not directly interact with each other. That is, the hard and collinear modes factorize.
Let us denote collectively the infrared (soft or collinear) particles by X IR and write the differential cross section for the process N 1 (P 1 ) + N 2 (P 2 ) → boson(p V ) + X IR as where |X IR is a product of n i -collinear states |X n i and soft state |X s , dX IR denotes the measure for the n−body relativistically invariant phase space of X IR and y V is the boson rapidity. Generically, the leading-order operators are built out of three collinear fields along the three collinear directions of the beams and the jet, and the effective Hamiltonian H takes the form where {dt} ≡ dt 1 dt 2 dt J is the integration measure and C is the Wilson coefficient. The field [φ n i ] α i a i represents an n i -collinear field carrying a color index a i and a Dirac or Lorentz index α i , and it can be either a collinear quark field or a collinear gluon field, which are given, respectively, by [68] where the n i -collinear covariant derivative is defined as D µ n i ⊥ ≡ ∂ µ ⊥ + gA µ n i ⊥ , and W n i is the n i -collinear Wilson line. The set of fields {χ n i ,χ n i , A µ n i ⊥ } form the building blocks to construct the effective operators.
Since the n 1 -and n 2 -collinear modes are not directly measured and will go along the beams, we need to sum over these collinear states. Using the fact that initial colliding hadrons are color neutral and keeping only the leading contribution inn i · P i , we have where = (4 − d)/2 in dimensional regularization and i = 1, 2 labeling the beam. The factor d i is the dimension of the color representation of the field φ f n i , and P is the projector defined as follows 3 , for quarks and antiquarks, (2.7) The function B f /N i is the beam function of the parton species f [70][71][72] in the x T space, which is the Fourier transform of the transverse-momentum dependent (TMD) parton distribution functions (PDFs). Next, we sum over the n J -collinear particles and perform multipole expansion so that the n J -collinear fields only depend on n J · x. Assuming m n J -collinear partons in the jet, we have where the four-momentum of the i-th collinear particle in the jet p µ J i is expressed in terms of the transverse momentum p J i T , the azimuthal angle φ J i and the pseudo-rapidity η J i of the particle. Also, the jet direction n J = (1, sin φ J / cosh η J , cos φ J / cosh η J , tanh η J ) with η J and φ J respectively the rapidity and the azimuthal angle of the jet. For reasons that will become clear later, we also assume that there are some n J -collinear particles radiated outside the jet with a total momentum p out t . Similar to the discussion of beam functions, a color-neutral jet function J k with the virtuality p 2 J and the parton species k can be defined as 3 For gluon beam functions, another projector needs to be included in the study of, e.g., the Higgs pT distribution in the gg → H 0 production channel [10,69]. However, for the process studied in this paper, one can show that the contribution from this projector vanishes at NLL level. Hence, it is neglected here and in the following sections.
for gluons. (2.10) After decoupling the soft fields from (2.4), we will have the product of three soft Wilson lines. Summing over the states of soft gluons gives The color structure of the soft function is the same as the gauge transformation of the amplitude squared for the process. From the gauge invariance, one can show that the above matrix element is always proportional to the unit color matrix for the processes studied in this paper. That is, Plugging (2.6), (2.9) and (2.11) into (2.3) (with X IR summed over as we have done), we have where the sum runs over all parton species i, j, k = q,q, g. The hard function is identified as with ξ 1 and ξ 2 completely determined by the conservation of the + and − components of the partonic momenta in the basis vectors n 1 andn 1 .

Separating coft modes from n J -collinear modes
In this step, we match the purely collinear theory along the jet direction onto an effective theory where the collinear field is split into two submodes as and, accordingly, We distinguish genuine collinear momenta from the coft ones, and the corresponding momentum scalings are shown in (2.2). Here the coft field describes low energy radiation which can resolve the substructure of the jet, and it is emitted from one of the collinear partons in the jet at an angle θ R. In effective theory languague one can take such coft radiation as being an independent mode [32,34]. In this effective theory the soft sector is a combination of soft radiation which can not resolve the detailed structure of the three collinear sectors, as well as coft radiation sourced by the collinear constituents of the jet 4 . In the limit q T p J T , the genuine n J -collinear particles are kinematically forbidden to be radiated outside the jet while coft modes are allowed to be either inside or outside the jet. Inside the jet, the contribution from the coft modes to the jet momentum can be neglected. Hence, the m n J -collinear particles introduced in the previous subsection are genunine n J -collinear particles while those outside the jet are coft particles with a total momentum p out t . The separation of the coft modes from the n J -collinear modes modifies the jet function (2.9) by organizing the coft radiation into the coft Wilson lines [32]. In the following discussion we use the notations adopted in [32,34] and write the amplitude squared for the n J -collinear particles as where the compact notation {p J } stands for the set of collinear parton momenta p J i . Then, the collinear Wilson line W n J in the definition of the n J -collinear field in (2.5) is replaced by with the coft Wilson line Un J along then J direction which organizes coft radiation emitted from the n 1 -and n 2 -collinear directions in the small R limit. Again, for brevity the collinear field φ n J dressed with coft radiation along then J -direction due to the replacement in (2.18) is written as φ n J → φ n J Un J . Also, each n J -collinear parton is dressed with a coft Wilson This means that separating out the coft modes is From the above two equations, one finally has where · · · ≡ 1 d J Tr[· · · ] denotes the trace over all the color indices divided by the dimension of the color representation of φ k n J , and ⊗ is a short-hand notation for 22) and the coft function U m takes the form The set of n J -collinear particles is defined by the anti-k t algorithm [74] which is used in jet reconstruction. The phase space constraint imposed by the sequential clustering can be quite complicated. Alternatively, here we require the angle ∆R ij between each pair of collinear particles be smaller than the jet radius R, In the small R limit, the above requirement is equivalent to imposing the following step functions, which collectively is denoted by Θ in ({p J }). The jet algorithm constraint for a coft gluon with momentum p t is then equivalent to a cone jet algorithm since collinear particles are clustered and define the jet direction n J , (2.26) By making the replacement in (2.21), (2.13) then gives the final factorized expression

Resummation of Large Logarithms
In this section, we discuss the renormalization of the bare functions in (2.27) and the resummation of large logarithms by solving the corresponding RG equations. We also calculate the anomalous dimensions relevant for the resummation at NLL level.

Renormalization and Resummation
The cross section is finite in the limit → 0 but all the bare functions in (2.27) are divergent. In this paper, these functions are renormalized in the MS scheme. The divergent pieces of the bare functions are removed by the renormalization constants, and the anomalous dimensions can be calculated from them according to (B.2). Then the resummation of large logarithms can be achieved by solving the RG equations.

Hard function
The Wilson coefficient C in (2.4) is determined order-by-order in perturbation theory by a matching calculation in QCD and in SCET. In dimensional regularization, the ultraviolet (UV) divergence in the Wilson coefficient is identical to the infrared (IR) divergence in the corresponding on-shell amplitudes in perturbative QCD. Hence, the singularities in the hard function can be subtracted by a multiplicative renormalization constant Z H ij→V k . From Z H ij→V k one can calculate the anomalous dimensions of the hard functions and resum large logarithms in µ h /µ by solving the RG equation with the initial condition H ij→V k (ŝ,t, m V , µ h ) at the hard scale µ h ∼ Q calculated in a matching calculation.

Soft function, beam function and collinear anomaly
The calculation of soft and beam functions involves extra complication which is not seen in the calculation of the hard function. Singularities unregularized by dimensional regularization arise in the calculation of these functions. However, such divergences are artificial because the product of the soft and beam functions is in fact finite, which is a result independent of the regulator. In this paper, we regularize such divergences by modifying the phase-space integrals as [75] Note that we have chosen the common factor of n 1 · k in the regulator. Moreover, the scale separation in (2.2) is broken due to loop corrections, and the hard scale shows up in the perturbative calculation of soft and beam functions. This is referred to as the collinear anomaly by the authors of [9] 5 . By refactorizing out the collinear anomaly, the product of beam and soft function can be written as [78] Alternatively, the collinear anomaly can be dealt with using the rapidity RG method [76,77].
where the hard scale dependence is factored out with the exponent F ⊥ only depending on x T and the scale µ. Here b 0 = 2e −γ E and C i is the Casimir operator of the color representation of the parton i. The divergence in the product of the soft and beam functions on the l.h.s. of (3.3) is to be removed by an overall multiplicative renormalization factor, denoted by Z SBB ( x T , µ, ). In order to construct a universal definition of the beam function (at least in the boson+jet processes considered in this paper), we take Z SBB as a product of the renormalization constants of the collinear anomaly Z CA , the soft function Z S , and the beam functions From these renormalization constants one can calculate the corresponding anomalous dimensions. The collinear anomaly exponent function F ⊥ (x T , µ), beam and soft functions B f /N and S ij→V k satisfy the following RG equations, respectively

Jet function, coft function and non-global logarithms
The calculations of jet and coft functions contain NGLs because of the restricted phase space due to jet definition. As discussed in [32,34], the RG running of the jet and coft functions in the factorized expression (2.27) automatically resums both global and nonglobal logarithms.
In the definition of the jet function in (2.22), the energy of the n J -collinear constituents are integrated over, which results in additional singularities. However, such singularities can be cancelled by the jet functions with lower parton multiplicity. Therefore, in general the renormalization constant of jet functions is a matrix [34], which is defined as Similarly, the renormalized coft function is written as where⊗ denotes the integration over the (m − l) additional directions. Note that Z U lm is defined in a reversed way as opposed to the other renormalization constants. By the RG invariance of the physical cross section, the renormalization matrix of the coft function satisfies In [34,35] one of the authors has explicitly verified that this matrix satisfies a renormalization group equation at two-loop level for non-global jet observables in electron-positron collisions.
For the coft function we specifically extract the global renormalization constant Z U which removes the divergence in the coft function U 1 , Accordingly, we define the non-global renormalization constant aŝ by separating out the global contribution. From (3.8), Z J lm can hence be expressed as the product of non-global and global renormalization constants where we also introduce a global renormalization constant Z J for the jet function. The evolution equations that resum both the global and non-global logarithms in the jet and coft functions can be obtained from (3.6) and (3.7). Differentiating both sides of these equations gives where the diagonal entry represents the global anomalous dimensions Γ J and Γ U , which can be calculated from Z J and Z U according to (B.2).

Resummed expression
Using the RG equations we can evolve each function from its characteristic scale where there are no large logarithms, and we get the following resummed expression with U ({n}, µ t , µ j ) = P exp µ j µt d log µΓ({n}, µ) , where P denotes the path ordering in log µ. This evolution matrix generates additional collinear partons with m ≥ l, therefore we define {n } = {n 1 , . . . , n l } and {n} = {n 1 , . . . , n l , n l+1 , . . . , n m } to distinguish these two configurations. According to the momentum scalings in (2.2), one should choose the hard scale µ h , the soft and beam scale µ b , the jet scale µ j and the coft scale µ t with the following typical values so that there are no residual large logarithms in the corresponding functions at these scales 6 .

Anomalous Dimensions for NLL resummation
To perform NLL resummation, one needs to include tree-level hard, jet, beam, soft, and coft functions, and evolves them using two-loop cusp anomolous dimension and one-loop regular anomolous dimensions. In this section we will provide all the regular one-loop anomalous dimensions relevant for the NLL resummation. In the calculation we neglect the difference between p T and p J T (recall that p T ≡ ( p J T − p V T )/2).

Anomalous dimensions of hard, soft and beam functions
The one-loop hard anomalous dimension is given by [79] with γ H ij→V k ≡ 2γ i (α s ) + 2γ j (α s ) + 2γ k (α s ), (3.17) where γ cusp is the cusp anomalous dimension, and γ f is the anomalous dimension of the parton species f (see appendix B). In dimensional regularization, at one-loop only real emission diagrams contribute to the soft function. Using the covariant gauge one has with δ + (k 2 ) = δ(k 2 )θ(k 0 ) andμ 2 ≡ µ 2 e γ E 4π . The evaluation of the soft function boils down to the calculation of the following master integrals, (3.19) 6 As shown in section 4, the hard function can have additional logarithms because it depends on two scales p J T and mV .
With the regulator we use, ω 12 involves a scaleless integral and hence vanishes. The divergent parts of the other two integrals are given by where L ⊥ ≡ log and φ x represents the azimuthal angle between x T and n JT . From the above expressions we can identify the soft anomalous dimension.
For 1/x T ∼ q T Λ QCD , one can calculate the beam functions from PDFs by an operator-product expansion [1,70,71] Since we are only resumming large logarithms of ln(µ/µ b ) at NLL level in this paper, we will neglect the non-logarithmic terms in I at O(α s ) in the following sections and only need the anomalous dimensions of the beam functions.
Let us focus on the non-PDF anomalous dimensions. The divergent pieces of the bare beam functions using the rapidity regulator in (3.2) take the form Since the soft and beam functions have the same characteristic momentum scale ∼ 1/x T , one can evolve the product of these functions from µ b ∼ 1/x T to µ instead of running each one of them individually. We write where the function W ij→V k satisfies the following evolution equations, with the anomalous dimension at one-loop level (3.25)

Anomalous dimensions of jet and coft functions
The global coft anomalous dimension Γ U can be derived from the one-loop calculation of the coft function U 1 . Explicitly, U 1 contains two Wilson lines, one along the n J direction and the other one along then J direction. After expanding Wilson lines in (2.23), at one-loop we have We only need the divergent terms in order to obtain the anomalous dimension, and we find which gives, From the definition (3.11), the global jet renormalization constant at one-loop is written as 29) and the anomalous dimension Γ J k has the form At one-loop level, it is the same as the one for the unmeasured jet function defined in [80]. In our framework, Z J is given by Z J 12⊗ 1 at this order, where Z J 12 removes the divergence in the jet function J 2 . However, beyond one-loop level a simple correspondence between Z J and the renormalization constants of the unmeasured jet function does not exist [34].
Finally, we will discuss the NGL resummation. At the NLL level, the non-global evolution matrix from the coft scale to the jet scale reduces to where we truncate the first sum in (3.14) at the tree-level jet function J 1 = 4πδ (d−2) ( n J 1 ⊥ )1, and we only include the tree-level coft function U m = 1. The non-global anomalous dimension defined in (3.12) has the following form [34] where we only need the one-loop results for the NLL resummation. The matrix elements are given by where T i,L are the color generators acting on the i-th particle in the amplitude and T i,R are the ones acting on the conjugate amplitude. The angular dipole factor W k ij is defined as The second line in (3.33) corresponds to the global anomalous dimension subtracted out from V m , where we define n 0 =n J . By expanding U k NG (µ t , µ j ) as a series of we have the evolution factor U NG (µ t , µ j ) = t U 1 + · · · with one-and two-loop coefficients as From this, one can show that the coefficient of the leading NGL at two loops is −4C k C A π 2 /3, which is the same as the results in [20]. As shown in [34,81], the coft function maps onto the hemisphere soft function under a Lorentz boost along the jet axis. Therefore, the evolution of the function U NG should be the same as Dasgupta and Salam's parametrization in [20]. Explicitly, in our numerical calculations we have Here u = 2t = 1 β 0 log αs(µt) αs(µ j ) , and the constants are given as a = 0.85C A , b = 0.86C A and c = 1.33.

NLL resummed expression
After plugging in with the above expressions, the all-order resummed expression (3.13) could be reduced to Let us compare our NLL resummed expression with those in [24,25,28]. The large logarithms log(Q/q T ) were resummed in γ+jet [24] and Z+jet [25] events at NLL level using the CSS formalism [1]. In these references, the calculations were carried out in the small R limit, but the terms with log R in the coefficients are not completely resummed. NGLs were also neglected. In the effective field theory language, this simply means that one does not distinguish the n J -collinear mode from the hard mode, and also not distinguishing the coft mode from the soft mode. By taking µ j = µ h , µ t = µ b and switching off the NGL resummation, our resummed expression reduces to those used in [24,25]. This can be shown more explicitly if one takes µ = µ b and On the other hand, in [28] the authors performed a resummation of log R and log(Q/q T ) without resumming non-global logarithms. If we take U k NG = 1, our resummed expression formally reduces to their results.

Analysis of Leading Logarithms
In this section we analyze the LL resummation. We shall study two cases of p T m V and p T m V , respectively. The first case is relevant in the studies of γ+jet production since photon is massless, or massive boson +jet production at high p T , while the second case is relevant in massive boson+jet production at low p T . Since leading logarithms are insensitive to scale choice, we choose µ b = b 0 /x T , µ t = R b 0 /x T and µ j = R p T in the following discussions.

Leading Logarithms for p T m V
In this case all the collinear particles typically carry an energy of order p T . Therefore one can simply make the following replacement Then at LL level, (3.13) reduces to the following form where dn(q T )/dq T is the differential probability of the boson+jet transverse momentum q T . We have only kept the LL terms in performing the Fourier transformation by using the relation We find that the resummation formula used in [24,25] give the same result at LL level. The double logarithms in (4.2) arise from soft radiation along the three collinear directions. Using a physical gauge, such as the light-cone gauge, the soft gluon spectrum is given by Let us first calculate the beam contributions to the integral distribution n(q T ). The phase space constraint for the soft gluon is as follows, k T q T , k T ω p T for real emissions, k T ω p T for virtual contributions. (4.5) The real and virtual cancellation yields the phase space shown as the shaded region in figure 2 (a) with Q = p T , and this gives On the other hand, soft radiation along the jet direction results in the log R-dependent terms. For simplicity, we assume that the jet is central. If a gluon is emitted inside the jet, there is no additional constraint on its phase space since it does not change the value of q T . If the gluon is emitted outside the jet, its energy has to satisfy ω q T . Combined with the virtual contribution, one can see that the phase space of the gluon is given by figure 2 (b) Figure 3. Illustration of the log R dependence in dn(q T )/dq T for the qg channel at LL. Here we take α s (q T ) ≈ 0.18 with q T ≈ 10 GeV around the peak region.
By including uncorrelated multiple soft gluon radiation, one obtains the Sudakov factor of the form It is easy to see that differentiating n(q T ) with respect to q T gives (4.2), which brings down from the exponent a factor given by single-gluon emission along the three collinear directions. Note that n(q T = p T ) = 1 which recovers the whole probability. Also, dn(q T )/dq T peaks at withᾱ = 4(C i + C j ) αs π . The peak location moves to a larger value of q T for smaller R because the probability becomes larger for a gluon to be emitted outside the jet. Also, the height of the peak becomes lower. All the above-mentioned features are illustrated in the left plot of figure 3, which shows dn(q T )/dq T with R = 0.4 and 0.8, and we set α s = 0.18 7 .

Leading Logarithms for p T m V
In this case the collinear radiation along the two beam directions typically has an energy ∼ m V , while the energy of the collinear particles inside the jet is of order ∼ p T . As a result, the phase space for soft radiation collinear to the jet direction is unmodified as in the previous case in figure 2 (b). On the other hand, m V , instead of p T , sets the phase space of soft radiation along the two beam directions, which is given by figure 2 (a) with Q = m V . Based on this physical argument, one expects the following logarithms to show up in the calculation: − αs π (C i + C j ) log 2 (m V /q T ) and 2αs π C k log(p T /q T ) log R. At LL accuracy, using (3.13) one can simply set Plugging them into (3.13), we get Note the additional, q T -independent logarithms log 2 (m V /p T ) appearing in the resummed result, which gives an overall normalization constant. In contrast, the LL result in [25] is given by (4.2) with p T replaced by m V , which is different from our result (4.12). Such differences can be seen in the right plot of figure 3, where we use the legend of "partial log R" to distinguish these two cases. Note that in the p T m V limit, the large logarithms of log(m V /p T ) need to be properly resummed which requires a factorization of the hard sector at the two scales m V and p T . We leave the study of constructing such an effective theory for future work.

NLL Resummation and Phenomenology
In this section, we study the Z+jet production in proton-proton collisions at √ s = 13 TeV in the high p J T case (p J T > 200 GeV) and the low p J T case (p J T > 30 GeV). We impose the constraint |η J | < 2.4 on the jet pseudo-rapidity and allow all values of boson rapidity. We then compare our theoretical predictions at NLL accuracy with Pythia simulations (version 8.2) [82] and the CMS data [56,59].

Characteristic Scales and Numerical Evaluations
We choose the following characteristic scales, Note that both µ b and µ t depend on x T , and one needs to include nonperturbative contributions when these scales approach Λ QCD . We focus on the effects of resummation in  Figure 4. Effects of log R resummation illustrated in the high p J T (left plot) and low p J T (right plot) cases. Here, NLL p stands for the NLL resummed results excluding NGLs. The label "partial log R" corresponds to the results from setting µ j = µ h and µ t = µ b .
perturbative QCD, and we simply impose an upper limit of x T < x max T = 1.5 GeV −1 in the x T -integral [11]. By varying x max T from 1 GeV −1 to 3 GeV −1 , we find that the dependence of the q T distribution on x max T is negligible compared to the uncertainties from scale variation. On the other hand, in the large q T p T region where µ t > µ j and µ b > min(µ h , p T ), the effective theory is no longer valid and we set µ t = µ j and µ b = min(µ h , p T ). In this region we need to switch off resummation and match the resummed results with the fixedorder predictions. However, different matching schemes will introduce additional source of uncertainties. We focus on estimating the theoretical uncertainty from scale variation since it is the dominant uncertainty at NLL accuracy. We leave the detailed studies of fixed-order matching and next-to-next-to-leading logarithmic (NNLL) resummation for future work.
The differential q T distribution dσ/dq T is calculated by numerically integrating over all the variables in (3.13). For the NLL resummation, we need all the one-loop anomalous dimensions in section 3.2 and the two-loop cusp anomalous dimension in appendix B. The exponential factors in (3.13) are evaluated analytically according to (B.10) and (B.11). We take the beam functions to be equal to the CT14 NLO PDF set [83] at the scale µ b according to (3.40). There is a constraint coming from requiring the φ x -integral to be convergent. Recall that both the soft and coft anomalous dimensions depend on cos φ x . The φ x -dependent terms can be combined and factored out as, αs(µ t ) . (5. 2) The φ x -integral is convergent only if  Figure 5. The effect of NGL resummation illustrated in the high p J T (left plot) and low p J T (right plot) cases. Here, NLL stands for the full results calculated from (3.13), and NLL p corresponds to the results without NGL resummation.
One encounters such a divergence when the coft scale approaches to the non-perturbative region. It would be intriguing to see how one can introduce nonperturbative functions to tame such a divergence. We instead only integrate x T over the region given by (5.3).

Effects of log R and Non-Global Logarithm Resummation
We study the effects of log R and NGL resummation at NLL accuracy.

log R resummation
Here we switch off the contribution from NGLs by setting U k NG = 1 in (3.13). We define the resummation accuracy without NGLs as NLL p where the subscript p means partial. Furthermore, we compare the NLL p result with the one by setting µ j = µ h , µ t = µ b which is denoted by "partial log R" resummation since part of the log R dependence is eliminated in the scale ratios. Note that, in the high p T case µ j = µ h ∼ p T while in the low p T case µ j = µ h ∼ m Z , and that the characteristic scale µ j = p T R. Therefore in the low p T case the "partial log R" results differ from the NLL p result by the missing contributions of the form log(m Z /(p T R)) as discussed in section 4.2. Figure 4 shows the effect of log R resummation in the high p T (left plot) and low p T (right plot) cases. The NLL p cross section is always larger than that with partial log R resummation. Note the significant effect on the overall cross section especially in the low p T case. As discussed in section 4, one can see that the overall factor e αs π C k log m Z p T log m Z p T +2 log 1 R (5.4) accounts for the cross section difference, which clearly comes from the running of the jet function between m Z and p T R.

NGL resummation
As discussed in section 3, NGLs arise from one coft gluon radiated outside the jet. The contribution at O(α 2 s ) takes the form − α 2 5) and the contribution increases as the ratio p T /q T increases. Therefore NGLs are expected to play a more important role at high p T . Figure 5 shows the cross sections calculated from (3.13) with (denoted by NLL) or without (denoted by NLL p ) the NGL resummation.
One can see that the NGL resummation lowers the peak of the cross section and pushes it to a larger value of q T .

Theoretical Predictions and Uncertainties
We compare our theoretical predictions with Pythia simulations and experimental data at the LHC. We estimate the theoretical uncertainties by varying each characteristic scale in (5.1) by a factor two and taking the envelope of all the results from scale variations. For both the p T > 30 GeV and p T > 200 GeV cases, we calculate dσ(q T )/dq T with R = 0.4, 0.6 and 0.8. As shown in figure 6, our theoretical predictions agree reasonably well with the Pythia partonic results within the uncertainty band 8 . However, some discrepancy in the overall cross section exists, especially for R = 0.4 with a smaller jet radius.
We then compare our theoretical calculation with experimental data. In order to impose the same cuts on kinematic variables as the experiments, we use the LO hard function including the leptonic decay of Z/γ * . We first compare with the data at √ s = 13 TeV in [59]. We impose the same kinematic cuts as The left plot of figure 7 shows the comparison between our prediction for dσ(q T )/dq T with the data 9 . Our result is consistent with the experimental data in the small q T region. We also show the result of the full LO distribution (black curve) calculated by MCFM program [84,85], and the one including only the logarithmic terms at LO (orange curve) predicted using SCET. In appendix C we give the expressions of LO singular terms. In the small q T region fixed-order expansion breaks down because of large logarithms of log(Q/q T ), and SCET can reproduce this singular behavior. For the large q T region, we need to include power corrections from fixed-order calculations. However, near q T ∼ 30 GeV where the p J T > p min T = 30 GeV selection is imposed, 8 We checked that the major difference between the partonic and hadronic results comes from multiparton interaction contributions. 9 Note that in experiment qT is defined as the sum of the transverse momenta of the Z boson and all the jets with p J T > 30 GeV and |ηJ | < 2.4 in the event [59], while in our calculation we only include the leading jet in defining the qT . From Pythia simulations, we find that using the leading jet to define qT brings down the first three bins of dσ/dqT (left plot of figure 7) Figure 6. Comparison between the NLL cross section calculations with Pythia simulations, in the high p J T case (top row) and the low p J T case (bottom row). In all the plots, the red curves are the theoretical predictions with the scale choice in (5.1), and the error bands are shown as the shaded regions. The histograms are the Pythia results at parton (dashed lines) and hadron (solid lines) levels.
the LO result has an artificial kink structure. The kink structure comes from the negligence of two jet events with p J T < 30 GeV due to such a kinematic cut. Explicitly, at LO p T and q T are the transverse momenta of leading and subleading jets, respectively. When q T > 30 GeV, the lower limit of the p T integral is q T . On the other hand, for q T < 30 GeV the lower limit is frozen at 30 GeV. Hence, we observe such kink structure near q T ∼ 30 GeV. The investigation of the kink and its treatment is beyond the scope of this paper and left for future work.
We also compare our theoretical calculation of the azimuthal angle decorrelation ∆φ between the boson and the leading jet with the experimental result at √ s = 7 TeV in [56]. In the numerical integration, we boost the tree-level partonic event such that the boson and the leading jet have total transverse momentum q T as q T = q T (sin φ q , cos φ q ). (5.7) After performing this transformation, the Z boson and the leading jet are not back to back in the transverse plane. Hence, we obtain the distribution of the azimuthal angle ∆φ(Z, j 1 ) In principle, one needs to perform a matching between the resummed result and the fixedorder calculation to calculate the azimuthal angle decorrelation (see, e.g., [22,24]). However, as we show in figure 6, at high p T our resummed result gives a good description even up to q T ∼ p T . We then use it to calculate the normalized distribution dσ/d∆φ by integrating out q T from 0 to 150 GeV. We find reasonable agreement with the experimental result.

Summary and Perspective
In this paper, we construct an all-order formalism in SCET for the systematic resummation of large logarithms of the form log(Q/q T ) when q T Q in boson+jet production in the small R limit. More precisely, the expression (3.13) resums the logarithms of log(Q/q T ), log R and non-global logarithms, at NLL accuracy. We first carried out an analysis of the leading logarithms. We find that in the the case of p T m V , the resummation of the logarithms log 2 (m V /p T ) is missing in the literature. In the case of p T m V the effect of log R resummation only comes in at NLL accuracy. At the end, we compare our theoretical predictions with Pythia simulations and available experimental data [56,59]. Within theoretical uncertainties, our results are consistent with the simulations and the data.
In the present work we obtained the resummed cross section at NLL accuracy. There are several issues that we leave for future studies. First, as shown in the plots at NLL accuracy, there are relatively large uncertainties both at small-q T and large-q T ∼ p T . The small-q T region shows the sensitivity to non-perturbative physics. In this case one needs to introduce non-perturbative functions to extend the x T -integral into the non-perturbative regime, and to regularize the singularity in the integration over φ x . On the other hand, the large uncertainties at large-q T is a signature of the breaking-down of the scale separation in (2.2). The improvement in this region is usually obtained only after matching the resummed result with a fixed-order calculation at higher-order in α s . In detailed phenomenological studies, one needs to include all these improvements and possibly perform the resummation at NNLL accuracy in order to reduce the overall theoretical uncertainties. Second, the possible breaking of the transverse momentum factorization in the processes studied in this paper is another intriguing issue. On the other hand, in Pythia simulations we see only small non-perturbative corrections. This suggests that the q T distribution in boson+jet production can be a clean and useful probe of factorization violation and Glauber contributions. Last, but not least, the observables studied in this paper can be used to measure the jet-quenching parameters in high-energy nuclear collisions [24,46,47]. The formalism presented in this paper can be used as a unified formalism to study boson-jet correlation in both proton-proton and high-energy nuclear collisions.
From the conservation of the + and − components of the momenta in the n 1 andn 1 basis, one has The amplitudes squared, averaged and summed over the color and spin indices in initial and final states are given by |M (qq → V g)| 2 = 16π 2 α s α em e 2 q (N 2 c − 1) N 2 where e q is the electric charge of the quarks in the case of photon production. For Z production we need to replace e q by e 2 q → 1 − 2 |e q | sin 2 θ W 2 + 4e 2 q sin 4 θ W 8 sin 2 θ W cos 2 θ W (A. 5) with θ W the weak mixing angle. In our numerical calculation the electroweak parameters we adopted are All the anomalous dimensions exceptΓ used in (3.13) consist of a cusp part, which gives the leading logarithms, and a non-cusp part, which only contributes to sub-leading logarithms. That is, these anomalous dimensions take the form Γ(α s ) = C Γ γ cusp (α s ) ln with r = α s (µ)/α s (ν).

C LO singular terms
In this appendix we give the analytical expressions of the LO singular terms. After expanding the resummed result (3.39) order by order in α s and performing the Fourier transform, we can obtain the singular terms of the q T distribution. The LO results are given by ij←ab (z 1 , z 2 , q T ) , (C.1) where the one-loop kernel Σ (1) ij←ab has the following form with the coefficients 3) The one-loop Alatrelli-Parisi splitting functions are given as follows,