High-energy EFT probes with fully differential Drell-Yan measurements

We study the potential of fully-differential measurements of high-energy dilepton cross-sections at the LHC to probe heavy new physics encapsulated in dimension-6 interaction operators. The assessment is performed in the seven-dimensional parameter space of operators that induce energy-growing corrections to the Standard Model partonic cross-sections at the interference level, and in the two-dimensional subspace associated with the W and Y parameters. A considerable sensitivity improvement is found relative to single-differential measurements, owing to the possibility of probing at the interference level more directions in the seven-dimensional parameter space. The reduction of parton distribution function uncertainties in the fully-differential fit is also found to play a significant role. The results are interpreted in the minimal Z' new-physics model, providing a concrete illustration of the advantages of the fully-differential analysis. We find that high-energy dilepton measurements can extend the Z' exclusion and discovery potential well beyond the reach of direct searches in a large region of the parameter space.


current-current quark-lepton operators
W and Y current-current operators appendices we discuss the definition of the kinematical variables beyond the tree-level approximation (Appendix A), we give the explicit expressions for the tree-level amplitude (Appendix B) and we report some selected results of our analysis for the LHC run 3 (Appendix C) 2 Fully-differential Drell-Yan We start our investigation of the fully-differential DY dilepton production and its sensitivity to new physics by developing a semi-analytic qualitative understanding based on the structure of the tree-level distributions. Quantitative estimates of the sensitivity are performed in Sections 2.2 and 2.3.

Tree-level distributions
Consider the neutral process qq → + − . The fully-differential cross-section is given by d 3 σ dm 2 dc * dy = τ 3 · 64 π m 4 q (1 + c * ) 2 L q (τ, y) + (1 − c * ) 2 L q (τ, −y) P q s (m ) where m = √ŝ is the dilepton invariant mass and τ =ŝ/S (with √ S the collider energy), while y is the absolute value of the rapidity (relative to the beam axis) of the dilepton system. We define c * = cos θ * as the cosine of the angle formed, in the rest frame of the dilepton pair, by the charge-minus lepton and the direction of motion of the dilepton rest frame relative to the lab frame. At tree-level, θ * as defined above is the angle between the − and the most energetic incoming parton. 3 The detailed definition of the kinematical variables beyond treelevel is reported in Appendix A.
The sum in eq. (1) spans over the light quarks q = {u, d, c, s, b}, and, for each quark species q, L q is the product of the corresponding q and q parton distribution functions (PDFs), namely L q (τ, y) = f q ( √ τ e y ; m 2 ) f q ( √ τ e −y ; m 2 ) .
3 It is essential not to define θ * with respect to a fixed beam-axis orientation. With that definition, the fullydifferential cross-section in eq. (1) would depend only on the combination (P q s + P q o ) like the single-differential cross-section in eq. (4), and all the advantages of the fully-differential analysis would be lost.
The coefficient functions P q s (P q o ) parametrize the contributions, including both SM and new physics, from the subprocesses where the chirality of the incoming quarks is the same (opposite) one of the outgoing leptons. Our target new physics operators are flavor-universal, like the SM contribution to the scattering amplitudes. Therefore the coefficient functions are the same for all the up-type and for all the down-type quarks, for a total of four independent functions P u,d s and P u,d o . In the high energy regime m m Z , and at the linear interference level in the new physics contribution, the coefficient functions read where G denotes the Wilson coefficients of the seven effective four-fermion operators defined in Table 1. In these expressions both the SM terms P u,d SM,s and P q SM,o , and the vectors V u,d s,o , are kinematics-independent numerical coefficients, reported in Appendix B. At the linear level, and up to tiny effects suppressed by m 2 Z /m 2 , the neutral DY cross-section depends on new physics only through the four linear combinations G u,d s,o = V u,d s · G of the seven Wilson coefficients. With the fully-differential analysis we can probe each of these four directions in the new physics parameters space independently, at least in line of principle.
Consider for comparison the single-differential cross-section dσ/dm 2 . By integrating eq. (1) over c * and y, we get where the parton luminosities are defined as dL q dτ = ymax −ymax dy L q (τ, y) , y max = − 1 2 log τ .
We see that the single-differential cross-section only depends on the sum of the "s" and "o" coefficient functions. At linear level, using eq. (3), it is thus only sensitive to G u s + G u o and G d s + G d o , i.e. to two combinations of the four directions in the EFT parameter space that the fully-differential analysis can probe. Actually it is not difficult to see that the single-differential analysis is not even sensitive to G u s +G u o and G d s +G d o independently, but only to the combination This is because the ratio between the up and the down quarks luminosities (that dominate over the one of the other quark flavors) is nearly constant in τ in the most sensitive energy range m ∼ 1 − 2 TeV. The ratio is approximately equal to 2 owing to the valence quarks content of the proton. The advantages of performing a fully-differential measurement can now be appreciated by analyzing the various regions in the (c * , y) kinematic space. In the kinematical regime with small center of mass rapidity (y 0), one has L q (τ, y) L q (τ, −y), and the cross-section in eq. (1) becomes proportional to (1 + c 2 * )L q (τ, 0)(P q s + P q o ). Hence this region provides sensitivity to the same combination of Wilson coefficients that can be probed through the dσ/dm 2 distribution. On the other hand, in the region with large y we have L q (τ, y) L q (τ, −y) for the (dominant) up and down quarks because the valence quarks are typically more energetic than the sea antiquarks. Therefore the cross-section is proportional to (1 + c * ) 2 P q s + (1 − c * ) 2 P q o and it is sensitive to both P q s and P q o for c * +1 and c * −1, respectively. Measuring the fully-differential distribution can also mitigate the degeneracy between the up and down quark contributions that is due, as previously discussed, to the similar shape of the parton luminosities. Indeed the dependence of L q (τ, y) on y is significantly different for the two quark species. In particular the up quark distribution is peaked at larger values of y than the one of the down quark. 4 The discussion above shows that the fully-differential cross-section measurement has the potential to disentangle the four G u,d s,o linear combinations of Wilson coefficients. This is a significant improvement relative to the single-differential measurement that is sensitive to one combination only. The quantitative assessment of this improvement is postponed to Section 2.3.
It should be stressed that our findings are based on the dependence of the cross-section on the Wilson coefficients at the linear order. At the quadratic level, all Wilson coefficients enter in the P q s,o functions with comparable coefficients (see the explicit expression for the amplitudes in Appendix B). Therefore also the combinations of parameters that do not enter or are suppressed in the linear term can be determined through their quadratic contributions. These combinations are all expected to be tested less effectively than the ones contributing to the linear terms, but with similar precision among them.
A similar analysis can be performed for the charged DY process qq → ν. In this case, however, a fully-differential measurement has a milder impact. The reason for this is twofold. First, the charged process is only affected by one operator, namely O (3) lq , so that no issue in disentangling various new physics contributions is present. Second, due to the presence of a neutrino, only two independent kinematic variables can be accessed, for instance the transverse momentum of the charged lepton p T, and its rapidity η . The new-physics contributions depend on the center of mass energy √ŝ , which is closely correlated with p T, , but has a very mild correlation with η . The additional benefit of considering both kinematical variables rather than only p T, is therefore expected to be small. However it should be taken into account that more differential information in the charged channel might help reducing the impact of PDF uncertainties in the combination with the fully-differential neutral DY measurements. Indeed some advantage of the doubly-differential measurement in charged DY will be observed in the analyses presented below.

Bounds on the Universal parameters W and Y
As a first quantitative analysis we focus on the specific set of dimension-6 operators related to the Universal parameters W and Y. As in Refs. [2,4], we define the W and Y parameters in terms of the coefficients of the four-fermion operators built from the SU(2) L and hypercharge U(1) Y currents In the DY processes, O 2W and O 2B are equivalent to lepton-quark operators. 5 The explicit correspondence is given in Table 1.
There are some crucial differences between the W and the Y parameter, which make the latter more difficult to test. While W can be probed in both the charged and the neutral DY channels, Y only affects the neutral DY process. Furthermore, if the neutral channel is analyzed by fitting only the invariant mass distribution, the single combination of Wilson coefficients that is probed at the linear level, in eq. (6), turns out to be proportional to W+0.6 Y (using Table 8). Therefore the sensitivity to Y is almost two times weaker than to W.
As we discussed above, a fully-differential analysis in the neutral channel can help to disentangle different new physics contributions. This happens also for the W and Y parameters. To illustrate this point we show in the left panel of Figure 1 6 the logarithmic derivatives of the tree-level differential cross-section with respect to W and Y, namely evaluated at the SM point W = Y = 0. For definiteness the dilepton invariant mass has been set to m = 1 TeV in the figure. The logarithmic derivative scales like m 2 as a function of the mass. As expected, in most of the kinematic space, i.e. for small rapidity and for c * 0, the cross-section dependence on W is roughly twice stronger than on Y. In particular this happens in the regions with larger cross-section, as can be seen from the plot in the right panel of Figure 1.
The behavior, however, drastically changes in the corner with c * −0.5 and y/y max 0.5. For these configurations the differential cross-section mostly depends on Y, while the sensitivity to W is small. This feature can be easily understood from the analysis we performed in the previous section. For large rapidity and c * ∼ −1 the differential cross-section is controlled by P q o , which gets contributions from subprocesses with opposite fermion chiralities. Since W corresponds to an operator with only left-handed fields, it can contribute only to the samechirality subprocesses and not to P q o . Exploiting the fully-differential distribution for the fit is thus expected to improve the determination of Y. It must be however noticed that the differential cross-section in the y ∼ y max and c * ∼ −1 corner is somewhat suppressed, and is an order of magnitude smaller than in the c * > 0 region. This means that a significant improvement in the Y determination can be obtained only when a high number of signal events are collected, so that the y ∼ y max and c * ∼ −1 region is sufficiently populated at high m . To give an idea, at the HL-LHC, out of ∼ 12000 SM events with m > 1.1 TeV, only 210 events are expected in the region with y/y max > 0.4 and c * < −0.6. We show in Figure 2 the comparison of the projected exclusion reach on the W and Y parameters obtained from a fit taking into account the fully-differential distribution or the single-differential (invariant mass or transverse momentum for neutral and charged DY, respectively) distributions. To obtain the bounds we considered the HL-LHC benchmark, with collider energy 14 TeV and L = 3 ab −1 integrated luminosity, and we assumed that the experimental measurements of the cross-section coincide with the SM predictions. 7 The fit of the charged DY process was obtained by considering a set of bins in the transverse momentum and rapidity of the charged lepton, whose boundaries are p T, : {150, 180, 225, 300, 400, 550, 750, 1000, 1300, 7000} GeV , where η max is the minimum between the acceptance cut of 2.5 and the maximal kinematically y/y max : 7 Results for the LHC run 3 benchmark are reported in Appendix C.
The cross-section predictions are obtained as in Ref. [4], at NLO in QCD combined with parton showering (based on POWHEG [27] and PYTHIA 8 [28]) and at the NLL order in the EW expansion. The effects due to the W and Y parameters (and the EW logarithms) are included through reweighting, which enables fast and accurate Monte Carlo predictions in the relatively large number of bins (a total of 234) that we consider in the fully-differential analysis. The O 2W , O 2B operators have been defined at the renormalization scale of 10 TeV.
The projected bounds take into account, following again Ref. [4], the PDF uncertainties estimated through the Hessian set PDF4LHC15 nlo 30 pdfas [29][30][31][32][33], and a 2% luminosity uncertainty. We considered an 80% reconstruction efficiency for each muon and 65% for each electron. The results do not take into account any additional experimental systematic uncertainty. This is because we expect that the size and the correlation of these uncertainties will strongly depend on the binning and it will be quite different in the fully-differential measurement and in the single-differential one. Since a quantitative estimate of the uncertainties is not available, we set them to zero for a fair comparison of the two analysis procedures. A qualitative assessment of their potential impact is presented in Section 2.4.
From Figure 2 we see that the fully-differential analysis gives a strong boost to the sensitivity of the neutral DY channel, improving in particular the sensitivity along the W + 0.6 Y = 0 line that is weakly probed by the single-differential analysis as previously discussed. The charged DY sensitivity also improves. However it should be taken into account that the single-differential analysis is performed (like in Refs. [2,4]) on the sum of the charge plus and charge minus crosssections in each p T, bin. The two charges are instead separately measured and combined in the fully-differential analysis, which is helpful to mitigate the impact of PDF uncertainties. The improvement we observe in the charged channel is partly due to this effect.
Interestingly, the improvement of the fully-differential analysis is quite significant for the combination of the neutral and charged DY channels. The 95%CL single-parameter bounds from the combined fit are given by where the numbers in brackets correspond to the single-differential fit. The constraint on W becomes nearly a factor 2 more stringent, whereas the determination of Y improves more mildly. It is interesting to notice that part of the improvement in the W determination does not come from the naive sum of the log likelihood for the neutral and charged processes, but is instead a consequence of the reduced impact of the PDF uncertainties. The PDF errors, in fact, are strongly correlated in the two channels, so that including both of them simultaneously in the fit allows one to distinguish their effects from the contributions due to new physics.

General quark-lepton interactions
We now consider the impact of the fully-differential analysis on the determination of the complete set of lepton-quark current-current operators listed in Table 1. In order to make the comparison with the single-differential analysis more straightforward, it is convenient to choose a basis in the space of Wilson coefficients which is aligned with the directions that appears in the invariant mass distribution for the neutral dilepton channel. As we discussed in Section 2.1, the dσ/dm 2 distribution depends at the linear level only on two particular combinations of parameters, . Moreover the ratio of up and down parton luminosities singles out one combination of Wilson coefficients (6) that is most effectively probed in the invariant mass distribution. We thus include in our basis the combination and the orthogonal one, which we denote by We further consider the two remaining combinations of parameters, which contribute to the fully-differential distribution at the linear level but not to the invariant mass distribution Finally, we complete the seven-dimensional basis with the G lq coefficient and two additional combinations, G ⊥ s and G ⊥ o , that are orthogonal to all the others. The explicit expressions are reported in Appendix B. The G ⊥ s and G ⊥ o coefficients contribute (at the quadratic level) to the same-chirality and opposite-chirality subprocesses, respectively. It is important to stress that the change of basis we are performing is not orthogonal. In particular this means that the G  lq (with the same value), but also gives a correlated contribution to the G ± E,O coefficients. Fully-differential measurements improve the determination of the different parameters defined above, relative to the single-differential analyses, to different extents. The G lq coefficient is mainly tested in the charged DY process, where the impact of the differential analysis is less pronounced. Its (single-operator) determination can thus improve only mildly from the combination with the neutral channel and the associated possible reduction of the impact of the PDF uncertainties. The G ⊥ s and G ⊥ o coefficients contribute only at quadratic order or through very small subleading terms in the m 2 Z /m 2 expansion, both to the single and to the fully-differential cross-section. Therefore they will be tested with lower accuracy and they will not improve significantly with the fully-differential analysis. We thus focus on the remaining four coefficients G ± E,O . Among those, G + E will not improve much, since it is already effectively probed in the invariant mass distribution. A significant improvement is instead possible for the other ones.
For a first assessment of the perspectives for progress we show in Figure 3 the logarithmic derivative of the tree-level fully-differential cross-section for the neutral DY channel with respect to the four G ± E,O parameters. In the low-rapidity region (y/y max 0.5) the cross-section is sensitive dominantly to G + E . This is not surprising since, as we saw in Section 2.1, the distribution in the low-rapidity region depends on the same combination of coefficient functions that enter in the invariant-mass distribution. The high-rapidity configurations, on the contrary, show very different sensitivity patterns to the G ± E,O coefficients. One can see, in particular, that the c * > 0 region, which has a high SM cross-section (see the right panel of Figure 1), shows a relatively large logarithmic derivative with respect to G + O . The fully-differential analysis is therefore expected to improve significantly the determination of this coefficient and to disentangle it from G + E , which gives a different dependence of the logarithmic derivative as a function of y.
On the other hand, the G − O coefficient affects the distribution mainly in the c * ∼ −1, y ∼ y max corner, in which the cross-section is rather small. For this reason we expect its determination to remain relatively poor. Finally the G − E coefficient is in an intermediate situation. The related logarithmic derivative is significantly smaller than for G + E and G + O , but nevertheless shows a distinctive pattern in the region with y/y max 0.5, which has a good cross-section. We thus expect that the fully-differential analysis could provide some improvement on its determination.
To estimate the sensitivity to G ± E,O we performed the same analysis presented in Section 2.2 for the W and Y parameters. The two-dimensional 95% CL contours for each pair of coefficients, setting the others to zero, are shown in Figure 4. 8 Different sets of bounds are compared in the plots. The solid contours correspond to the 95% CL constraints from the fully-differential analysis, whereas the dashed ones are obtained exploiting the invariant-mass distribution in the neutral channel and the transverse-momentum distribution (summed over the two charges as discussed in Section 2.2) in the charged channel. The blue shaded regions are obtained by considering the full dependence on the Wilson coefficients in the cross-section, while the orange shaded regions are found by taking into account only the linear terms. The axes of the ellipses for the fully-differential analysis at the linear level are aligned with the reference axes of each plane, owing to our judicious choice of the basis.
We also report, in Table 2, the expected sensitivity to all the seven parameters lq . We list both the single-parameter bounds obtained by setting all the others to zero and the bounds profiled over the other parameters. In the case of the fully-differential analysis we also report the results of the linearized fit.
We see from the figure and the table that the single-operator determination of G + E is only marginally modified, with a modest improvement or order 10%. This was expected, as previously discussed, since the determination of G + E from the invariant-mass distribution is already quite good. Since G + E contributes at the linear level and it is well probed, no significant difference is present between the full fit and the linearized one in this direction. Figure 4: Allowed regions at 95% CL on the six coordinate planes along the four G ± E,O coefficients. Solid contours correspond to the fully-differential analysis, while the dashed ones are obtained with the single-differential measurements. The blue shaded regions include the full dependence on the Wilson coefficients in the cross-section, while only the linear terms are retained in the orange shaded regions.
A strong improvement is instead found in the sensitivity to G + O , as anticipated. The bound from the full fit (i.e. including both the linear and quadratic dependence on the Wilson coefficient in the cross-section) improves roughly by a factor of 3. The improvement in the linearized fit is even more dramatic, since G + O does not contribute to the invariant mass distribution at the linear level up to small effects, as previously discussed. Correspondingly, an approximate flat direction is present for G + O (see for instance the middle plot on the top row of Figure 4) in the single-differential linearized contour. The fully-differential analysis is instead strongly sensitive to G + O at the linear level and the linearized and the full fit agree very well. The impact of the fully-differential analysis on the G − E and G − O parameters follows a slightly different pattern. In the full fit a mild improvement of the bounds, of order 15%, is found. The results, however, change drastically at the linearized level. In this case the fully-differential analysis is able to significantly improve the constraints on both parameters (see for instance the middle plot on the second row of Figure 4).
The profiled bounds reported in Table 2 are more difficult to interpret. They significantly differ from the single operator ones, signaling the presence of non-negligible correlations among the various parameters. We notice that for many parameters the fully-differential analysis improves the profiled bound more than the single-operator one. This pattern is particularly visible for the G lq , G + E and G − E parameters, and, to a lesser degree, for G + O and G − O . The origin of this behavior can be traced back to the reduction of flat directions in the fully-differential fit, which helps in reducing the correlations among the various Wilson coefficients.
For completeness, we report in Table 3 the bounds on the four-fermion operators in the Warsaw basis. In this basis we find that G lq is expected to be determined with much higher precision than the other parameters. Moreover its determination is only mildly affected by profiling, differently from the bounds on the other coefficients that significantly degrade in the lq , G ± E,O and G ⊥ s,o . The first three bounds correspond to single-operator fits, in which all other parameters are set to zero, while the last three are profiled over the other parameters. For each set of bounds the three columns correspond to the complete fully-differential fit, the linearized one and the single-differential measurement fit.

95%CL
single parameter profiled  profiled fit. This behavior is clearly due to the fact that G lq is tested with high precision in the charged DY channel, which is not affected by the other effective operators. The impact of a fully-differential analysis is quite large for many Warsaw operators. In particular the bounds on G lq , G qe and G lu become roughly 40% tighter, while the constraints on the other operators improve by an amount of order 10 − 20%.
As a last point, we investigate the dependence of our results on the maximal energy scale, Λ cut , of the measurements included in the fit. This gives useful indications on the measurements that contribute more to the final sensitivity and on the energy range of validity of the EFT description of new physics that is theoretically required for the bounds to apply. Following Refs. [2,4], we show on the left panel of Figure 5 how the single-operator bounds change by retaining in the fit only the bins where m < Λ cut in the neutral channel and p T, < Λ cut /2 in the charged one. We see that the bounds on the Wilson coefficients with more stringent constraints, namely G lq , G + E and G + O , saturate around Λ cut ∼ 2 TeV. In particular, removing the last bin (starting at 2600 GeV) has an extremely mild impact. The constraints on the other coefficients, on the contrary, receive sizable contributions from the events in the last bin. This behavior can be explained by recalling that the bounds on the G  in the fully-differential analysis are mainly driven by the linear interference terms in the crosssection prediction. These terms grow linearly with the partonicŝ so that the energy region with good sensitivity, which we find to be √ŝ ∼ 1 − 2 TeV, is where this growth starts being balanced by the decrease of the quark luminosity. For the other parameters, instead, the bounds are driven mostly by the square of the BSM contributions, which grow likeŝ 2 . The faster growth pushes the sensitive region to higher energies.
To appreciate better this point, we show in the right panel of Figure 5 the relative change in the bounds on the various parameters when the last bin is removed. For the fully-differential fit (darker shadowing), a variation below around 10% is observed for G  ( 30%). Finally, the G − E coefficient shows an intermediate behavior, which is explained by the fact that for its determination the linear interference terms and the quadratic terms have comparable weight. It is interesting to notice that the sensitivity to the last bin of G + O is quite lower in the fully-differential fit than in the single-differential one (displayed with lighter shadowing in the figure). This is because in the former case the bound is driven by the linear terms, while in the latter it is mainly driven by the quadratic terms. This difference can be also seen in the single-parameter bounds on G + O , which improve by roughly a factor 3 with the fully-differential analysis (see Figure 4 and Table 2).

Running Effects
Our cross-section predictions include EW corrections at the single-log accuracy, among which the ones associated with the Renormalization Group evolution of the EFT operators [4]. 9 Therefore our results depend, in line of principle, on the operator renormalization scale. This has been set to 2 TeV because the measurements at that scale dominate the sensitivity as previously shown. However the running effects are extremely small and our results do not depend on this choice in practice. This has been verified by repeating the fit in two ways. In one case we switched off completely the running, while in the other one we fixed the values of the Wilson coefficients at an energy scale E = 10 TeV. In both cases the bounds on the Wilson coefficients, both single-operator and the profiled ones, change at most by few %. Our results can thus be safely applied even to EFT operators defined at several tens of TeV. 9 The tools in Ref. [34] have been employed for the implementation.
It is important to keep in mind that only quark-lepton current-current operators are included in our calculation. Other EFT operators do not produce growing-with-energy effects in highenergy DY, therefore their contribution is very suppressed relative to the quark-lepton ones and completely negligible if their size is not anomalously large. In particular this means that operators induced by the quark-lepton ones through running are completely negligible because running is itself a small effect. The contribution to the running of the quark-lepton operators by the other ones is also negligible, for the same reason. On the other hand, one cannot firmly exclude the presence of other EFT operators with anomalously large coefficients that are not already excluded or that can not be probed with other LHC measurements. Such operators, if found to exist after a more systematic global exploration of the LHC EFT potential, should be included in the predictions.

Experimental uncertainties
This section is devoted to uncertainties that are not included in the analysis presented above, namely the presence of systematic uncertainties in the experimental cross-section measurements other than the luminosity uncertainty that was already taken into account.
Like in Ref. [4], our fits are based on the Poisson likelihood, which takes automatically into account the statistical component of the cross-section measurement errors. The systematic component of the experimental error is instead incorporated by nuisance parameters on the expected Poisson countings. Only the nuisance corresponding to a 2% luminosity uncertainty has been included in the analysis and its effect is very small as expected. The dominant experimental errors are indeed those that, unlike the luminosity, distort the shape of the differential distributions [4]. We cannot rely on any estimate of the size of these uncertainties, nor of their correlations across different analysis bins which on the other hand are expected to have a major impact on the sensitivity to new physics. In order to get a feeling of their possible impact, we adopt a crude parametrization of these effects by introducing a fully-uncorrelated 2% error in all bins, both for the fully-differential and for the single-differential fit.
Including these uncorrelated systematic errors, the combined 95% CL bound on the W and Y parameters become The numbers in parentheses refer to the single-differential fit. Comparing with the results in eq. (16), we see that the bounds from the fully-differential analysis become roughly 25% weaker.
In the case of the single-differential analysis, the bounds on Y suffer from a similar change, while the ones on W are less affected and are only 10% weaker. The advantage of a fully-differential analysis is however still evident also in these results.
In Tables 4 and 5 we give the bounds on the lepton-quark four-fermion operators including the uncorrelated systematic uncertainty. The impact of the uncertainty is relatively large on G (3) lq and G + E , whose determination becomes roughly 25% weaker. The reduction in sensitivity on the other coefficients is instead milder, at most of order 10%. Similar results are found for the operators in the Warsaw basis (Table 5). In this basis the most affected operators are G (3) lq , G eu and G ed , with a loss of sensitivity of order 25%, while the bounds on the other operators are quite stable. One can also see that the impact of the systematic uncertainty on the fully-differential fit and on the single-differential one is comparable.

Sensitivity to minimal Z models
For a concrete assessment of the benefits of the fully-differential analysis, we consider in this section a minimal BSM scenario featuring a single additional vector boson that gauges a generic linear combination of the hypercharge U(1) Y and B − L. The Lagrangian describing the new vector boson is with the current where f denote the SM fermions and H is the Higgs doublet. In the above formula g Y and g BL are free parameters, while Y (f ) and B(f ) and L(f ) are the hypercharge, the baryon and the lepton numbers of the various fermions, respectively. This model has been studied extensively in the literature and in particular in Ref. [35], where a first projection of LHC direct searches sensitivity was given, and compared with the indirect constraints from precision measurements In both panels the blue shaded region corresponds to the fully-differential fit, while the orange shaded region is obtained with the single-differential one.
(EWPT) performed at LEP and other experiments. 10 When integrated out at tree-level, the massive Z produces all the flavor-universal leptonquark operators in Table 1 lq , with Wilson coefficients that are readily computed in terms of the three free parameters g Y , g BL and M . Clearly the Wilson coefficients are quadratic polynomials in the ratios g Y /M and g BL /M , which are therefore the only two parameter combinations that can be probed by indirect searches 11 . Furthermore the indirect constraints are symmetric under an overall change of sign of the couplings (g Y , g BL ) → (−g Y , −g BL ). The 95% CL reach on the model at the HL-LHC is displayed in Figure 6 on the (g Y /M, g BL /M ) plane (left panel) and on the (g Y , g BL ) plane for a fixed mass M = 7 TeV (right panel). The bounds are obtained from the fully and single-differential analyses described in the previous section, but including in this case a 2% uncorrelated experimental uncertainty in the measurements, aiming at a more conservative result.
The advantage of the fully-differential analysis over the single-differential one is mainly in the region g BL −g Y . This region is particularly difficult to probe as it entails the cancellation of the Z coupling to the right-handed electrons as well as the suppression of the couplings to the left handed quark doublets, the right-handed up-type quarks and the left-handed lepton doublets. Therefore in this region G qe , G eu and G ed vanish and G (1) lq , G lu and G ld are suppressed (and G (3) lq is always zero). We notice in passing that the suppression of the couplings to quarks also determines a reduction of the direct production cross-section at the LHC, which makes direct searches less effective. The fully-differential analysis not only improves the sensitivity along the g BL = −g Y direction, it also mitigates the impact of the quadratic terms in the crosssection prediction. This is shown by the dashed lines in the left panel of Figure 6, reporting the results of the linearized fits. The single-differential linearized analysis possesses two very pronounced flat directions that correspond to directions in the (g Y , g BL ) plane where the G + Figure 7: HL-LHC 95% CL (1 d.o.f) exclusion reach in the mass/coupling plane for three different Z benchmark models, namely g BL = 0, g BL = −g Y and g BL = g Y . The blue shaded region can be excluded through the fully-differential di-lepton DY analysis, while the orange one can be probed with the invariant-mass fit. The green shaded region corresponds to the exclusion from direct searches. coefficient cancels. The fully-differential analysis linearized contour is instead quite close to the full fit thanks to the improved sensitivity to G + O at the linear level. In the left panel of Figure 6 we also compare our result with existing EWPT constraints, extracted from Ref. [35]. With the fully-differential analysis, the progress of the HL-LHC is of a factor around 3 in g/M in most of the directions in the (g Y /M, g BL /M ) plane, which corresponds to an improvement of one order of magnitude in the sensitivity to the Wilson coefficients that scale like (g/M ) 2 . Furthermore, notice that the EWPT bounds in the figure are based on actual experimental measurements whose central value, while compatible with the SM, disfavors the Z model. This is easily verified in the direction g BL = 0, where integrating out the Z produces only the O 2B operator with negative coefficient, that corresponds to a positive Y parameter. The central value of Y measured at LEP is instead negative (see e.g. Ref. [2]) making the EWPT exclusion on the model stronger. Our HL-LHC projections assume instead a central value at the SM point. Depending on the sign of central value that will be eventually observed the actual sensitivity to the model could be stronger or weaker than the projection.
We turn now to the comparison of our findings with the projected HL-LHC sensitivity for direct searches of the Z particle, which are most effectively performed in the dilepton final state. The exclusion on the resonant production cross-section times branching ratio is obtained from the projections in Ref. [36], slightly improved to take into account the more recent and refined results in Ref. [37]. The Z production cross-section is obtained by two MadGraph [38] simulations (at each Z mass) with the Z coupling only to up-or to down-type quarks, rescaled based on the analytical calculation of these couplings as a function of g Y and g BL . The branching ratio is also computed analytically. The results are reported in Figure 7, in the mass/coupling plane for three benchmark models (g Y = g * , g BL = 0), ( . Notice that the plot extends up to the maximal g * coupling for which, depending on the model, the width over mass ratio Γ/M of the Z is reasonably small (< 0.3) enabling a perturbative treatment. The indirect reach from our analyses, and from EWPT, is also reported in the plots. We find a substantial improvement of the mass reach for relatively large g * , up to around 30 TeV in the first and in the second benchmark model. Finally, in the right panel of Figure 6 we compare direct and indirect searches in the (g Y , g BL ) plane at a fixed mass M = 7 TeV, slightly below the threshold of around 8 TeV after which direct searches become ineffective. The direction g BL = −g Y is difficult to probe also directly, as anticipated. The sensitivity improvement of the fully-differential analysis along this direction is significant.

Discovery and characterization
High-energy measurements have the potential to discover the Z . This is shown on the left panel of Figure 8 by comparing the HL-LHC 5σ discovery reach with the current exclusion bound from EWPT in the (g Y /M, g BL /M ) plane. For M of several TeV or more, direct searches are ineffective and high-energy measurements will provide the only evidence for the existence of the Z . While "indirect", i.e. not based on the detection of a resonant peak, this evidence would be a conclusive and convincing proof of the existence of new physics thanks to the peculiar behavior (growing with energy) of the observed signal and to the possibility of getting confirmations on its nature by the study of angular distributions. The fully-differential analysis would clearly play a major role in this context, on top of course of enabling the discovery itself in a larger region of the parameter space.
We illustrate the benefits of the fully-differential analysis for the characterization of a putative signal by picking up a point (g BL /M = 0.12 TeV −1 and g Y = 0) which is discoverable at the HL-LHC, but close enough to the boundary of the discovery region to make characterization more difficult. We assume the presence of the corresponding signal in the data and we obtain the 95% CL likelihood contours on the right panel of Figure 8. A simple question related to characterization is whether we can establish that the underlying Z couples to the B −L current, rather than for instance to the hypercharge current. The figure shows that this is possible only with the fully-differential analysis.
correlations in the fully-differential fit. The augmented sensitivity at the interference level makes the fully-differential results generically more stable when the quadratic new physics terms are excluded from the predictions. This is beneficial for considerations related to the validity of the EFT, as it lowers the scale of the measurements that drive the sensitivity and reduces the impact of removing the highest-energy cross-section bins as shown in Figure 5.
The observed sensitivity improvement is due to two distinct factors. The one with the strongest impact is the extended linear-level sensitivity mentioned above and explained in Section 2.1 in details. The second, which is also quite a strong effect, is the reduction of the impact of PDF uncertainties due to their correlations across different analysis bins. These correlations are typically different from the ones of the EFT differential cross-section predictions, making harder for the PDF nuisance parameters to mimic the signal and to reduce the sensitivity. Of course the effect is quantitatively so important because the PDF are among the dominant sources of uncertainties in our fit, compatibly with the findings of Ref. [4]. Correspondingly, the benefits of the fully-differential analysis are (mildly) reduced when other sources of systematic uncertainties are assumed to be present, lowering the relative impact of PDF uncertainties in the total error budget. We have verified this fact in Section 2.4 by including a 2% systematic uncertainty uncorrelated across all bins, on top of the fully correlated luminosity uncertainty that is present (but has a totally negligible impact) in the results of Sections 2.2 and 2.3.
The dependence (at the 10% or 20% level) of our results on the assumed patterns of experimental systematic uncertainties outlines the need of detailed experimental projections for DY measurements. Experimental uncertainties that are fully correlated in all bins as in Sections 2.2 and 2.3 are definitely unrealistic. However assuming the uncertainties to be fully uncorrelated, as in Section 2.4, is equally unrealistic. We do expect correlations, especially in the fully-differential measurements, whose impact could be beneficially for the sensitivity analogously to what we have found happening for the PDF. The final HL-LHC sensitivity could thus be closer to the one in Sections 2.2 and 2.3 than to the one in Section 2.4. Furthermore, our findings are based on the statistically sub-optimal strategy of comparing cross-section measurements with EFT predictions, rather than comparing directly the EFT with the observed data. More sophisticated and unbinned strategies could be considered to further improve the sensitivity.
A significant improvement in the sensitivity is also found in the LHC run 3 projections (see Appendix C). The gain is however much milder than for the HL-LHC, mostly because the number of expected events is too low to efficiently reconstruct the full angular distributions at high energy (i.e. to populate enough all the bins required for a fully-differential analysis).
The sensitivity improvement of fully-differential measurements has a direct impact on concrete putative new physics scenarios, as we discussed in Section 3 for a simple minimal Z model. The point is that in models where the charged current O lq operator is absent, the single-differential DY analysis is mostly sensitive to a single EFT parameters combination: G + E . In the new physics model it will be generically possible to suppress G + E without particular finetuning, making the single-differential analysis loose sensitivity in a large region of the parameter space as in Figure 6. The fully-differential analysis will boost the sensitivity in that region.
The results of Section 3 also outline the effectiveness of high-p T probes on a well-established new physics benchmark that has been investigated since the beginning of the LHC program. High-p T probes by high energy DY measurements extend (see Figure 7) the projected HL-LHC exclusions well beyond the reach of direct searches in a large region of the parameter space, with sensitivity to masses from 10 to 30 TeV. Discovery is also possible up to around 20 TeV. In the event of a discovery, fully-differential measurements will play a crucial role in the characterization of the observed signal as illustrated in Figure 8. Table 6: Explicit expressions for the new physics coefficients K 0 qχ q , lχ l defined in eq. (25) for the generic quark-lepton operators in the Warsaw basis (Table 1). with a SM contribution and new physics effects encoded in K 0 qχ q , lχ l . The explicit expressions for the latter are reported in Table 6 for the quark-lepton interaction in the Warsaw basis, and in Table 7 for the W and Y parameters.
From the previous equations one can easily derive the expressions for the vectors V q s,o defined in eq. (3) as where G are the Wilson coefficients in the Warsaw basis. Analogously we defined the V ⊥ s,o vectors We recall that the G ⊥ s,o coefficients are chosen to have components only along the operators contributing to the neutral DY process and to be orthogonal to the coefficient combinations G u,d s,o . Moreover V ⊥ s and V ⊥ o contribute to same-chirality and opposite-chirality subprocesses, respectively. We report the explicit values of the V components in Table 8.
The explicit expressions for the G ± E and G ± O coefficients can be derived by substituting the above expressions in the definitions in eqs. (17), (18) and (19).

C LHC projections
In this appendix we report the projections for the LHC run 3, assuming an energy of 14 TeV and an integrated luminosity L = 300 fb −1 . The where the two sets of intervals correspond to the fully-differential and single-differential combined bounds (the latter in brackets). These bounds are obtained including only a 2% luminosity Generic quark-lepton operators W and Y  uncertainty, but no additional experimental systematic error. Comparing with the analogous HL-LHC bounds in eq. (16), we see that the all bounds are roughly a factor 2 weaker. The fullydifferential analysis provides an improvement in the sensitivity, although significantly milder than in the HL-LHC case. This pattern is not unexpected, since the benefits of a fully-differential analysis tend to be larger when the number of expected events is large enough to allow the fits to be dominated by the linear interference terms. At the LHC with L = 300 fb −1 the importance of the quadratic terms in the fits is still high, thus explaining the much milder gain in sensitivity. As we will see in the following, a similar pattern is found for the fits on the seven 4-fermion operators and on the minimal Z model. In Tables 9 and 10 we report the LHC run 3 bounds for the single parameter and the profiled fits in the basis introduced in Section 2.3 and in the Warsaw basis, respectively. Also in this case the bounds are derived assuming no uncorrelated experimental systematic uncertainty. As for the HL-LHC, the main benefit of the fully-differential fit is in the determination of the G + O parameter, which improves by roughly 40% in the single-parameter fit. Mild improvements in the sensitivity to the other parameters (of order 10 − 20%) are also found.
Finally in Figure 9 we report the projected LHC exclusion reach in the mass/coupling plane    for the three benchmark Z models considered in Section 3. The mass exclusion reach is roughly 25% weaker than the projections for the HL-LHC (compare Figure 7). The improvement due to the fully-differential analysis is clearly visible. In particular it plays a significant role in the difficult-to-test benchmark with g Y = −g BL , allowing one to cover part of the parameter space not yet excluded by the EWPT, which is not accessible with the single-differential analysis.