Precision boson-jet azimuthal decorrelation at hadron colliders

The azimuthal angular decorrelation of a vector boson and jet is sensitive to QCD radiation, and can be used to probe the quark-gluon plasma in heavy-ion collisions. By using a recoil-free jet definition, the sensitivity to contamination from soft radiation on the measurement is reduced, and the complication of non-global logarithms is eliminated from our theoretical calculation. Specifically we will consider the $p_T^n$ recombination scheme, as well as the $n\to \infty$ limit, known as the winner-take-all scheme. These jet definitions also significantly simplify the calculation for a track-based measurement, which is preferred due to its superior angular resolution. We present a detailed discussion of the factorization in Soft-Collinear Effective Theory, revealing why the transverse momentum $\vec q_T$ is more complicated than the azimuthal angle. We show that potential glauber contributions do not spoil our factorization formalism, at least up to and including order $\alpha_s^3$. The resummation is carried out using the renormalization group, and all necessary ingredients are collected or calculated. We conclude with a detailed phenomenological study, finding an enhanced matching correction for high jet $p_T$ due to the electroweak collinear enhancement of a boson emission off di-jets. We also compare with the Pythia event generator, finding that our observable is very robust to hadronization and the underlying event.


JHEP02(2023)256
Our main result consists of a phenomenological study, consisting of numerical results for the analytic resummation at NNLL accuracy. We include a discussion of the resummation, perturbative uncertainties, nonperturbative effects, and the matching to fixed order. We furthermore study features of this observable using Pythia, and compare to our NNLL result. As promised by the recoil-free axis, we find that the observable is insensitive to soft effects (hadronization, underlying event), a rather small sensitivity to the jet radius, and negligible differences when measuring the decorrelation on tracks. The uncertainty band of our resummed predictions are reduced when going from NLL to NNLL. Our resummed results are consistent with Pythia, except for the matching corrections from the NLO cross section (as we didn't match Pythia to NLO). These matching corrections are substantial for high jet p T , and arise from the boson being emitted from a leading order dijet configuration. Though this is formally power suppressed in the back-to-back limit, it is enhanced by an electroweak logarithm. While we focus on the WTA axis, we also explore p n T -weighted recombination schemes, and these conclusions also hold there (if n > 1). The paper is structured as follows: We introduce the kinematic setup and discuss different observables measuring the transverse momentum decorrelation of the V +jet pair in section 2. Section 3 establishes the factorization formula, gathers the available ingredient functions, and includes a brief discussion of factorization violation induced by a Glauber mode. Section 4 contains various calculations relating to the jet function: The linearly polarized jet function appearing in our factorization formula, as well as the changes to the jet function induced by the use of a different recoil-free jet axis or track-based measurements, respectively. We establish our resummation strategy in section 5 (with certain ingredients relegated to appendix A) which is then carried out to derive the results in section 6. We conclude in section 7.

Geometry of the collision
We begin by describing the geometry of the collision and defining the observable for which we perform the resummation. It will be instructive to contrast our target observable, the azimuthal decorrelation, with the closely related radial decorrelation. (These correspond to the two components of the difference in transverse momentum between vector boson and jet.) This comparison will demonstrate in detail where simplifications due to the choice of a recoilfree recombination scheme arise and how non-global effects are suppressed: While using a recoil-free axis removes non-global logarithms for the azimuthal decorrelation, the radial variety still suffers from NGLs and needs to include effects related to the technical definition of this axis. For concreteness we discuss the case of the Winner-Take-All (WTA) axis.
The geometric setup is illustrated in the left panel of figure 1, where we choose to align the y-axis with the reconstructed jet axis (its projection onto the transverse plane, to be precise). Our starting point is momentum conservation in the transverse plane, which reads p T,a + p T,b + p T,S + p T,c + p T,V = 0, The azimuthal angle between the vector boson (green) and jet axis (blue) is related to the momentum of the vector boson p x,V transverse to the colliding protons (red) and jet. Collinear initial (purple) and final-state (blue) radiation and soft radiation (magenta) is also shown. Right: Schematic of the transverse plane: the angle φ c between WTA axis (along y axis) and collective collinear momentum, the angle φ ij between two generic collinear emissions i and j, and the angular decorrelation observable δφ. These quantities are of the same parametric size. and expresses that the vector boson (with transverse momentum p T,V ) recoils against emissions off the beams ( p T,a , p T,b ), soft radiation ( p T,S ), and the total collinear radiation in the jet ( p T,c ). Note that p T,c is non-trivial, i.e. possesses both x and y-components, as we chose to align our coordinate axes with the reconstructed jet axis p T,J , which in the WTA case does not follow the collective momentum of all collinear emissions in the jet: As detailed in ref. [34], the WTA axis always follows the orientation of one input particle (not necessarily the most energetic one, which would not be collinear safe), whereas the total momentum in general is not aligned with any individual particle direction.
To illustrate this point, we revisit the WTA recombination scheme and the associated axis in detail: in sequential jet clustering algorithms, particles are pairwise recombined and assigned a joint 4-momentum. For the WTA recombination scheme, this joint momentum is taken massless, and points along the direction of the "harder" of the two particles in the pair. The criterion to determine which of two particles is the harder is typically the hierarchy in either the particle energy ("WTA-E-scheme") or the transverse momentum ("WTA-p T -scheme"). 2 The latter choice is of course more suitable for hadron colliders and shall therefore be used here. The magnitude of the newly combined pseudo-particle is the sum of the energies or the scalar sum of transverse momenta -the "winner" absorbs the JHEP02(2023)256 loser, so to speak. The WTA axis is then simply the final pseudo-particle after clustering. It is in principle possible that many soft emissions are clustered into one highly energetic pseudo particle that emerges as the winner, but this is extremely unlikely because they would have to be clustered together first to stand a chance against collinear emissions, contrary to the expectation that they are spread out over the jet. Thus the WTA axis will be aligned with one of the collinear emissions, rendering its direction free of soft recoil.
The WTA axis is a particularly powerful tool when combined with SCET, because it can exploit the parametric hierarchy of the modes in SCET to simplify the calculation: As stated, the magnitude of the jet axis will be a sum of transverse momenta where the sum runs over all emissions in the jet. One consequence of the WTA recombination is then that soft emissions inside the jet can affect the magnitude as participants in this sum, but never the orientation, as soft emissions always lose against collinear emissions when determining with which emission the axis should be aligned. The dependence of the axis on soft emissions can be expanded away and represents a subleading effect in the effective theory -and in some cases even the dependence in the magnitude can be expanded away, with more details laid out in the following two sections. The azimuthal angle is now directly related to the (dimensionful) offset between vector boson and reconstructed jet momentum, defined as In the absence of QCD radiation, q T vanishes, and the limit of small q T is of interest for the resummation. The angular decorrelation can be e.g. written as δφ = arcsin q x /|p T,V | ≈ q x /|p T,V |, and is essentially a dimensionless version of the dimensionful q T .

Azimuthal and radial decorrelation
Besides the azimuthal decorrelation, which represents the tangential offset q x as shown in figure 1, we can define a second quantity of interest here, the radial decorrelation, using q y . Using eq. (2.1) these can be written as where we used that p x,J = 0 due to the alignment of the y-axis with the WTA axis. We expect the contributions from the beams to be isotropic in the transverse plane (and their xor y-dependence therefore to be similar) and the soft function to capture the y-dependence via its dependence on the geometry of the beam-jet system, so to proceed we need to understand the relation between p T,J and p T,c , i.e. the reconstructed WTA axis and the -5 -

JHEP02(2023)256
collective collinear momentum of the jet, and their components. 3 Right from the start it is clear that the radial decorrelation is more complex, as it involves a large cancellation between the boson and jet axis magnitude. For p T,c and its decomposition into x-and y-components, this yields the situation illustrated in the right panel in figure 1. Thus where p T,c (which is of the order of the center-of-mass energy Q) is the vector sum of all collinear emissions off the parton initiating the jet, and φ c measures the angular separation between the WTA axis and the collective collinear momentum flow.
To obtain δφ or q x , we still need to include the collinear and soft contributions in eq. (2.4). As the x-direction is perpendicular to the plane containing both the jet and the beams, we expect the appropriate SCET modes to scale -with notation p µ = (n b · p, n a · p, p ⊥ ) -as where the scaling for the jet -with subscript J -is to be understood as scaling of lightcone components based on the jet direction, i.e. p µ = (n J · p, n J · p, p ⊥ ) J . We thus find a SCET II situation, as expected for a de facto transverse momentum measurement. Finally, note that there is no special dependence on soft in-jet radiation: It is always expanded away when the jet axis direction is determined, but it can also be expanded away in the magnitude in this case, as the azimuthal decorrelation only depends on the direction of the jet axis, which determines the coordinate system we use to define it.

Properties of the radial decorrelation
We will now contrast this surprisingly simple observable with the radial decorrelation q y , where the details of the recombination become explicitly relevant. We emphasize that we do not wish to actually perform the resummation for this observable, but merely include this discussion to highlight why the azimuthal decorrelation is so amenable to precision calculations.
A key insight is that while p y,J in eq. (2.4) includes a sum over collinear (and soft) momenta, it makes use of a scalar sum and aligns it with a collinear emission's direction, while p y,c is a component of a vector sum. Iterating the law of cosines/binomial theorem over all relevant momenta, the vector and scalar sum of collinear momenta can be related, to Here φ ij are the pairwise opening angles between emissions i and j, which satisfy φ ij 1 because i and j are both collinear. Thus we can expand the above expression to obtain 3 This is also true for qx, which depends implicitly on pT,J , through our choice of coordinate system.

JHEP02(2023)256
We can now derive the expression for the radial decorrelation: As p x,J = 0 due to our choice of conventions, p y,J is the scalar sum of all momenta inside the jet, which for R δφ includes all collinear radiation, as well as the soft emissions that end up in the jet: where the leading scalar sum of collinear momenta cancelled, and the corresponding functions in the factorization are indicated. The expected power counting of the modes involved is also non-trivial: The jet power counting is determined by its ⊥-component, which must scale as O(Qφ c ), as the angle between WTA axis (i.e. one of the emissions) and the collective collinear momentum has to be parametrically similar to the opening angle of collinear splittings. The transverse components p y,a and p y,b for the beams must scale as O(Qφ 2 c ) to contribute at equal footing with the jet contributions in (2.9). The soft contributions are isotropic and a similar argument applies. It thus follows that which is an interesting combination of SCET I and SCET II characteristics. For the beams this is a transverse momentum measurement, and so they share the virtuality of the soft emissions, as typical in SCET II . The jet, however, measures essentially an invariant-mass-like observable, a SCET I situation.
We also note that the differing treatment of soft radiation inside and outside the jet means this observable exhibits non-global logarithms. This is a persistent feature: while it is e.g. possible to remove the term depending on pairs of collinear emissions from eq. (2.9) by determining the magnitude of p T,J using a vector-sum-based recombination scheme, the inclusion of soft in-jet emissions will also be present in such a case.
From the above discussion we can now understand why the azimuthal decorrelation is so simple, by comparison: The non-trivial effects present in its radial sibling, including the appearance of non-global logarithms, are O(φ 2 c ) effects, and are therefore power suppressed compared to the dominant contributions to azimuthal decorrelation in (2.5), which are O(φ c ).

Extensions
We conclude this section, by discussing other recoil-free axes, double differential measurements in the azimuthal and radial decorrelation and tracks for the radial decorrelation: First, we note that the choice of a recoil-free axis other than the WTA-axis (in essence the FastJet [41] recombination scheme of equation (4) with w i = p n T i and 1 < n < ∞, see also eq. (4.20)) constitutes a small change for the azimuthal recombination: as long as it is recoil free, only the jet function will be changed, and we will calculate the corresponding  (For the radial decorrelation this may be a more intricate affair,  as subleading effects are elevated by virtue of a large cancellation.) Second, we point out that measuring any quantity in the q-plane other than the elementary q x and q y is a de facto double-differential measurement, and requires the specification of the relative scaling of the q x and q y components, to cleanly establish which contributions to the observable from the different regions (or ingredient functions) are subleading and can be neglected. This ultimately is a consequence of the different scaling of q x and q y themselves, in terms of power counting. In particular, one could consider q x ∼ q y or q 2 x ∼ Qq y , which correspond to making the factorization of the azimuthal or radial decorrelation differential in both q x and q y , or something in between (see e.g. ref. [42] for a factorization description of double differential measurements).
Finally, we point out that using track-based measurements is a possibility for the azimuthal decorrelation, and will be discussed in 4.3, but is firmly ruled out for the radial decorrelation. The reason is simply that the vector boson will always be fully reconstructed with all of its transverse momentum, but the jet axis is only established based on the particles included in the jet recombination. If all collinear particles are clustered into the jet, the jet axis magnitude is roughly equal to the total collinear transverse momentum, and the cancellation against the vector boson transverse momentum establishes the radial decorrelation as a small quantity, which vanishes in the singular limit. A resummation program can then be set up. If tracks are used, however, only the subset of charged particles in the final state contributes to the jet axis. The difference to the vector boson p T is therefore a large quantity, finite even in the singular limit, and given by the transverse momentum of the neutral particles in the jet. The azimuthal decorrelation is unaffected, as it does not involve a large cancellation.

Factorization
In the previous section we performed an analysis of the contributions to the transverse momentum q of the boson+jet system from soft and collinear emissions, indicating that these contributions are dominant at leading power. In the limit where q x is small, the infrared structure of QCD results in large logarithms of ln(Q/q x ) caused by soft and collinear emissions from initial state beam partons and final state jet partons (and similarly for q y ). To obtain reliable predictions in this limit, we will derive a factorization formula in section 3.1, that enables us to resum these large logarithms, which will be discussed in section 5. Resummation to next-to-next-to-leading logarithmic accuracy requires the ingredients of the factorization theorem at one-loop order, given in section 3.3 or calculated in section 4, 4 as well as their anomalous dimensions at two-loop order (and the three-loop 4 The calculation of the one-loop gluon jet function is detailed in section 4.1, showing for the first time how to obtain the linearly-polarized contribution. Generalizing the WTA recombination scheme to the recoil-free p n T scheme (with n > 1) only results in a modification of the jet function, which is calculated in section 4.2. Switching to track-based measurements also only affects the jet function constant, this calculation is outlined in 4.3.

Factorization formula for azimuthal decorrelation
As the sum of momenta in (2.4) shows, the component q x of the vector boson perpendicular to the jet axis receives contributions from a linear sum of four terms. Three of these terms represent the characteristic emissions from one of the hadronic directions, and by choosing the WTA axis, the other term is the momentum component of all soft radiation in the collision.
We will consider the case where the azimuthal decorrelation δφ = arcsin(q x /p T,V ) ≈ q x /p T,V R, so that contributions from out-of-jet emissions are suppressed by powers of δφ/R. This implies that the observable δφ is not sensitive to the non-global correlation between radiation inside and outside the jet. Starting with these assumptions, we can factorize the cross section, which allows us to calculate the δφ distribution at high logarithmic accuracy. The corrections to this factorization are suppressed by powers of δφ, and can be included for the region where δφ is not small by matching to the fixed-order prediction for the cross section, as discussed in section 5.2.
A factorization formula for the boson-jet transverse momentum imbalance q T was derived by some of us in ref. [25], for the case where jets are defined using the anti-k t algorithm with the standard recombination scheme. In that case, the jet axis is along the direction of the total jet momentum, and is sensitive to recoil from soft radiation enclosed within the jet boundary. As a consequence, the cross section receives nonglobal contributions that involve soft radiation inside and outside the jet, preventing a simple factorization of collinear and soft contributions (though there has been substantial progress in resumming nonglobal logarithms, see e.g. [29,[43][44][45][46][47][48][49][50][51][52]). Here we instead use the WTA recombination scheme for jet clustering, which is only sensitive to the total amount soft emissions. The WTA axis is sensitive to the distribution of collinear radiation in the jet, but the contributions from collinear emissions off the beams are the same as for the standard jet axis case.
We start from a general differential momentum distribution of all the particles in an event and project it onto the observable q x , is the complete phase space of final state particles except the boson, whose transverse momentum and rapidity are denoted by p T,V and y V , respectively. 5 Similarly, the transverse momentum and pseudo-rapidity of the jet are labelled p T,J and η J . The delta function restricts the value of the observable q x , through the measurement functionq x , which is a function of final state particle momenta {p i }. The leading power 5 As the jet and vector boson recoil in the transverse plane, we could choose to keep pT,J differential, instead of pT,V . However, pT,J is the reconstructed jet momentum, and as such loses its sensitivity to the hard scattering kinematics if we measure the angular decorrelation using tracks: the neutral particles' contribution will then be missing, and pT,J will no longer be parametrically similar to pT,V . We therefore pick the more robust pT,V .
where dX a , dX b , dX c and dX S correspond to the differential phase space of collinear emissions of beam a, beam b and jet c, as well as soft emissions from these collinear directions. From (2.4) we see that the measurementq x simplifies, as with the contribution from each factorized IR sector given by 6 Therefore, at leading power, the differential cross section can be organized as follows, dσ dq x dp T,V dy V dη J = dp x,a dp x,b dp x,c dp x,S δ q x + p x,a + p x,b + p x,c + p x,S × dσ dp x,a dp x,b dp x,c dp x,S dp T,V dy V dη J . (3.5) In the next section we will discuss how the soft-collinear decoupling in the SCET Lagrangian allows us to factorize the multi-differential cross section dσ/( i dp x,i dp T,V dy V dη J ) to obtain dσ dq x dp T,V dy V dη J = dp x,a dp x,b dp x,c dp x,S δ q x + p x,a + p x,b + p x,c + p x,S (3.6) Here the indices i, j, k label the partonic channels of the hard scattering processes, encoded in the hard function H, producing a high-transverse momentum boson V recoiling against the jet. The contribution to q x of collinear initial and final-state radiation is described by the beam functionsB and jet functionJ , and the soft functionS accounts for the contribution from soft radiation. (The tilde signifies that these functions are defined in momentum space, not impact parameter space.) The treatment of Lorentz and color indices in deriving eq. (3.6) will be discussed in the next section, and leads to a linearly-polarized contribution from gluon beams and jets, for which there is a corresponding change to the hard function H. The Bjorken variables x a and x b are determined by the boson and jet kinematics, 6 Remember that px,c sensitively depends on the WTA scheme, since the y-axis is along the jet direction.

JHEP02(2023)256
We can eliminate the convolution in the above factorization formula by switching to the impact parameter variable b x , which is the Fourier conjugate of q x , The factorization formula in impact parameter space thus has a product form. In the next section we will present some details of the derivation of the above factorization formula, and discuss each factorized component.

Derivation of factorization formula
We will derive the factorization theorem using Soft-Collinear Effective Theory (SCET) [35][36][37][38][39]. For an introduction to SCET, see e.g. refs. [53,54]. The leading operators describing the hard scattering for Z+jet production are [55] O The superscripts S and T on the operators label different quark spin structures. The four-vector b denotes the position of the operator and the t i will be integrated over in the matching in eq. (3.10) below. The lightcone directions n a , n b and n c are along the two beams and jet directions, and χ n i = W † n i ξ n i (χ n i =ξ n i W n i ) is the collinear quark field and B µ n i ⊥ = 1 g W † n i iD µ n i ⊥ W n i the collinear gluon field, which include collinear Wilson lines W n i for gauge invariance. They appear in combination with soft Wilson lines S n i in the pattern outlined by the square brackets. Multipole expansions simplify the coordinate dependence of the soft Wilson lines, and are implied for the b-dependence of the collinear ones, as well. We have here labelled the fields according to the qq → Zg channel, and the other channels can be obtained by a permutation of the lightcone directions. To keep the notation compact in the following, we will denote these SCET operators as O ν j , with the index j labeling the partonic channels and quark spin structures, and denote the collinear fields collectively as [φ n i ] α i a i with Lorentz and color indices α i and a i . The contribution from operators with additional collinear fields is power suppressed by O(δφ).
The SCET operators are matched to the full electroweak current which produces the boson, where C j 's are the Wilson coefficients. In the case of direct photon production, the cross section differential in the infrared degrees of freedom has the following form after performing -11 -

JHEP02(2023)256
the integrals over the t i : where p µ is the polarization vector of the boson with the index p labeling the polarization states, and |P a and |P b are the initial hadron states with momenta P a and P b . Here α em is the electromagnetic coupling constant, e q is the electric charge of the quark (in terms of elementary charge units), and for Z production one replaces e 2 q by With the capability of tagging the boson polarization, one will be able to assign a specific polarization tensor p * µ q ν . For observables inclusive in boson polarization, the polarization sum and the Ward identity give Together with the factorization of infrared degrees of freedom, |X IR = |X a |X b |X c |X S , the above matrix element factorizes into a product of three collinear matrix elements and one soft matrix element. For particular quark spins or gluon polarizations either in the initial or final state, 7 one needs to project the fields onto the corresponding components. The collinear matrix elements associated with the incoming protons give the bare quark and gluon beam functions, Here [P α i α i n i ] j f is the Dirac or Lorentz structure of the beam function for parton flavor f and spin or polarization structure j, d i is the dimension of the color representation of a generic collinear field φ n i , and = (4 − d)/2 is the dimensional regulator.
For the quark beam function (i.e. i = a, b, f = q,q), only a single Dirac structure contributes at leading power, Another possible Dirac structure, γ µ ⊥ , would need to be combined with a transverse momentum, and is therefore power suppressed. For the gluon beam function (i = a, b, f = g),

JHEP02(2023)256
there are two independent Lorentz structures: transverse (T ) and linear (L) polarization, for which the projectors are [59,60], The linear polarization contributes, because the collinear gluon splitting is intrinsically polarized.
The collinear matrix elements associated with the outgoing jet are the quark and gluon jet functions. Before introducing any jet definition, one has One can, then, make the following decomposition Below, we shall drop the first phase factor on the right hand side of this equation, which specifies the jet momentum that enters the hard function and the conservation of p ± with respect to the beam directions, respectively. As discussed in section 2.1, the jet definition respects factorization, meaning that the jet functions are only dependent on n J -collinear modes. Effectively, a jet definition simply defines the transverse momentum and rapidity of the jet: Here the momentum p x,c corresponds to the x-component of the momenta of the jet constituents, and the operatorsˆ p T,J andŷ J give the jet momentum transverse to the beam and jet rapidity, respectively. The decomposition of the Dirac and Lorentz structures is 8 As shown in section 4.1, it is more convenient to evaluate jet functions in a frame in which px,J = 0.
Therefore, we keep the dependence of px,J explicit here.

JHEP02(2023)256
similar to the beam functions (with a different normalization), where this is now perpendicular to the jet (not beam) direction, as indicated by the ⊥ (instead of T ). The soft matrix element associated with the soft operator O s j give the soft function, where T (T) denote (anti-)time ordering. The overall factor is chosen such that at tree-level the soft function equals 1. The qqg color space is one-dimensional, but there is a different soft function for each partonic channel j, because it matters whether the jet or beam are in the adjoint representation for the gluon.
In the remainder of this paper we will label the soft function with the full partonic channel ijk, where i, j, k are now parton flavors.
The hard functions are the square of the matching coefficients of SCET operators and observable independent. These are at tree-level given by the corresponding QCD matrix-elements M and including factors that account for averaging over spin/color states of incoming partons (indicated by the bar in M ) and the flux factor, At one-loop order, this equation still holds after renormalization and dropping infrared divergences (which cancel in the matching between QCD and SCET). Note that the indices of the projectors in eqs. (3.16) and (3.21) will be contracted with the matrix elements of the hard scattering. If an observable is not sensitive to the spin or polarizations states, the typical spin-averaged hard function is sufficient. However, in our case we will get a different expression when the initial or final gluon is linearly polarized, and we need to sum all polarizations in eq. (3.8), which affects both the gluon beam/jet function and the hard function.

One-loop ingredients
Here we provide one-loop expressions of the beam, jet, soft and hard functions. Our observable leads to rapidity divergences in the beam, jet and soft function that are not regulated by dimensional regularization. We use the η-regulator [61,62] to regularize rapidity divergences, and the resulting evolution in the rapidity scale ν can be used to resum the corresponding rapidity logarithms. This is discussed in section 5.1, where also the anomalous dimensions needed for next-to-next-to-leading logarithmic resummation are collected.
The beam functions are matched to collinear parton distribution functions (PDFs) with perturbatively calculable matching coefficients [62][63][64][65] The matching coefficients I have been calculated up to three-loop order [66][67][68][69][70][71][72][73][74][75][76], and the linearly-polarized contribution at two-loop order [77]. Up to one-loop order, they are given by 26) and the splitting functions The one-loop jet functions for the WTA recombination scheme are [30,78,79] The result for the linearly-polarized gluon jet function was first quoted in our letter [30], and we therefore provide a calculation of the gluon jet functions in section 4.1. There we also calculate jet functions for other recoil-free axes and for track-based measurements. Up to order α 2 s , the soft function can be obtained [30] from the standard TMD soft function, which is known up to three loop order [70,80,81]. The contribution involving exchanges between three Wilson lines vanishes due to color conservation [82]. For exchanges involving only two Wilson lines, we can perform a boost to make them back-to-back. Since our observable is perpendicular to the boost, only the rapidity regulator is affected, which can be taken into account in a straightforward manner [83]. The resulting one-loop soft functions are where the color factor C i is C F if parton i is an (anti-)quark and C A if it is a gluon, and The hard function for the process ij → V k is given by for which the tree-level matrix elements are (3.32) and g L denotes a linearly-polarized gluon. For Z production one replaces e q in analogy to (3.12), and the partonic Mandelstam variables arê The loop corrections that enter the hard function have been calculated at one-loop order in refs. [84][85][86][87][88][89], and we give their expressions in the appendix B. The resulting one-loop hard function for a transversely polarized gluon can be obtained from the appendices of [55,90,91]. Since the beam and jet function for a linearly-polarized gluon only start at one-loop, the tree-level hard function suffices in this case.

Glauber interaction and factorization breaking
We conclude this section by briefly commenting on the appearance of a Glauber mode and its impact on the validity of the factorization theorem in (3.6) and (3.8), following the framework laid out in ref. [92], supplemented by the results of refs. [93,94].
Ref. [92] explains in detail that Glauber rungs connecting active lines to other active or to spectator lines can be absorbed 9 into the soft (for active-active) or collinear Wilson lines (for active-spectator), and only the pure spectator-spectator Glauber exchanges require special attention. Before we discuss our specific application, we first summarize and clarify a few points: First, we note that in line with the expectation that "spectator gluons or quarks may be created by collinear radiation from active lines", according to ref. [92], we take only the parton lines directly connected to the hard scattering vertex as "active". This implies that collinear splittings in the initial state, as well as the branching of the jet progenitor parton generate spectators, and that thus there are three collinear sectors containing spectators (two proton remnants, one jet sector).
Second, we note that for single-scale observables, refs. [93,94] demonstrate that pure Glauber exchanges (i.e. without soft emissions off Lipatov vertices) do not lead to factorization violation.
Third, we resolve the ambiguity of the Glauber mode for three distinct rapidity sectors by observing that Glauber modes only ever connect two of them, and that the frames in which these are pairwise back-to-back are connected by -generically -O(1) boosts. Together with the collapse rule this means that the different Glauber dipoles communicate only via soft emissions from Lipatov vertices, whose power counting is isotropic, and therefore unchanged by the boosts relating the different Glauber frames.
Lastly, we follow the diagrammatic conventions of ref. [92] for the discussion below, where auxiliary interactions are introduced that create active and spectator lines from incoming hadron fields, to be able to focus on the topology of the appearing diagrams.
Moving on to the discussion of our case: We begin with Glauber exchanges involving Lipatov vertices off Glauber rungs connecting to at least one active parton, as e.g. in 9 Alternatively, the Glauber mode is then a true subset of the corresponding soft or collinear mode, so a dedicated Glauber contribution does not add anything that is not already covered by the naive soft or collinear mode.  [92], the Glauber exchange is an element of the soft or collinear sector, the Lipatov vertex therefore represents soft radiation off either another soft, or off a collinear emission. In the former case (the soft overlap for active-active Glauber exchange), this simply is the soft dynamics encoded in the soft Wilson line structure and Lagrangian. For the latter case, it corresponds to a soft emission off a collinear field, which does not appear in the SCET II Lagrangian. Instead it has been integrated out and is thus accounted for in the matching of QCD to SCET, where a soft Wilson line is introduced for every collinear Wilson line. This soft Wilson line then accounts for the Lipatov emissions off the Glauber subset of the corresponding collinear Wilson line (besides other soft effects). This leaves us with Glauber exchanges between spectators only, which also have at least one Lipatov emission. It is clear that the lowest order that such diagrams could appear in the perturbative expansions is O(α 2 s ). Focusing on Glauber exchanges between the proton remnants, we naively expect that a contribution of the form of figure 2(b) could appear at O(α 2 s ), when interfering with a suitable tree-level diagram. However, such interference is prohibited by momentum conservation: The proton remnants have transverse momenta of O(Λ QCD ), the typical scale of intra-proton dynamics, while the Glauber exchange represents a recoil between the involved spectators, causing them to have transverse momentum of O(q x ). This does not match up with any tree level conjugate diagram. Accordingly, the conjugate amplitude also requires some perturbative exchange (Glauber or otherwise), which pushes the contribution to O(α 3 s ). In addition, such diagrams do not contribute to our observable, as we only measure properties of the jet: Any recoil of beam remnants is never measured, and accompanying soft emissions are not relevant, because we use the WTA axis. For a diagram to have any effect it must impart recoil on the active partons or the jet constituents: The WTA axis follows a collinear emission in the jet, which can only pick up Glauber effects by direct recoil, or by inheriting recoil from the partons initiating the hard scattering.
An example for an allowed diagram would be figure 2(c), where a perturbative splitting of an initial state gluon creates spectators that already have transverse momentum of O(q T ), such that a Glauber exchange does not fall prey to momentum conservation. Such diagrams appear, by virtue of the initial state branching, at the earliest at O(α 4 s ). Lastly, we turn to diagrams involving the jet, like e.g. figure 2(d). As stated above, spectators in the jet can arise from perturbative splittings of the parton initiating the jet, which means such topologies start at O(α 3 s ) at the lowest. Momentum conservation does not pose an obstacle here, as the soft emission across the cut can attach to the beam spectator in the conjugate amplitude, pushing its transverse momentum to the right scaling. However, an explicit calculation of the loop diagram in figure 2(d) shows that it evaluates to zero, as may have been expected due to its similarity with deep-inelastic scattering. Explicitly, as only two rapidity sectors are involved, we can boost to a frame in which they are back-to-back. Of the four non-Glauber propagators involving the loop momentum k, two will then depend only on n · k, and two only onn · k, after expansions around the momentum scaling for Glauber and collinear modes. These two subsets correspond to the two collinear sectors, i.e. two propagators arise from the jet, the other from the beam remnant and active line. The two propagators arising from the jet have poles on the same side of the real line, and -18 -

JHEP02(2023)256
the loop integral thus evaluates to zero by residues. 10 Numerator factors and the rapidity regulator do not change this outcome.
We thus conclude that effects of perturbative Glauber exchanges are suppressed by at least O(α 3 s ) compared to the accuracy achieved in this paper.

Gluon jet function calculation at order α s
Next we present the calculation of the one-loop gluon jet functions. According to eq. (3.19), they are defined through where the transverse momentum of the jet with respect to the initial parton is encoded in In our original coordinates p x,J = 0, but we will instead perform our calculation in a frame where the total collinear momentum p x,c = 0. In this subsection we drop the subscript J on light-cone coordinates for brevity. As discussed in section 2.1, the difference between p T,J and p T,c is O(λ), so we can replaceˆ p T,J → p T,c , and similarly we can replace y J by y c . Then, following the same steps as for the standard jet axis in ref. [25], we can switch to coordinates along the momentum of the initial parton, to obtain Using the projectors in eq. (3.21), one finds at tree-level All one-loop diagrams for the gluon jet functions are shown in figure 3. The relevant Feynman rules for B µ n⊥ = 1 g W † n iD µ n⊥ W n are given by

5)
10 See in particular the discussions surrounding equations (11.9) and (B.6) in [92]. where g is the renormalized coupling withμ 2 = µ 2 e γ E /(4π), and we use the η-regulator [61,62] to regularize rapidity divergences. The w in eq. (4.5) is an (artificial) coupling used to derive the corresponding rapidity evolution. Diagram (d) vanishes in Feynman gauge, and the others contribute to the jet function according to, where the 1/2 is an identical particle factor and with p 1 · p 2 = p 2 ⊥ /(2z (1 − z)). The two-body collinear phase space entering in eq. (4.6), is defined as which we rewrite in terms of the transverse momentum p ⊥ and momentum fraction z, The amplitudes in eq. (4.7) are also expressed in terms of these variables. For the WTA scheme, one has (4.10) -20 -

JHEP02(2023)256
The evaluation of J i g involves the following two integrals For J T g , one has (4.12) By using eq. (4.8) and eq. (4.11), this leads to and the jet function renormalization is From this, one can easily obtain the one-loop anomalous dimensions which agree with their all-order form in eq. (5.3). For J L g , one only needs to include the terms ∝ p µ ⊥ p ν ⊥ in eq. (4.7). It can then easily be evaluated by using the two integrals in eq. (4.11), yielding -21 -

18)
J L g is finite at this order (as required, since it vanishes at tree-level). It is given by

Recombination scheme dependence
The WTA algorithm is not the only recombination scheme that can be used to construct a recoil-free jet axis. In this subsection we will employ a more general recombination scheme. It turns out that this only changes the finite part of the jet function, and at the end of this section we give explicit results. A recombination scheme dictates how particles are merged during the clustering procedure. The simplest is to add the momentum four-vectors, which is known as the E-scheme. In the rest of this paper we focus on the WTA-scheme, specifically the WTA-p Tscheme, described in section 2.1. Here, we consider a generalization in which the momenta p i and p j of two particles (or pseudojets) are recombined into p r as follows 20) in terms for the transverse momentum p T , azimuthal angle φ and rapidity y. The azimuthal angle and rapidity are combined in a way that favors the direction of the harder particle, with the weight of the factors p n T,i controlled by the power of n ≥ 1. For n > 1 this recombination scheme is recoil free, implying that the same factorization theorem holds as for the WTA scheme (which corresponds to n → ∞). Since only the jet function depends on n, consistency of the factorization implies that only the constant term of the jet function depends on it, which is borne out by an explicit calculation. Now, let us calculate gluon jet functions at one loop using the general scheme in eq. (4.20). For two collinear partons, using coordinates along the parton that initiates the jet, one has φ i ≈ p x,i /p T,i and p x = p x,1 = −p x,2 such that 11 in terms of the kinematic variables in eq. (4.9). As a consistancy check, one can easily see that δ x reduces to the expression for the WTA axis in eq. (4.10) in the limit n → ∞. With the following replacementẑ one can straightforwardly evaluate the gluon jet functions for the general recombination scheme by extending the calculation in the previous section. For J L g , one can easily see that eq. (4.19) is valid for all the values of n because the dependence onẑ n drops out in eq. (4.18). J T g is given by the first two lines of eq. (4.13) with the replacement forẑ n in eq. (4.22). The 1/η pole arises only from the following expansions (4.23) in the integrand of eq. (4.13). Since theẑ n -dependence drops out and the 1/η term is the same as for the WTA Scheme. Similarly, the only source for the 1/ pole in J T g comes from the pole in I 1 with Hence, the jet function renormalization in eq. (4.15) and the anomalous dimensions in eq. (4.16) are valid for any value of n > 1, as required by consistency of the factorization. By using the expansion in eq. (4.25), we obtain the constant term of J T g (i.e. taking L b = 0)

JHEP02(2023)256
For n → ∞, it reproduces that for the WTA scheme: One can also find analytic results for specific values of n. E.g. for the p 2 t scheme (n = 2), we have . (4.28) The two constants J C A and J n f are plotted as function of n in figure 4. The analytic result for the WTA scheme is also shown, and one can see that it is approached in the limit n → ∞.
For n = 1 the jet is no longer recoil-free. This is reflected in a diverging constant in the n → 1 limit, indicating that the poles of the jet function, and thus the entire factorization are different for n = 1. Similarly, the finite part of J q in general depends on n and takes the form In figure 4, J C F is also shown. For the p 2 t scheme and for the WTA scheme,

Track-based jet function
Next, we will consider the case where the jet is measured using only charged particles, exploiting the superior angular resolution of the tracking system. Since this does not modify the effect of soft radiation on the measurement, which still contributes through the total recoil, this only affects the jet function. This is similar to the different recombination schemes considered in section 4.2, and we can reuse much of the calculation, though we will restrict ourselves to the effect of a track-based measurement for WTA scheme. In particular, the jet function anomalous dimensions are not modified but only the constant, which should be contrasted with the complicated jet function encountered for a track-based measurement of thrust in ref. [32]. We will account for the conversions of the partons to charged particles using the track function formalism [31,32]. At one-loop order, the jet consists of (at most) two partons and for q x Λ QCD they can be treated as fragmenting independently into charged hadrons moving in the same direction as the original partons. Denoting the total momentum fractions of charged hadrons produced by each of the two partons by z 1 and z 2 , the only change due to a track-based measurement is that the condition which parton "wins" gets modified:ẑ . (4.31) We also need to take into account the nonperturbative distribution of z 1 (and z 2 ) which is described by a track function T f (z 1 , µ), where f is the flavor of the parton. For example, for

JHEP02(2023)256
the gluon jet function in eq. (4.6), the corresponding track-based measurement (indicated by the bar) is where |δ x,ch | = |p x |/ẑ ch , and the track functions depend on the flavor of the partons in the final state. Note that different quark flavors have different track functions, but T q = Tq due to charge conjugation invariance. The subsequent steps directly parallel those in section 4.2, so we find again that the one-loop jet function for the linearly-polarized gluonJ L g is not modified, while for the transversely-polarized gluon and quark we havē (4.33) in direct analogy to eqs. (4.26) and (4.29).

Renormalization group evolution
In section 3 we have given the derivation of the factorization formula and the explicit expressions for the one-loop ingredients. Within the EFT framework one then uses the renormalization group (RG) evolution equations to resum the large logarithms between different scales. In addition to the standard UV divergences regularized by the dimensional regularization, the jet, beam and soft functions also involve rapidity divergences, as these modes are not separated in invariant mass but rapidity, see eq. (2.6). In order to resum the corresponding rapidity logarithms we apply the rapidity RG method developed in refs. [61,62]. Different regulator choices (e.g. [80,[95][96][97][98][99]) or resummation via the collinear anomaly framework [65,100] are also possible and (up to the possibility of scale variation) equivalent. To achieve next-to-next-to-leading logarithmic resummation, we need the oneloop fixed order ingredients in section 3.3 and the two-loop anomalous dimensions (except for the cusp term in the anomalous dimensions, which is required at three loop order). Generally, the RG equations for a function F are given by where Γ F µ and Γ F ν denote the standard and rapidity anomalous dimensions, respectively. For the beam, jet and soft function, this multiplicative form of the evolution equation only -25 -

JHEP02(2023)256
holds in impact parameter space, and the anomalous dimension depends on b x . In addition the anomalous dimension may e.g. depend on the hard kinematics or the jet direction, which we omit.
The anomalous dimension of the hard function are 2) where C i and γ i µ are the color factor and non-cusp anomalous dimension for parton i, i.e. C A and γ g µ for a gluon, and C F and γ q µ for an (anti-)quark. As we always have (up to permutation) two quarks and a gluon in our process, the non-cusp anomalous part simplifies to a γ a µ = 2γ q µ + γ g µ . The expressions for the partonic Mandelstam variables in terms of the kinematics of the boson and jet are collected in eq. (3.33). The perturbative expansion for Γ cusp and γ a µ have been collected in appendix A. The anomalous dimensions of the beam functions are where ω i is explicitly defined in eq. (3.26), , and the cusp anomalous dimension Γ cusp , non-cusp anomalous dimension γ B i µ , and rapidity anomalous dimension γ ν in appendix A. The function A Γcusp is obtained by replacing γ i → Γ cusp in A γ i in eq. (A.10). The jet function J satisfies the same anomalous dimension as the beam function, where now ω J =n J · p J = 2p T,J cosh η J .
The anomalous dimensions of the soft function are From eqs. (5.2), (5.3), (5.4) and the non-cusp anomalous dimensions in appendix A, it is straightforward to verify that the factorized cross section in (3.8) is independent of the renormalization scale µ and rapidity renormalization scale ν: Figure 5. We evolve all ingredients from their natural (µ, ν) scale to the scale of the soft function along the indicated paths. We take the different ν i into account though they are of the same parametric size.
Using the renormalization group equations to evolve all ingredients in eq. (3.8) from their natural µ and ν scales to a common scale, the all-order resummation formula can be written as where we use i to label the parton flavor but also the rapidity scales of the beams and jets, and understand Γ B k ν to refer to the jet function rapidity anomalous dimension. We chose to evolve the beam and jet function from their natural rapidity scales ν i to the rapidity scale ν S of the soft function at the common invariant mass scale µ B , and evolve the hard function from µ H down to the scale µ B , as summarized in figure 5. The evolution factor can be evaluated analytically as The functions A and S are given in appendix A. The natural scales for the various ingredients in the resummation formula are taken to be where we avoid unphysical results in the large-b region, by applying the b * -prescription [64] b In figure 6 we show the resummation results at NLL and NNLL accuracy, separately displaying the various contributions to the perturbative uncertainties. Specifically, we assess the uncertainties by varying µ H , µ B and ν S up and down by a factor of two around their default values in (5.9). The fact that the uncertainty band is smaller at NNLL than at NLL, and that the bands overlap over almost the entire range, suggest that this is a reasonable estimate. In our full result we will combine these uncertainties by taking the envelope.
To model non-perturbative corrections, we furthermore will include the multiplicative function e −S NP (b) We take the same nonperturbative function S NP as for transverse momentum distributions in [101], assuming that the nonperturbative contribution to the rapidity anomalous dimension (with coefficient g 2 ) can be obtained for a gluon beam or jet function by Casimir scaling. It is not clear whether this should also be the case for the nonperturbative model (with parameter g 1 ) and we take it the same for quark and gluon beams/jets, finding that it has a negligible effect on numerics anyway. We will use the results for the nonperturbative parameters obtained in ref. [101]: Q 2 0 = 2.4 GeV 2 , g 1 = 0.212 GeV 2 and g 2 = 0.84. The sensitivity of our predictions to these nonperturbative parameters is explored in figure 7, finding minimal sensitivity to g 1 but sensitivity to g 2 at the percent level in the region of ∆φ ∼ π.

Matching to fixed-order MCFM
In the back-to-back region where δφ → 0 (∆φ → π), the resummation formula will cure the divergence behavior of fixed-order results. However, if δφ is not small, the factorization formula receives large corrections of powers of δφ. In this region, the resummation should be switched off, since ln δφ is no longer large, and we therefore need to use fixed-order calculations that include these power corrections.
We use an additive matching scheme, in which the "naive" matched result of NNLL resummed prediction and the fixed-order can be obtained by the following relation dσ naive (NLO + NNLL) = dσ(NNLL) + dσ(NLO) − dσ(NLO singular) dσ(NLO non−singular) , (5.12) where we use MCFM [102,103] to calculate the NLO results. The NLO singular distribution removes the overlap between the first two terms and can be obtained by expanding the resummation formula (5.7) to O(α s ). The NLO non-singular distribution is given by the difference between NLO and NLO singular results, as indicated in the above definition.
In principle, as δφ becomes large, the NNLL resumed cross section reduces to the NLO singular, leading to a cancellation between the first and third term in eq. (5.12). However, as is clear from (5.7), the numerical Fourier transformation from b-space to the momentum space rapidly oscillates when the resummation turns off and the evolution factor approaches 1. To avoid the corresponding numerical instability, we apply a transition function t(∆φ) as follows: The transition function we use is defined as  figure 9. The nonsingular is in principle much smaller than the singular, but because the singular is resummed, there is less of a difference between them. Therefore we can see a sizable effect on how the nonsingular correction is treated in the resummation region, particularly when the jet p T is large. We will discuss the origin of this large nonsingular correction in section 6.2, arguing that it should not be Sudakov suppressed in the back-to-back limit, which is why we use an additive rather than multiplicative matching. We conclude this section by presenting expressions for the NLO singular at the level of the integrated cross section, At order α s we have ij←ab . The first term on the second line is the contribution from the unpolarized gluon beam/jet functions, while the second term is the linearly-polarized contribution, as indicated by the superscript L. The one-loop coefficients are given by where L = ln δφ cut , the splitting functions are given in eq. (3.27), and For the different partonic channels, the coefficients A ij are given by where the partonic Mandelstam variables are given in terms of the kinematics of the hard scattering in eq. (3.33).

JHEP02(2023)256 6 Results
We start in section 6.1 with a detailed study of the azimuthal angle between a recoilfree jet and Z boson using the Pythia Monte Carlo parton shower. This allows us to investigate the effect of hadronization and underlying event, and corroborate conclusions of our factorization analysis with regards to the dependence on the jet radius, recombination scheme and track-based measurements. Our resummed predictions are shown in section 6.2.
We also explain sizable non-singular corrections, particularly at large p T,J , which may be largely removed by boson isolation cuts.

Monte Carlo analysis
In this subsection we present a phenomenological study of recoil-free boson-jet correlation using the Pythia 8.3 [104] Monte Carlo parton shower. The Z+jet events in 13 TeV proton-proton collisions at the LHC are simulated with the decay of the Z boson turned off. In experiments the clean, leptonic decay channels of Z boson are reconstructed with suitable cuts on the lepton kinematics. In these Monte Carlo studies we sum over all the Z boson polarization states to match our analytic calculation. In all events, jets are reconstructed using the anti-k t algorithm [105] with R = 0.5 (also R = 0.8 or R = 1.0 when studying the jet radius dependence) using FastJet 3 [41] and |η J | < 2. The azimuthal angle is defined as the one between the Z boson and the leading jet in each event. 12 We consider two kinematic regions: p T,J > 60 GeV and p T,J > 200 GeV, to study the dependence of this observable on the hard energy scale. Two million events for each region are simulated, providing sufficient statistics to obtain smooth distributions.
We first examine the sensitivity of the azimuthal decorrelation to hadronization and underlying event in the left panel of figure 10, which shows the ∆φ distributions with or without hadronization or underlying event contributions. In Pythia the underlying event is modeled as multi-parton interactions (MPI). We see that the shape of the ∆φ distribution is remarkably insensitive to hadronization and MPI, which suggests that it is dominated by perturbative contributions. This is expected: due to our recoil-free jet definition, these soft contributions do not interfere with the jet finding, and only provide a total recoil of the V +jet system. Since the azimuthal angle is a vector quantity, the net effect of this recoil tends to be (close to) zero. There is a change in the normalization of the absolute cross section, because the additional radiation affects the number of events having a jet with sufficient transverse momentum.
In order to exploit the high angular resolution of charged particle tracking, we study the case where the recoil-free axis is determined only by charged tracks, and the difference compared to using all charged and neutral particles is examined. This is shown in the right panel of figure 10. As can be seen, the azimuthal decorrelation distributions using the recoil-free axis determined by charged particles or all the particles within jets are almost identical, in line with the our conclusion in section 4.3 that this difference is beyond NLL accuracy. To contrast the WTA axis choice, we provide distributions for jets defined using the more general p n t recombination scheme in eq. (4.20). In the left panel of figure 11 we examine the case where n = 1 (the standard), n = 2 and n → ∞ (WTA). Since the WTA axis is sensitive to momentum-conserving collinear splittings within jets and following the energetic branch, the ∆φ distribution is broader near the back-to-back region. The difference between the recoil-free n = 2 and WTA axis is beyond NLL order, as discussed in section 4.2, and is indeed small. The difference with the recoil-sensitive case n = 1 is larger. We also compare the distributions for jets with different radii as a way to highlight again the insensitivity to soft, typically wide angle, radiation, in the limit where the jet radii are much larger than the azimuthal decorrelation, δφ R. The right panel of figure 11 shows the ∆φ distributions for jets reconstructed using different jet radii, namely R = 0.5, R = 0.8 and R = 1.0. The shapes of the distributions are very similar as expected from our factorization analysis. The differences for ∆φ in the vicinity of 180 • arise because these distributions are normalized. However, some normalization is necessary when comparing different jet radii, because the jet radius has a non-negligible (5-10%) effect on the normalization of the cross section through the cut on the jet transverse momentum.

Resummed predictions
In this section we present numerical results from the resummation formula in eq. (5.7), which we match to MCFM [102,103]  m V there is also a collinear singularity at ∆φ 12 → 180 • from a Z boson emitted from a dijet partonic configuration. singular 1/(δφ) terms from expanding around the back-to-back limit of the boson and jet, missing an important contribution to the power corrections that arises from a dijet configuration where a Z boson is emitted from one of the jets. This contribution enters at the same order in the coupling and is enhanced for p T,J m Z , as will be discussed below. To illustrate the difference between the contributions in the factorization theorem and this important nonsingular correction, consider the diagrams for the gg → qqZ process in figure 13: Diagrams (a) and (b) contribute to our factorization theorem in the region of phase-space where a final-state quark is collinear to one of the incoming gluons. However, diagram (a) and the new diagram (c) are also enhanced for the region of phase space describing a Z emission of a dijet configuration. This contribution is included by the matching procedure in section 5.2 and discussed more below. . Accordingly, at low p T , the Z emission off dijets is expected to be independent of the azimuthal angle based on the above expansion. Indeed, the non-singular in figure 12 is almost constant at small δφ. This is also confirmed by looking at the azimuthal angle ∆φ 12 of the two final-state partons shown in figure 14. We find that for p T,J m V (e.g. the 65 GeV> p T,J > 55 GeV bin) the ∆φ 12 distribution is flat near ∆φ 12 = 180 • . This flatness is the reason that we do not expect the nonsingular to go to zero at ∆φ → 180 • , and don't employ a multiplicative matching that would enforce that. Indeed, collinear and soft emissions will smear the nonsingular δφ distribution, but since it is (almost) constant it will not change (much). Figure 14 also shows that at high p T,J there are large contributions from the Z emission off back-to-back dijets, in addition to soft-collinear QCD radiation for back-to-back boson-jet production. The latter is the focus of this paper, incorporated in the resummation formula eq. (5.7), while the former is included by matching and needed at NNLL and beyond. When p T,J m V , it is insufficient to include this contribution by matching in the region m V /p T,J δφ 1, as it now also contains large logarithms of δφ that require resummation. Alternatively, one can also remove the contribution of a Z boson emitted from a dijet by introducing an isolation cone around the boson (or its decay products).
We conclude by showing in figure 15 our resummed predictions at NNLL+NLO accuracy, and comparing to Pythia results at hadron level with MPI contributions for the two p T,J regions. We show the resummed predictions with and without matching, in view of the large power corrections we just discussed, and also include the NLO cross section as a separate curve. The NLO cross section divergences as ∆φ → 180 • , which is remedied by 13 Including numerator factors it only gives rise to a NLO term ∝ 1/δφ in the cross section. the resummation. Our resummed results without matching agree with the shape of Pythia simulations. By including the matching, our predictions smoothly approach NLO for smaller values of ∆φ, where resummation is not needed. The large nonsingular correction (larger than our uncertainty bands) is important to include and cannot be neglected as ∆φ → 180 • . It is not accounted for in the Pythia simulations we used. Simply attempting to include it through a K-factor does not yield the correct shape, particularly at high p T,J .

Conclusions
In this paper we describe our calculation of the cross section for a vector boson and jet, differential in the azimuthal decorrelation δφ, at next-to-next-to-leading logarithmic order. We provide substantially more details than our earlier work [30], an expanded phenomenological analysis, and also include a Pythia study. We discuss for the first time why the azimuthal angle is simpler than the total transverse momentum q T , potential Glauber contributions, different recoil-free recombination schemes, the jet radius dependence, and the non-singular matching.
Our focus is on the region δφ 1, which dominates the cross section and requires the resummation of logarithms of δφ to obtain reliable predictions. We carried out a detailed study of its factorization in SCET, deriving a factorization formula. We investigated potential factorization-violating Glauber contributions, finding that they could first appear at order α 4 s . Many of the ingredients in the factorization are available, and we present calculations of the jet functions, including the linearly polarized jet function, as well as jet functions for p n T -weighted recombination schemes. We find that these different recombination schemes only change the constant term in the jet function, having a minimal effect on the prediction. Furthermore, we verify the independence on the jet radius predicted by the factorization (for R δφ) in Pythia. By using the (rapidity) renormalization group we achieve the resummation of logarithms of δφ.
The key to obtaining predictions beyond NLL is a recoil-free recombination scheme, which reduces the effect of soft radiation to a total recoil of the V +jet system, eliminating non-global logarithms. This is in contrast to the standard recombination scheme, or other transverse momentum measurements such as the radial decorrelation, which involve NGLs. A recoil-free recombination scheme is also interesting for studying the properties of the medium in heavy-ion collisions, because it is more sensitive to collinear splittings inside the jet and less sensitive to contamination of soft radiation. Indeed, we observe negligible effects of hadronization and MPI contributions in Pythia, as well as compatibility with track-based measurements. Our resummed predictions are matched to MCFM in the region where δφ is no longer small. Here we find that the non-singular contributions are quite large, particularly for high jet p T , increasing the sensitivity to the details of the matching procedure. It would be interesting to investigate the resummation of these power corrections.
We have established the azimuthal angle between a vector boson and a recoil-free jet, as a robust observable for which high precision is possible. In this paper we achieve NNLL order, but NNNLL is within reach. As the effect of the underlying event, hadronization and from performing the measurement using charged particle tracks is minimal, it is attractive experimentally.

A Anomalous dimensions
In general, for the (next-to) i -leading logarithmic (N i LL) resummation, one needs up to (i − 1)-loop fixed-order ingredients, i-loop non-cusp anomalous dimensions and the (i + 1)loop cusp anomalous dimension and QCD beta function. In this appendix, we collect the beta function and all the anomalous dimensions that are needed for this paper.
For the beta function, one has up to three loops [108,109] where n f is the number of active quark flavors. Identifying standard anomalous dimensions by their associated functions, and using the universality of the rapidity anomalous dimension, we write the perturbative expansions of cusp, non-cusp, and rapidity anomalous dimensions as The cusp anomalous dimension is, up to three loops, given by [110,111] Γ 0 = 4 , T F n f ,

JHEP02(2023)256
For the qq → V g channel, we have C qq→V g (t, u) = C A π 2 6 + C F −16 + 7π 2 3 + 2C A ln 2 s m 2 and for the qg → V q channel we have

JHEP02(2023)256
with the function T 0 defined as Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited. SCOAP 3 supports the goals of the International Year of Basic Sciences for Sustainable Development.