Two-loop mixed QCD-electroweak amplitudes for $Z+$jet production at the LHC: bosonic corrections

We present a calculation of the bosonic contribution to the two-loop mixed QCD-electroweak scattering amplitudes for $Z$-boson production in association with one hard jet at hadron colliders. We employ a method to calculate amplitudes in the 't Hooft-Veltman scheme that reduces the amount of spurious non-physical information needed at intermediate stages of the computation, to keep the complexity of the calculation under control. We compute all the relevant Feynman integrals numerically using the Auxiliary Mass Flow method. We evaluate the two-loop scattering amplitudes on a two-dimensional grid in the rapidity and transverse momentum of the $Z$ boson, which has been designed to yield a reliable numerical sampling of the boosted-$Z$ region. This result provides an important building block for improving the theoretical modelling of a key background for monojet searches at the LHC.


Introduction
The production of a Z boson in association with hadronic jets is a key standard candle at the Large Hadron Collider (LHC).Thanks to its large production rate and its relatively clean dilepton plus jet final state, it allows for a multitude of investigations like detector calibration and luminosity monitoring, validation of Parton Shower Monte Carlo tools, extraction of fundamental Standard Model (SM) parameters, as well as high-precision scrutiny of the structure of the SM and searches for new physics.
Despite the fact that the bulk of the Z+jet cross section at the LHC comes from a region where the jet's transverse momentum is not very large, p T,j ≲ 100 GeV, the number of events where the Z boson is accompanied by a highly-energetic jet is still quite sizable.This allows for very precise measurements all the way up to the TeV scale.Indeed, existing experimental analyses [1,2] already have a total uncertainty of just a few percent in the p T ∼ 200 GeV region and around 10% in the highly-boosted p T ∼ 1 TeV region.With ever more data being recorded and analysed, the situation is only going to improve: at the High-Luminosity LHC, few-percent experimental precision is expected up to transverse momenta of the order of 1 TeV, and O(10%) precision up to 2.5 TeV [3].
A good control of Z production in the boosted region allows for interesting physics explorations.For example, precise data in the dilepton channel can constrain otherwise elusive dimension-8 Standard Model Effective Theory (SMEFT) operators [4], provided that adequate theoretical predictions are also available.Also, boosted Z-boson production in the Z to invisible channel provides a key background for monojet searches at the LHC, where one looks for a high-p T jet recoiling against missing energy [5,6].Such a signature is very interesting because it is quite common in many new physics models, ranging from weakly coupled dark-matter, to leptoquarks models, supersymmetric scenarios, or large extra dimensions.
In the recent past, there has been a large community-wide effort to improve the theoretical description of boosted vector-boson production.In particular, in ref. [3] the authors provided theoretical predictions that include state-of-the-art NNLO QCD results [7][8][9][10][11][12][13][14][15] 1 and NLO electroweak (EWK) ones [17][18][19][20].In the boosted region, the latter are crucial.Indeed, despite being suppressed by the weak coupling constant α, they are enhanced by large Sudakov logarithms of the form α , where s w is the sine of the weak mixing angle, s is a large scale of the process and m V is the vector boson mass, see e.g.[21].At large scales, EWK corrections then become as large as QCD ones.For this reason, ref. [3] also included the dominant two-loop electroweak effects coming from Sudakov logarithms [22][23][24][25].
Given the size of QCD and EWK corrections, it is also mandatory to properly control mixed QCD-EWK ones.At present, O(α s α) corrections to dilepton+jet production are not known.In ref. [3], the size of these corrections was estimated by essentially multiplying NLO QCD and NLO EWK ones.This prescription is very reasonable at asymptoticallyhigh scales, since Sudakov logarithms and QCD corrections mostly factorise.However, at large but finite energies this approximation is bound to receive corrections.A more rigorous assessment of mixed QCD-EWK corrections then becomes important.In fact, the lack of exact O(α s α) corrections is now a major bottleneck towards highest-precision theoretical predictions in the boosted region [3].
Computing O(α s α) corrections to dilepton+jet or missing-energy+jet production at the LHC poses significant challenges.First, such processes involve a non-trivial mixed QCD-EWK radiation pattern, which requires proper regularisation.Achieving this is complicated if one wants to retain differential information on the final state.This problem has only recently been solved, and only for the simplest processes [26][27][28][29][30][31].Although the techniques used for the calculation [31] could be extended to deal with more complex processes, this remains a non-trivial task.Second, mixed QCD-EWK corrections for dilepton+jet or missing-energy+jet production require non-trivial two-loop scattering amplitudes.These involve both a complex final state and massive internal virtual particles, which makes the calculation notoriously difficult.
In this article, we perform a first important step towards the calculation of mixed QCD-EWK two-loop scattering amplitudes relevant for boosted dilepton+jet production.To make the problem manageable, we focus on the production of an on-shell Z boson rather than on the production of the dilepton final state, with the idea that this is going to provide the dominant contribution for all observables which are dominated by the Z-pole region (and in particular for the observables relevant for the boosted region).This allows us to only consider amplitudes for 2 → 2 scattering rather than the much more complicated ones relevant for the 2 → 3 process.Also, we only target the boosted region, where vector-boson resonance effects are not present and hence use of the complex-mass scheme [32][33][34][35][36] for EWK corrections becomes less critical. 2Finally, as a first non-trivial step towards the full result here we only consider bosonic corrections, i.e. we systematically neglect closed fermion loop corrections.This allows us to set up a framework for computing mixed QCD-EWK corrections without the additional complication of corrections involving top-quark virtual effects, which cannot a-priori be neglected in the boosted region. 3We believe that our framework could be extended to cover this case as well, but this warrants an investigation on its own.Even in this somewhat simplified setup, an analytical calculation of the amplitude remains challenging.Because of this, we decided to adopt a semi-numerical approach.Our main result is then an evaluation of the two-loop amplitudes over a twodimensional grid that parametrises the 2 → 2 kinematics.The grid is designed to provide an adequate coverage of the boosted region.
The remainder of this paper is organised as follows.In sec. 2 we provide details of our notation and of the kinematics of the process.In sec. 3 we describe some methods used in our work i.e. the Lorentz tensor structure of our scattering amplitude in sec.3.1, and its relation to the leptonic current at fixed helicity in sec.3.2.In sec.3.3 we describe our calculation of the bare amplitudes.In sec. 4 we discuss the ultraviolet and infrared structure of our result and we also define one-and two-loop finite remainders.The latter are the main result of our paper.In sec.5 we document the checks that we have performed on our calculation and illustrate our results.Finally, we conclude in sec.6.Our results for the finite remainders are available in a computer-readable format in the ancillary material that accompanies this submission.

Notation and kinematics
We consider virtual O(αα s ) corrections to the production of the Z boson in association with one hadronic jet.We focus on bosonic corrections, i.e. we neglect contributions stemming from closed fermion loops.We then consider the channels where q is either an up-or a down-type (anti) quark.We take the CKM matrix to be diagonal, and neglect b-quark induced contributions.This way, no virtual top-quark contributions are present.All the processes in eq.2.1 can be obtained by crossing the following master amplitude 2 For a discussion of the complex-mass scheme and of its importance for EWK radiative corrections, see e.g. the recent review [37] and references therein. 3For the same reason, we do not-consider b-quark induced contributions.
where q is either an up or a down quark.In this symmetric notation, momentum conservation reads All external particles are on-shell (2.4) The three kinematic Mandelstam invariants of this process are related by the momentum-conservation relation For physical kinematics, one Mandelstam invariant is positive and two are negative.Results in the Euclidean region where m 2 Z < 0 and s ij < 0 can be analytically continued to the physical Riemann sheet by giving a small positive imaginary part to all the invariants, see ref. [38] for a thorough discussion.
For definiteness, we now focus on the s 12 > 0, s 13 , s 23 < 0 channel.In the partonic center-of-mass frame, the kinematics can be parametrised as follows (1, 0, 0, −1) , p with where s had is the collider center-of-mass energy squared.In what follows, we will either use {s 13 , s 23 } or {p t,Z , y Z } as independent variables.
To deal with the ultraviolet (UV) and infrared (IR) divergences of the amplitude, we work in dimensional regularisation and set d = 4 − 2ϵ.In particular, we adopt the 't Hooft-Veltman scheme [39], i.e. we treat all external particles as purely four-dimensional and the internal ones as d = (4 − 2ϵ)-dimensional.We write the unrenormalised amplitude as where i 1 , i 2 , and a 3 are colour indices of the quark, antiquark, and gluon, respectively, ϵ 3,µ and ϵ 4,ν are the polarization vectors of the massless gluon and of the massive Z boson, m i,b , i ∈ {Z, W } are the (bare) vector boson masses, and α b , α s,b are the bare electromagnetic and strong couplings, respectively.The dependence on other electroweak parameters like the quark isospin or electric charge is implicitly assumed.Finally, T a ij is the SU (3) generator in the fundamental representation, rescaled such that (2.11) We find it convenient to express our results in terms of the quadratic Casimir invariants C A and C F , which for SU (N c ) read In QCD, C A = 3 and C F = 4/3.The amplitude in eq.2.10 can be written as a double perturbative series in the strong and electromagnetic couplings where we have made explicit the dependence on the dimensionful reference scale µ 0 but have kept the dependence of A b and A (i,j) b on the vector boson masses and on the external kinematics implicit.
To obtain the renormalised amplitude A, we multiply A b by the external wave-function renormalisation factors and express the generic bare parameter g i,b in terms of its renormalised counterpart g i .Schematically, We discuss the renormalisation procedure in more detail in sec.4.Here we only note that we renormalise the strong coupling in the MS scheme, and all the EWK parameters in the onshell scheme.Also, we adopt the G µ input parameter scheme, i.e. we choose {G µ , m Z , m W } as independent parameters.For our results, we use G µ = 1.16639 × 10 −5 [40].Similarly to eq. 2.10, we write the renormalised amplitude as and expand A as In these equations, m i are the renormalised masses, α is the renormalised electromagnetic coupling and α s = α s (µ R ) is the renormalised strong coupling with µ R the renormalisation scale.The main result of this article is the computation of A (1,1) .

Details of the calculation
For further convenience, we provide here some details on the formalisms used for our calculation.

Tensor decomposition
In full generality, the bare scattering amplitudes in eq.2.10 can be written as where are scalar form factors and , where Γ µν k are n t independent Lorentz tensors.We now discuss a basis choice for these tensors.First, we note that in this paper we do not consider corrections coming from closed fermion loops.As a consequence, there are no anomalous diagrams and we can take γ 5 as anticommuting.
Because of this, it is enough to study the tensor structure of a vector current, and the axial case will follow.For a vector current, by simple enumeration one finds n t = 39 independent Γ µν i structures.Using the Dirac equation for the external on-shell quarks v2 / p 2 = / p 1 u 1 = 0, as well as the transversality condition for the massless gluon ϵ 3 • p 3 = 0, n t decreases by 26.Finally, if one also uses the Z-boson transversality condition ϵ 4 • p 4 and fixes the reference momentum q 3 for the polarisation vector ϵ 3 to the momentum of one of the external fermions, one is left with n t = 7 independent d-dimensional tensor structures.
Following ref. [16], we set q 3 = p 2 , such that and define the first six tensors as follows4 These six structures span the the whole space of purely four-dimensional tensors [41,42].We can then construct the last structure to only have components in the unphysical (−2ϵ)-dimensional space.This can be done through a simple orthogonalisation procedure [41,42], which yields Since in strict d dimensions the seven tensors are independent, there cannot be any cancellation of IR and UV singularities among the different T i .Hence, the renormalised and IR-regulated form factor F (i,j) 7 cannot have any poles.Since the corresponding tensor T 7 is evanescent by construction, the contribution of F (i,j) 7 drops from the final result if one works in the 't Hooft-Veltman scheme and sets d = 4 after the renormalisation and IRregularisation procedure.Therefore, in the 't Hooft-Veltman scheme, all physical results can be obtained from the first six tensor structures alone [41,42].We note that the number of independent tensors coincides with the number of independent helicity states of the external particles, i.e. it matches the independent four-dimensional degrees of freedom.A relation between the six tensor structures in eq.3.3 and the independent helicity amplitudes is discussed in the next sub-section.
Since the six tensors in eq. ( 3.3) span the physical space, it is always possible to find coefficients c ik such that In fact, solving these equations for c ik is a matter of trivial algebra.For convenience, we report them in appendix A. The projector operators P i defined in eq.3.5 can then be used to extract the F (i,j) k;b form factor from the bare amplitude Finally, we stress that the tensors 3.3 have been determined under the transversality condition p i • ϵ i = 0 and with the reference momentum choice q 3 = p 2 .It is then important to use the following expressions for the sum over polarisations of the external gluon and gauge boson We close this section by discussing how to generalise the above construction to the case of different left-and right-handed interactions.As we already mentioned, if we neglect closed fermion loops there are no anomalous diagrams and one can treat γ 5 as anticommuting.This makes such generalisation straightforward.The amplitude eq.3.1 can be decomposed as with and Γ µν i defined in eq.3.3.Projectors onto the left/right form factors F (l,m) i,b,L/R can be performed using the same procedure explained above, see eqs 3.5, 3.6, with the replacements i,L/R , and c ik → c ik /2.In practice, we always work with vector-current tensors and projectors, and then dress them with relevant left-and right-handed couplings.This allows us to effectively halve the various tensor manipulations on the amplitude.

Helicity amplitudes
The form factors F i introduced in the previous section can be used to obtain the helicity amplitudes for the process where l and l are a pair of massless leptons.In the pole approximation [43,44] and neglecting O(α 2 ) corrections, the amplitude can be schematically written as where "prod"/"dec" refers to the amplitude for the production/decay of an on-shell Z boson (p 2 Z = m 2 Z ) with polarisation λ Z .The numerator in eq.3.11 can be written in terms of the tensor structures introduced in sec.3.1 as with p Z = p 5 + p 6 = −p 1 − p 2 − p 3 , and p 2 Z = m 2 Z .In eq.3.12, c l,L , c l,R are generalised leftand right-handed couplings that parameterise the Zl l vertex.At LO, they are given by where I 3 l = ±1/2 is the weak isospin of the lepton l, Q l is its electric charge in units of e (Q e = −1 for an electron), and s w = sin θ w , c w = cos θ w with θ w the weak mixing angle.For mixed QCD-EWK corrections, one needs c l,R/L up to O(αα s ).These are well known, see e.g.refs [45,46], so we won't discuss them further.
The amplitude M in eq.3.12 can be easily projected to helicity states.To do so, we use the spinor-helicity formalism (see e.g.[47]) and write left-and right-handed currents as We also define the helicity of a particle/antiparticle to be equal/opposite to its chirality, and write the polarisation vector of the incoming gluon as with q 3 the reference momentum the gluon, q 3 • ϵ 3 = p 3 • ϵ 3 = 0. We remind the reader that our tensor T i have been constructed using the choice q 3 = p 2 .With these assignments, the helicity amplitudes M ⃗ λ , ⃗ λ = {λ 1 , λ 2 , λ 3 , λ 5 , λ 6 } for the all-incoming process eq.3.10 read with Because of this, we will focus on the latter in what follows.We conclude this section by pointing out that eq.3.16 is in agreement with eqns (5.6-5.17) of ref. [16] once all the small differences in notation are accounted for.

Computation of the bare amplitude
We now describe the calculation of the bare two-loop mixed QCD-EWK factors F. We start with generating, using Qgraf [48], all the relevant Feynman diagrams contributing to the bare amplitude A  we neglect closed fermion loops and bottom-induced contributions.We are left with 462 non-vanishing Feynman diagrams (FD) Comparing to lower orders, (see tab. 1) the complexity grows significantly, thus making the efficiency of the calculation crucial.For this reason, we first perform a useful parallelisation criterion, and only later the projection of the amplitude onto the tensor structures introduced in sec.3.1.To this end, we split the Feynman diagrams by their graph structure.
In the two-loop four-point amplitude, each Feynman diagram can have at most 7 virtual propagators.We will refer to diagrams containing a set of 7 different virtual propagators as top sector diagrams.Given the kinematics of the problem, there are 9 independent Lorentz invariants involving loop momenta k i , i = 1, 2. Out of these, 3 are of the type k i • k j and 6 of the type p i • k j .These 9 invariants can be linearly related to the 7 top sector propagators complemented with 2 additional irreducible scalar products (ISPs), which for convenience we also choose to be of the propagator type.We will refer to this set of 9 propagators as an integral topology.One can find a minimal set of integral topologies onto which all Feynman diagrams can be mapped.For our case, we require 18 basic topologies, as well as their 61 crossings.We report the definition of the 18 basic topologies in appendix B, and in computer-readable format in an ancillary file provided with this submission.Feynman diagrams with a smaller number of virtual propagators can be obtained by pinching top sector diagrams, which often makes their mapping onto an integral topology not unique.In practice, we perform this mapping by finding an appropriate loop momenta shift with Reduze 2 [49,50].As a result, we split the whole set of Feynman diagrams contributing to A Common structures of different integral topologies are typically revealed only after loopmomenta integration and/or expansion in ϵ.Hence, grouping Feynman diagrams by topology and performing preliminary manipulations separately for each topology provides an efficient parallelisation strategy.
To actually compute the amplitude, we work in the Feynman gauge.After substituting the Feynman rules, we immediately perform the colour algebra and write our result in terms of C A and C F , see eq. 2.12.Next, we project on the 12 form factors F 1,...,6,L/R;b defined in sec.3.1.To do so, we apply the P i projectors (see sec.3.1) and evaluate all required traces of Dirac gamma matrices in d-dimensions with Form [51], treating γ 5 as anticommuting.At the end, each form factor can be written as , (3.20) where the first sum runs over all the relevant topologies, FD t is the set of Feynman diagrams mapped onto a given topology t, N is a polynomial in the space-time dimension d, the Z and W masses and scalar products involving both loop and external momenta, D t,i is the i-th denominator factor of the topology t (see appendix B for their explicit definition), and n f,i ∈ {0, 1, 2}.For convenience, we linearly relate all the 9 independent kinematic invariants involving loop momenta, which are present in numerators N , to the complete set of propagators ⃗ D t , fixed by the integral topology t, such that eq.3.20 assumes the form with now It is well-known that the I t,⃗ n integrals are not independent.Indeed, they satisfy so-called Integration-By-Parts (IBP) identities [52].Typically, only after the form factors are expressed in terms of a minimal set of independent master integrals (MIs) one sees a significant reduction in the complexity of the N coefficients, as redundancies tend to be minimised.
Since in our problem there are many mass scales, satisfactorily performing a fullyanalytic reduction of the amplitude in terms of MIs is complicated.In particular, the step of expressing the various integrals in terms of MIs typically introduces complex rational functions.Although these are expected to greatly simplify when all the pieces of the amplitude are combined together, achieving such a simplification is non-trivial.Fortunately, we can simplify our calculation by setting the W and Z masses to numerical values.For numerical efficiency, we choose the m 2 W /m 2 Z ratio to be a simple rational number, while being close enough to its actual value.In particular, we choose Such a choice for the W mass differ from its actual value by few permille, and it is fully adequate for our phenomenological purposes.Setting the Z and W mass to numbers vastly reduces the parametric complexity of the integral structure.Keeping both the dimensional-regularisation parameter ϵ and the kinematics invariant {s 13 , s 23 } symbolic, we generate IBP identities by LiteRed [61,62] and solve them by the finite field arithmetic [63,64] implemented in FiniteFlow [65].We also employ the method described in ref. [66] to exploit linear relations among the reduction coefficients, which can effectively utilize their common structures to reduce the number of finite-field samples.Comparing to the traditional reconstruction strategy, this method can improve the computational efficiency by approximately an order of magnitude in our computation.
After reducing the amplitude to MIs, we evaluate the latter using the power-expansionbased differential equations method [67,68].To this end, we construct a differential equation for master integrals with respect to the kinematics invariants s 13 and s 23 .Schematically, where ⃗ I is a vector containing all the MIs and A x , are matrices.Such a writing is possible since the MIs provide a basis not only for integrals but also for all their possible derivatives.To solve the differential equations, we first compute the boundary conditions by numerically evaluating the MIs at the following regular point in s 12 scattering region: To do this, we use AMFlow [69] which implements the auxiliary mass flow method [70][71][72], which has already proven very successful in many complex calculations, see e.g.[73][74][75].With differential equations and boundary conditions in hand, the master integrals have already been fully determined.What we need to do in practice is just to use the differential equations solver in AMFlow to solve them from the boundary point to any desired phase space points.We compute the ϵ expansion of our results by fitting 10 numerical samples obtained using specific numerical values for ϵ [69].To make sure that we retain enough numerical precision even in case of large cancellations among terms, at this stage we evaluate all the MIs with about 120-digit precision.From these 10 numerical samples, we reconstruct the ϵ expansion of the result directly at the level of the form factors, without first reconstructing the analogous expansion for the MIs.After this procedure, we expect our leading pole to have very high precision of about 50 digits or better.We then lose a bit of precision for each subsequent order in ϵ, but still we estimate that our final part retains about 20-digit precision or better.
In this way, we are able to numerically evaluate all the F form factors at a given phasespace point in a robust and relatively efficient way.Contrary to a fully-analytic result, such an evaluation cannot be kinematically crossed.Therefore, we need in principle to provide 6 separate results for the independent channels In practice, we evaluate our result on a two-dimensional {p t,Z , y Z } grid that is symmetric under p 1 ↔ p 2 exchange, see sec. 5.As a consequence, we only need to compute results for the four channels

UV renormalisation and IR regularisation
The bare two-loop mixed QCD-EWK amplitude A (1,1) b has pole singularities in ϵ starting at ϵ −4 , stemming from UV and IR divergences.As we mentioned in sec.2, we renormalise the strong coupling in the MS scheme, all the EWK parameters in the on-shell scheme, and we adopt the G µ input-parameter scheme.To renormalise our result, one has to write the bare parameters in terms of the renormalised ones and to multiply by the relevant Z-factors, see eq. 2.14.To be more concrete, we expand our bare form factors as with i = 1, ..., 6 and c = L, R, see eq. 2.13.We also define renormalised form factors through the expansion with now α s = α s (µ).Similarly to eq. 2.14, bare and renormalised form factors are schematically related by At O(α s ), all the Z i factors are equal to one and the QCD renormalisation procedure amounts to writing the bare coupling α s,b in terms of the MS one as with S ϵ = (4π) ϵ e −γ E ϵ and α s = α s (µ).If we neglect contributions coming from closed fermion loops, β 0 reads The O(α) EWK renormalisation procedure is more complicated.However, since at this order neither the strong coupling nor Z g receive corrections, it is formally identical to the renormalisation of the q qZ vertex.To this order, we can write with [45] 6 In this equation, where Q is the electric charge of the external quark in units of e w .Note the γZ kinetic mixing term δ Z γZ in eq.4.7, which can be accounted for by a simple refactoring of the overall coupling.Note also that, contrarily to the MS renormalization scheme, in the on-shell scheme, the counterterms have a non-trivial ϵ expansion.Their explicit expression up to O(ϵ 0 ) can be found e.g. in ref. [76].
At O(αα s ), there is a non-trivial interplay between QCD and EWK renormalisation.However, the situation drastically simplifies if one neglects contributions coming from closed fermion loops.Indeed, in this case still Z g = 1.Also, the QCD and EWK renormalisation procedures almost decouple from each other.Precisely, only the quark wave-function renormalisation factor Z q receives O(αα s ) contributions.For left-handed quarks, it reads [26,27] The analogous result for right-handed quarks can be obtained from eq. 4.9 by substituting g q,L → g g,R and removing the W contribution.
Combining everything together, we can then write the renormalised form factors in terms of their bare counterparts as We note that since F (1,0) i,c;b contains poles, in principle we require δ (0,1) UV,c at higher orders in ϵ.In practice, this is not the case since they always decouple from physical quantities.To see how this comes about, we first need to discuss the structure of IR divergences.
The soft and collinear structure of UV-renormalised one-loop amplitudes can be immediately extracted from Catani's formula [77].In our case, we write = I (0,1) where are finite remainders and . (4.12) In these equations, s ij = (p i + p j ) 2 , for all-incoming kinematics eq.2.2.At mixed QCD-EWK order, we follow refs [26,31] and write where again F (1,1),fin i,c is finite and7 The finite remainders F (1,1),fin i,c are the main results of this work.We report them on a numerical grid in {p t,Z , y Z } in the ancillary files.For convenience, we also include the lower order results F where ellipses stand for terms that do not contain δ (0,1) is finite, higher orders in the ϵ expansion of δ (0,1) UV,c decouple from the finite remainder F (1,1),fin i,c for ϵ = 0.

Checks and final results
As we have mentioned in sec.3.3, we decided to compute the mixed QCD-EWK amplitudes numerically.In order for our result to be useful, we have performed the numerical evaluation on a dense-enough grid in {p t,Z , y Z } (see sec. 2).We focus on the boosted-Z region, where there are no thresholds.Because of this, we choose our grid to be logarithmically uniform in p t,Z and linearly uniform in y Z , with dimension 40x41.Specifically, The kinematic ranges are chosen to allow for studies both at the LHC and at future colliders.The grid 5.1 is symmetric under y Z → −y Z , which allows us to only consider a subset of the relevant partonic channels, see the discussion at the end of sec.3.3.Note that both the grid parametrization as well as the change of variables to Mandelstam invariants {p t,Z , y Z } → {s 13 , s 23 } are not rational.In order to control our numerical precision, we evaluate all the required Feynman integrals at the resulting numerical {s 13 , s 23 } grid rationalized within 8-digit agreement.The rationalised version of the grid 5.1 can be found in an ancillary file.To check the suitability of our grid for phenomenological studies, we have compared LO predictions for a) the total cross section for Z + j production with p t,Z > 200 GeV, b) the differential distributions dσ/dp t,Z , dσ/dy Z , and c) the double differential distribution dσ/dp t,Z /dy Z obtained from our grids against MCFM 6.8 [78] and found satisfactory agreement.Before presenting our final results, we describe various checks of our calculation that we have performed.First, as a byproduct of our computation we have re-computed the tree-level, one-loop QCD, and one-loop EWK amplitudes.We have benchmarked them at the level of amplitude squared against OpenLoops 2 [79].We found at least a 12digit agreement on the whole grid in all partonic channels.Second, at the mixed QCD-EWK order, we have checked that all the UV and IR poles disappear from the finite remainders F (1,1),fin i,c defined in sec 4. With our high-precision numerical evaluation of the Feynman integrals described in sec 3.3, we found ϵ-poles cancellation to ∼50 digits on the whole grid in all partonic channels.Finally, to validate our computational framework beyond the universal UV and IR structure, we have applied it to the computation of twoloop massless QCD corrections and compared our findings against results available in the literature [16,80].For this check, we picked a representative point in the uū → gZ channel and targeted a precision for the finite remainders of about 16 digits.We found perfect agreement. 8ll the tree-level, one-loop QCD, one-loop EWK, and two-loop mixed QCD-EWK finite remainders F (i,j),fin i,c evaluated on the grid 5.1 with the renormalisation scale set to can be found in a ancillary files.For illustrative purposes, here we show results for the helicity amplitudes in the uū → gZ channel at a typical point.In particular, we choose the grid point number 805, i.e.We also set the azimuthal angle ϕ of eq.2.7 to π/4.We parametrise the leptons 5 and 6 in the Z rest frame, with polar angle θ l such that cos θ l = −0.754and azimuthal angle ϕ l = π/3.The full kinematics then reads ( To present the helicity amplitudes, we define the finite part of eq.3.16 as where α s = α s (µ) and M (i,j),fin are computed from eq. 3.16 using the corresponding F (i,j),fin i,c tensors defined in sec.4. We re-absord the tree-level amplitude in the "pref ⃗ λ " terms, which are given in tab. 2 for half of the helicities.Results for the {λ 1 , λ 2 , λ 3 , −, +} helicities are obtained from the ones for {λ 1 , λ 2 , λ 3 , +, −} by exchanging 5 ↔ 6 and replacing c l,L → c l,R .We report the results for all the amplitudes in eq.5.5 in tab.3, for the scale choice eq.5.2.Finally, to illustrate the behaviour of the two-loop mixed QCD-EWK amplitudes across the whole kinematic coverage that we consider, in fig. 2 we plot col pol 2ℜ A (0,0) * A (1,1),fin , ( Table 2: Tree-level prefactors in the uū channel, see eq. 5.5 and text for details.see eq. 2.15, as a function of the Z-boson transverse momentum and rapidity, for all the independent partonic channels.The (i, j) superscript in A (i,j) indicates that only the A (i,j) contribution to eq. 2.15 is kept, see eq. 2.16.As it was the case for M fin , the "fin" superscript implies that we only consider the finite part of the amplitude, as defined in sec.4.

Conclusions and outlook
In this paper, we have performed a first important step towards the calculation of mixed QCD-EWK corrections for boosted Z production in association with one hard jet.In particular, we have computed for the first time the bosonic contributions to the two-loop mixed QCD-EWK amplitudes.Our calculation relies on a recently proposed tensor decomposition method which reduces the redundancy stemming from unphysical dimensional regularization remnants.Moreover, we exploited the modern AMFlow method [69] for highly efficient numerical evaluation of Feynman integrals.We have numerically evaluated the finite part of the two-loop amplitudes on a two-dimensional grid in {p t,Z , y Z } designed to offer good coverage for phenomenological investigations at the LHC and future colliders.We have performed extensive checks on our calculation, and expect the precision our final results to be more than sufficient for phenomenological applications.Our numerical evaluation of the tree-level, one-, and two-loop finite remainders in all the relevant partonic channels are provided in ancillary files.Our framework is easily generalisable, and we expect it to be able to cope with the case of fermionic contribution as well, which we plan to investigate in the future.This would complete the calculation of the mixed QCD-EWK two-loop amplitudes for this process and will open the way for interesting investigations.In particular, one could apply these results to phenomenological studies at the LHC.This requires devising appropriate IR subtraction schemes for mixed QCD-QED real-emission, e.g.along the lines of ref. [31].A successful completion of this programme would lead to a significant decrease on the theoretical uncertainty for important LHC analysis, like boosted Drell-Yan studies or monojet searches.It would also allow for a thorough study of the onset of the Sudakov regime, where EWK corrections are dominated by large logarithms.Insight in the transition region can provide important clues on the structure and size of subleading terms, and inform Sudakov-based approximations to mixed QCD-EWK corrections for more complicated processes.We leave these interesting avenues of exploration to the future.

A Projector coefficients
For convenience, we provide here an explicit form of the coefficients c ik required to define in eq.3.5 the projectors P i onto our Lorentz tensor basis T k .

F
section by explicitly illustrate how higher orders in ϵ in the one-loop UV EWK renormalisation counterterm δ (0,1) UV,c decouple from our finite remainder F (1,1),fin i,c .Substituting eqs 4.10 in the finite-remainder definitions eqs 4.11, 4.13, it is straightforward to find that the δ

Figure 2 :
Figure2: Absolute value of the virtual NNLO finite remainders summed over colour and polarisations, see eq. 5.6, in all partonic channels as functions of transverse momentum p t and rapidity y of the Z boson.For simplicity, in this plot we set α s = 0.118.

Table 1 :
Complexity comparison at different loop orders.

Table 3 :
Finite amplitudes in the uū channel, see eq. 5.5 and text for details.