1 Introduction

The physics program centering on electroweak (EW) precision observables receives essential inputs from measurements of W and Z bosons at the LHC. Owing to the cancellation of many systematic uncertainties, the forward–backward asymmetry \(A_{\textrm{FB}}\) in NCDY lepton-pair production is a crucial component of this program. The \(A_{\textrm{FB}}\) asymmetry is employed for determinations of the weak mixing angle \(\theta _W\) from LHC measurements at the Z-pole [1,2,3,4], complementing LEP/SLD [5] and Tevatron [6] results.

Fig. 1
figure 1

The predicted \(A_{\text {FB}}\) as a function of the invariant mass of the dilepton system in different rapidity intervals (left) and rapidity of the dilepton system in different invariant mass intervals (right) at LO in the SM

Fig. 2
figure 2

The partial derivatives of the predicted \(A_{\text {FB}}\) with respect to \(\delta g_R^{Zu}\) (upper left), \(\delta g_L^{Zu}\) (upper right), \(\delta g_R^{Zd}\) (lower left) and \(\delta g_L^{Zd}\) (lower right) couplings as a function of the invariant mass of the dilepton system in different rapidity intervals at LO

Fig. 3
figure 3

The partial derivatives of the predicted \(A_{\text {FB}}\) with respect to \(\delta g_R^{Zu}\) (upper left), \(\delta g_L^{Zu}\) (upper right), \(\delta g_R^{Zd}\) (lower left) and \(\delta g_L^{Zd}\) (lower right) couplings as a function of the rapidity of the dilepton system in different invariant mass intervals at LO

Fig. 4
figure 4

Same as in Fig. 2 for the axial and vector couplings

Fig. 5
figure 5

Same as in Fig. 3 for the axial and vector couplings

Fig. 6
figure 6

The average error in the fitted couplings as a function of the invariant mass (left) and rapidity (right) bin widths

Table 1 The binning scheme used in our analysis
Fig. 7
figure 7

The partial derivatives of the predicted \(A_{\text {FB}}\) with respect to \(\delta g_R^{Zu}\) (upper left), \(\delta g_L^{Zu}\) (upper right), \(\delta g_R^{Zd}\) (lower left) and \(\delta g_L^{Zd}\) (lower right) couplings weighted by the inverse of the statistical uncertainty as a function of the invariant mass of the dilepton system in different rapidity intervals at LO

Fig. 8
figure 8

The partial derivatives of the predicted \(A_{\text {FB}}\) with respect to \(\delta g_R^{Zu}\) (upper left), \(\delta g_L^{Zu}\) (upper right), \(\delta g_R^{Zd}\) (lower left) and \(\delta g_L^{Zd}\) (lower right) couplings weighted by the inverse of the statistical uncertainty as a function of the rapidity of the dilepton system in different invariant mass intervals at LO

Fig. 9
figure 9

Same as in Fig. 7 for the axial and vector couplings

Fig. 10
figure 10

Same as in Fig. 8 for the axial and vector couplings

Given that parton distribution functions (PDFs) constitute one of the dominant uncertainty sources in the precision EW physics program at the LHC, it is especially relevant that the \(A_{\textrm{FB}}\) asymmetry has been shown to provide us with new sensitivity to PDFs [7,8,9,10,11]. This sensitivity is currently not exploited in global PDF extractions [12,13,14,15,16], and could potentially lead to dramatic improvements in our knowledge of PDFs. This applies, in particular, in kinematic regions which are relevant for new physics searches in the multi-TeV region at the LHC, for instance in the context of Beyond-Standard-Model (BSM) heavy \(Z^\prime \) [17, 18] and \(W^\prime \) bosons [19], and photon-induced di-lepton production processes [20,21,22].

Furthermore, as in Ref. [11] the impact of the forward–backward \(A_{\textrm{FB}}\) asymmetry in the neutral current sector may be combined with that of the lepton-charge \(A_W\) asymmetry in the charged current sector. This points to strategies which are alternative to those taken in experimental analyses such as in Refs. [1, 3, 23], and aim at exploiting new measurements, capable of providing sensitivity to PDFs with low theoretical and experimental systematics while controlling correlations. Related investigations of the \(A_{\textrm{FB}}\) asymmetry in Ref. [24] focus on the behaviour induced by the NNPDF4.0 set [14]. See also the studies [25,26,27,28,29,30] based on the package ePump [31, 32].

To systematically investigate the role of the asymmetry in precision EW measurements, searches for BSM phenomena, and determinations of PDFs, a well-established framework is provided by the Standard Model Effective Field Theory (SMEFT) [33]. Details can be found in recent reviews [34, 35] and SMEFT fitting packages [36,37,38]. Recent SMEFT studies of precision electroweak observables in di-lepton channels at the LHC have been performed in Refs. [39,40,41,42,43,44] and analogous studies on the role of PDFs in BSM searches in Refs. [45, 46].

In this paper we will concentrate on \(A_{\textrm{FB}}\) asymmetry measurements in NCDY production in the region near the Z-boson mass scale. The analysis will be performed in the framework of the SMEFT Lagrangian, including operators up to dimension \(D = 6\) [33, 47],

$$\begin{aligned} \mathcal{L} = \mathcal{L}^{\mathrm{(SM)}} + {1 \over \Lambda ^2} \sum _{j=1}^{N_6} C_j^{\mathrm{(6)}} \ \mathcal{O}_j^{\mathrm{(6)}} \; , \end{aligned}$$
(1)

where the first term on the right hand side is the SM Lagrangian, consisting of operators of mass dimension \(D=4\), while the next term is the EFT contribution containing \(N_6\) operators \( \mathcal{O}_j\) of mass dimension \(D = 6\), each weighted by the dimensionless Wilson coefficient \(C_j\) divided by \(\Lambda ^2\), where \(\Lambda \) is the ultraviolet mass scale of the EFT.

In the di-lepton mass region near the Z-boson peak, four-fermion operators and dipole operators coupling fermions and vector bosons can be neglected [39, 43] in Eq. (1), and the whole effect of the \(D = 6\) SMEFT Lagrangian is a modification of the vector boson couplings to fermions. Using LEP constraints [5], corrections to Z-boson couplings to leptons can also be neglected [43]. We will thus focus on the SMEFT corrections to Z-boson couplings to u-type (including c) and d-type (including s, b) quarks, that are least constrained by LEP and have not comprehensively been studied at the LHC.

To explore these SMEFT couplings, we will extend the implementation of the \(A_{\textrm{FB}}\) asymmetry provided in Ref. [9], using the quantum chromodynamics (QCD) fit platform xFitter (formerly known as HERAFitter) [48,49,50]. As a check, we will recover the results of Ref. [9] on PDF extracted from \(A_{\textrm{FB}}\) pseudodata, and in addition we will obtain new constraints on Z-boson vector and axial couplings. We will examine the projected luminosity scenario of 3000 \(\hbox {fb}^{-1}\) for the High-Luminosity LHC (HL-LHC) [51,52,53].

The paper is organized as follows. In Sect. 2 we present the SMEFT treatment of the NCDY di-lepton production process in terms of EFT corrections to the SM couplings, and its implementation to make predictions for the \(A_{\textrm{FB}}\) asymmetry in xFitter. In Sect. 3 we describe the \(A_{\textrm{FB}}\) pseudodata generation. In Sect. 4 we carry out the main analysis within xFitter, leading to the determination of the SMEFT couplings. In Sect. 5 we give conclusions.

2 \(A_{\textrm{FB}}\) within SMEFT in xFitter

In this section we start by describing the \(D = 6\) SMEFT Lagrangian for Z-boson interactions with fermions, and introduce the SMEFT couplings for left-handed and right-handed u-quarks and d-quarks. Next we define the SMEFT vector and axial couplings, and express the forward–backward asymmetry \(A_{\textrm{FB}}\) in terms of these couplings. We discuss the extension of the xFitter implementation [9] for \(A_{\textrm{FB}}\) to the SMEFT case.

The SMEFT Lagrangian for the coupling of the Z-boson to fermions is given by

$$\begin{aligned} \mathcal{L}_Z^{\mathrm{(SMEFT)}}= & 2 M_Z \sqrt{ G_\mu \sqrt{2} } \ Z^\alpha \left\{ {\overline{q}}_L \gamma _\alpha \left( g_{L ({\textrm{SM}})}^{Zq} \right. \right. \nonumber \\ & + \left. \left. \delta g_L^{Zq} \right) q_L + {\overline{u}}_R \gamma _\alpha \left( g_{R ({\textrm{SM}})}^{Zu} + \delta g_R^{Zu} \right) u_R \right. \nonumber \\ & + \left. {\overline{d}}_R \gamma _\alpha \left( g_{R ({\textrm{SM}})}^{Zd} + \delta g_R^{Zd} \right) d_R \right. \nonumber \\ & + \left. \{ {\textrm{leptonic}} \, {\textrm{terms}} \} \right\} . \end{aligned}$$
(2)

Here \(q_L\) is the left-handed quark SU(2) doublet, while \(u_R\) and \(d_R\) are the right-handed quark SU(2) singlets. The left-handed and right-handed quark SM couplings are expressed in terms of the weak mixing angle \(\theta _W\) as follows,

$$\begin{aligned} & g_{R ({\textrm{SM}})}^{Zu} = 1/2 - 2/3 \sin ^2 \theta _W , \nonumber \\ & g_{L ({\textrm{SM}})}^{Zu} = - 2/3 \sin ^2 \theta _W , \nonumber \\ & g_{R ({\textrm{SM}})}^{Zd} = -1/2 + 1/3 \sin ^2 \theta _W , \nonumber \\ & g_{L ({\textrm{SM}})}^{Zd} = 1/3 \sin ^2 \theta _W . \end{aligned}$$
(3)

The SMEFT couplings are obtained from the SM couplings via the corrections \(\delta g\), i.e., \(g_{ ({\textrm{SMEFT}})} \equiv g_{ ({\textrm{SM}})} + \delta g\):

$$\begin{aligned} & g_L^{Zu} \equiv g_{L ({\textrm{SMEFT}})}^{Zu} = g_{L ({\textrm{SM}})}^{Zu} + \delta g_L^{Zu} , \nonumber \\ & g_R^{Zu} \equiv g_{R ({\textrm{SMEFT}})}^{Zu} = g_{R ({\textrm{SM}})}^{Zu} + \delta g_R^{Zu} , \nonumber \\ & g_L^{Zd} \equiv g_{L ({\textrm{SMEFT}})}^{Zd} = g_{L ({\textrm{SM}})}^{Zd} + \delta g_L^{Zd} , \nonumber \\ & g_R^{Zd} \equiv g_{R ({\textrm{SMEFT}})}^{Zd} = g_{R ({\textrm{SM}})}^{Zd} + \delta g_R^{Zd} . \end{aligned}$$
(4)

In our analysis, we assume \(\delta g_{R,L}^{Zd} = \delta g_{R,L}^{Zs} = \delta g_{R,L}^{Zb}\) and \(\delta g_{R,L}^{Zu} = \delta g_{R,L}^{Zc}\). The contributions from heavy c- and b-quarks are small, of the order of 10%. The vector and axial couplings of the Z-boson are defined by taking the combinations \( L \pm R\) of the left-handed and right-handed fermion couplings. So the SMEFT vector and axial couplings are given by

$$\begin{aligned} & g_{V }^{Zu} = g_{R }^{Zu} + g_{L }^{Zu} , \quad g_{A }^{Zu} = g_{R }^{Zu} - g_{L }^{Zu} , \nonumber \\ & g_{V }^{Zd} = g_{R }^{Zd} + g_{L }^{Zd} , \quad g_{A }^{Zd} = g_{R }^{Zd} - g_{L }^{Zd} . \end{aligned}$$
(5)

In order to maximize the sensitivity, we consider the DY triple-differential cross section in the di-lepton invariant mass \(M_{\ell \ell }\), di-lepton rapidity \(y_{\ell \ell }\) and angular variable \(\theta ^{*}\) between the outgoing lepton and the incoming quark in the Collins–Soper (CS) reference frame [54]. In this frame, the decay angle is measured from an axis symmetric with respect to the two incoming partons. The expression for the angle \(\theta ^{*}\) in the CS frame is given by

$$\begin{aligned} \cos \theta ^{*} = \dfrac{p_{{Z},\ell \ell }}{M_{\ell \ell }|p_{{Z},\ell \ell }|} \dfrac{p_{1}^{+}p_{2}^{-}-p_{1}^{-}p_{2}^{+}}{\sqrt{M^{2}_{\ell \ell }+p^{2}_{{T},\ell \ell }}}, \end{aligned}$$
(6)

where \(p_{i}^{\pm }=E_{i}\pm p_{{Z},i}\) and the index \(i = 1,2\) corresponds to the positive and negative charged lepton respectively. Here E and \(p_{{Z}}\) are the energy and the z-components of the leptonic four-momentum, respectively; \(p_{{Z},\ell \ell }\) is the di-lepton z-component of the momentum and \(p_{{T},\ell \ell }\) is the di-lepton transverse momentum. At leading order (LO) in QCD and EW theory, this cross section can be written as

$$\begin{aligned} \frac{d^3 \sigma }{dM_{\ell \ell }dy_{\ell \ell }d\cos \theta ^*}= & \frac{\pi \alpha ^2}{3M_{\ell \ell }s} \sum _q P_q \left[ f_q (x_1, Q^2) \right. \nonumber \\ & \times f_{\bar{q}} (x_2, Q^2) + \left. f_{\bar{q}} (x_1, Q^2) \right. \nonumber \\ & \times \left. f_q (x_2, Q^2)\right] , \end{aligned}$$
(7)

where s is the square of the centre-of-mass energy of the collision, \(x_{1,2} = M_{\ell \ell } e^{\pm y_{\ell \ell }}/\sqrt{s}\) are the momentum fractions of the initial-state partons, \(f_{q,\bar{q}} (x_i, Q^2)\) are their PDFs, \(Q^2\) is the squared factorization scale, and the factor \(P_q\) contains the propagators and couplings of the Z-boson, photon, and Z-\(\gamma \) interference,

$$\begin{aligned} P_q&= e^2_\ell e^2_q (1 + \cos ^2\theta ^*) \nonumber \\&\quad + \frac{M^2_{\ell \ell }(M^2_{\ell \ell } - M^2_Z)}{2 \sin ^2\theta _W \cos ^2\theta _W\left[ (M^2_{\ell \ell } - M^2_Z)^2 + \Gamma ^2_Z M^2_Z\right] } \nonumber \\&\quad \times (e_\ell e_q) \left[ g_V^{Z \ell } g_{V }^{Zq} (1 + \cos ^2\theta ^*) + 2 g_A^{Z \ell } g_{A }^{Z q} \cos \theta ^*\right] \nonumber \\&\quad + \frac{M^4_{\ell \ell }}{16 \sin ^4\theta _W \cos ^4\theta _W\left[ (M^2_{\ell \ell } - M^2_Z)^2 + \Gamma ^2_Z M^2_Z\right] } \nonumber \\&\quad \times \left\{ [ ( g_A^{Z \ell } )^2 + ( g_V^{Z \ell })^2 ] [ ( g_{A }^{Z q } )^2 + (g_{V }^{Z q } )^2 ] \right. \nonumber \\&\quad \times (1+\left. \cos ^2\theta ^*) + 8 g_A^{Z \ell } g_V^{Z \ell } g_{A }^{Z q} g_{V }^{Z q} \cos \theta ^* \right\} . \end{aligned}$$
(8)

Here \(M_Z\) and \(\Gamma _Z\) are the mass and the width of the Z boson, \(e_\ell \) and \(e_q\) are the lepton and quark electric charges, \(g_V^{Z \ell }=-1/2 + 2 \sin ^2 \theta _W\) and \(g_A^{Z \ell }=-1/2\) are the vector and axial couplings of leptons, and \(g_V^{Z q }\) and \(g_A^{Z q }\) are the SMEFT vector and axial couplings of quarks in Eqs. (4), (5). The first and third terms on the right hand side of Eq. (8) are the square of the s-channel diagrams with photon and Z-boson mediators respectively, while the second term is the interference between the two.

The forward–backward asymmetry \(A_{\textrm{FB}}^*\) is defined as

$$\begin{aligned} A_{\textrm{FB}}^*= & \left( d^2 \sigma / d M_{\ell \ell } d y_{\ell \ell } [\cos \theta ^*>0] \right. \nonumber \\ & - \left. d^2 \sigma / d M_{\ell \ell } d y_{\ell \ell } [\cos \theta ^*<0] \right) \nonumber \\ & / \left( d^2 \sigma / d M_{\ell \ell } d y_{\ell \ell } [\cos \theta ^*>0] \right. \nonumber \\ & + \left. d^2 \sigma / d M_{\ell \ell } d y_{\ell \ell } [\cos \theta ^*<0] \right) . \end{aligned}$$
(9)

We will consider the measurement of the \(A_{\textrm{FB}}^*\) asymmetry differentially in \(M_{\ell \ell }\) and \(y_{\ell \ell }\) according to Eqs. (7), (9).

To perform this study, we extend the implementation [9] of the \(A_{\textrm{FB}}^*\) asymmetry in the xFitter platform [48, 49] to (i) include the SMEFT couplings described above in Eqs. (4), (5), and (ii) upgrade the calculations to double-differential distributions in both invariant mass \(M_{\ell \ell }\) and rapidity \(y_{\ell \ell }\) of the di-lepton final-state system. The collider energy, acceptance cuts and bin boundaries in \(M_{\ell \ell }\) and \(y_{\ell \ell }\) are adjustable parameters in the present computation. Fiducial selections are applied to the leptons, by requiring them to have a transverse momentum \(p_T^\ell > 20 \) GeV and pseudorapidity \( | \eta _\ell | < 5 \). The mass effects of charm and bottom quarks in the matrix element are neglected, as appropriate for a high-scale process, and the calculation is performed in the \(n_f = 5\) flavour scheme. The input theoretical parameters are chosen to be the ones from the EW \(G_\mu \) scheme, which minimizes the impact of NLO EW corrections, see e.g. Ref. [55]. The explicit values for the relevant parameters in our analysis are the following: \(M_Z = 91.188\) GeV, \(\Gamma _Z = 2.441 \) GeV, \(M_W = 80.149 \) GeV, \(\alpha _{em} = 1/132.507\).

The predicted \(A_{\textrm{FB}}\) as a function of the invariant mass and rapidity of the di-lepton system at LO in the SM is shown in Fig. 1. The \(A_{\textrm{FB}}\) crosses zero around \( M_{\ell \ell } \approx M_Z \). Also, due to its definition using the longitudinal boost of the di-lepton system, it approaches zero at \(y_{\ell \ell } = 0\). For this calculation, we used the HERAPDF2.0 [16] PDF set, however, its general features do not depend on the PDF set.Footnote 1 We do not include any QED effects in our calculation since the experimental data are typically corrected for QED effects, and the uncertainties in these corrections are much smaller than the PDF uncertainties in \(A_{\text {FB}}\)  [3].

We next investigate the dependence of the predicted \(A_{\textrm{FB}}\) on the couplings. In Figs. 2 and 3 we show the numerically-calculated partial derivatives of the \(A_{\textrm{FB}}\) with respect to each coupling as a function of the invariant mass and rapidity of the dilepton system. Furthermore, in Figs. 4 and 5 these derivatives are shown with respect to the axial and vector couplings. It is instructive to see from Fig. 2 that the partial derivatives as functions of \( M_{\ell \ell }\) cross zero at values of \( M_{\ell \ell }\) which are almost independent of \(y_{\ell \ell }\). As a result, the partial derivatives as functions of \(y_{\ell \ell }\) vanish after integrating over the \( M_{\ell \ell }\) regions which contain such turnover points near their centers (e.g., \( \partial A_{\textrm{FB}} / \partial \delta g_{R }^{Zd} \approx 0\) for 85 \(< M_{\ell \ell } < 95 \) GeV). This is an important observation for experimental analyses which aim to measure \(A_{\textrm{FB}}\) in bins of \( M_{\ell \ell }\) and \(y_{\ell \ell }\): in particular, in order to retain sensitivity to the couplings, the binning scheme should be chosen carefully, preferably such that the points where the derivatives vanish are placed at the bin boundaries, rather than at their centers. Furthermore, the magnitudes of the derivatives give an idea of which phase-space regions are expected to be most sensitive to the couplings. However, one needs to take into account the expected statistical uncertainties also. Therefore, we will come back to this after introducing the pseudodata in the next section.

3 Generation of pseudodata sets

Suitable data files which mimic future measurements at the HL-LHC have been generated for the analysis. Namely, we used the expected HL-LHC luminosity, SM theoretical predictions and our assumption of 20% for the detector response to predict the number of events and statistical uncertainties for the future AFB measurement at the HL-LHC. The central values of the pseudodata points are set to the SM theoretical predictions. An important piece of information contained in the data files is the statistical precision associated to the \(A_{\text {FB}}\) experimental measurements in each bin. It is given by:

$$\begin{aligned} \Delta A_{\text {FB}} = \sqrt{\frac{1-{A_{\text {FB}}}^2}{N}}, \end{aligned}$$
(10)

where N is the expected total number of events in a specific invariant mass interval. We use the number of events with electron pairs from Z decays as predicted at LO with the acceptance cuts \( | \eta _\ell | < 5 \) and \(p_T^\ell > 20 \) GeV and introduce a further correction factor of \(20\%\) to model a realistic detector response [56]. The choice of LO accuracy for the expected number of events provides a conservative estimation of the statistical uncertainty. The higher-order QCD corrections through next-to-next-to-leading order (NNLO) for the DY process are, generally, moderate and do not distort much differential distributions, see, e.g. Ref. [57]. We have checked that the usage of NNLO QCD predictions would increase the expected number of events by factor \(1.1\text {--}1.4\) depending on the phase space region, so the statistical uncertainties does not change by more than 20% [9]. Furthermore, we have tested our approach by comparing the statistical uncertainties from the ATLAS measurement of \(A_{\text {FB}}\)  [3] with the ones produced using our pseudodata scenario, and found a reasonable agreement within a factor of two.Footnote 2

The pseudodata have been generated for the collider centre-of-mass energy of 13 TeV and integrated luminosity of 3000 \(\hbox {fb}^{-1}\), the designed integrated luminosity at the end of the HL-LHC stage [51]. To explore different proton PDF sets, several data files have been generated adopting the recent NNLO variants of the PDF sets CT18 [12], NNPDF4.0 [14], ABMP16 [15], HERAPDF2.0 [16] and MSHT20 [13] along with their respective uncertainties as provided by each fitting group.

Theoretical uncertainties arising from the choice of factorization and renormalization scales have been assessed. For this purpose, we used theoretical predictions at NLO obtained using MadGraph5_aMC@NLO [58] interfaced to APPLgrid [59] through aMCfast [60]. We have found that the impact of scale variations by the conventional factor of two in the theoretical predictions at NLO is small compared to the statistical uncertainties of the pseudodata, and thus we do not include it in our analysis. A similar study (focused on the impact of the \(A_{\text {FB}}\) on PDFs) was performed in Ref. [9]. Also, it is worth mentioning that even smaller theoretical uncertainties could be expected at NNLO in QCD. Furthermore, while it would be important to include the NLO EW effects in an analysis of experimental data aiming to obtain accurate central values, they are not expected to bring significant modification to the uncertainties.

Another important ingredient of the pseudodata is the binning scheme. As discussed in Sect. 2, bins with the turnover points of the partial derivatives of \(A_{\text {FB}}\) with respect to the couplings should be avoided, because in such bins the sensitivity to the couplings is washed out after integrating over the bin. In general, one needs as fine as possible bins in order to maximize the sensitivity to the parameters of interest. However, due to limited detector resolution too fine bins cannot be used. We have optimized the bin widths based on the precision of the fitted couplings which we extract from the pseudodata. As a figure of merit, we have used the geometrical average error of the couplings, i.e. the fourth root of the product of the resulting uncertainties on each of the four determined couplings, \(\root 4 \of {\Delta \delta g_R^{Zu}\Delta \delta g_R^{Zd}\Delta \delta g_L^{Zu}\Delta \delta g_L^{Zd}}\). In Fig. 6 we show this quantity as a function of the invariant mass and rapidity bin widths. Based on this study, we chose the bin width of 5 GeV in the invariant mass and 0.6 in the rapidity of the dilepton system, since a further reduction of the bin width does not improve the sensitivity to the couplings significantly. Namely, using the 10 GeV bin width in \(M_{ll}\), one will get about 1% larger uncertainties on the fitted couplings than using the 5 GeV bin width, which we find already sizeable. On the other hand, using a smaller bin width \(<5\) GeV provides only a marginal further improvement at permil level. These bin widths are feasible given the resolution of the existing detectors [61, 62]. The minimum expected number of events in a bin amounts to \(\sim 10{,}000\), thus the probability distribution function of the statistical uncertainty is well approximated by a normal distribution. Our binning scheme is given in Table 1.

It is illustrative to look at the magnitudes of the partial derivatives as a function of the invariant mass and rapidity of the dilepton system weighted by the inverse of the statistical uncertainty of pseudodata in the corresponding phase-space region, as shown in Figs. 7, 8, 9 and 10. This quantity is proportional to the sensitivity to the couplings which can be extracted from such a phase-space region. The largest sensitivity to the couplings is expected in the region \(55 \lesssim M_{ll} \lesssim 110\) GeV (see Figs. 7, 9), and in our analysis we have adopted a slightly wider range \(45< M_{ll} < 145\) GeV. Possible contributions from four-fermion operators in this kinematic range are expected at the \(\lesssim 1\permille \) level [41]. As a simplifying assumption, we neglected them in our analysis. A study of their possible impact is described in Appendix A.

We have limited the rapidity region \(|y_{ll}|<3.6\) assuming the extension of the detector acceptance up to pseudorapidity \(|\eta _l|<5\) (while potentially the region \(|y_{ll}|>3.6\) could provide even further improvement for the sensitivity to the couplings).

4 Results

The pseudodata sets are fitted with the four modifications to the couplings \(\delta g_L^{Zu}\), \(\delta g_R^{Zu}\), \(\delta g_L^{Zd}\), \(\delta g_R^{Zd}\) being free parameters. The fit is performed using the MINUIT [63] library by minimizing a \(\chi ^{2}\) expression:

$$\begin{aligned} \chi ^{2} = \Sigma _i\frac{(m^i-\Sigma _j\gamma ^i_jm^is_j-\mu _i)^2}{\delta _i^2}+\Sigma _js_j^2, \end{aligned}$$
(11)

which follows the one used in Ref. [16]. Here, \(\mu ^i\) is the measured value in bin i, \(\delta _i\) is its experimental uncertainty, \(m^i\) is the theoretical prediction, \(\gamma ^i_j\) is its relative uncertainty due to the PDF eigenvector j which is shifted in units of sigma by \(s_j\). The \(\gamma ^i_j\) relative uncertainty is calculated for each bin i by taking the difference between the theoretical prediction obtained using the j-th PDF eigenvectorFootnote 3 and the nominal theoretical prediction obtained using the central PDF set, and dividing this difference by the nominal prediction. This treatment of the PDF uncertainties follows the so-called profiling

Fig. 11
figure 11

Allowed regions for all pairs of corrections to the Z couplings to u- and d-type quarks obtained using different PDF sets

Fig. 12
figure 12

Same as in Fig. 11 for the axial and vector couplings

Fig. 13
figure 13

Allowed regions for all pairs of corrections to the Z couplings to u- and d-type quarks obtained using the ABMP16 and HERAPDF2.0 PDF sets with and without PDF uncertainties

Fig. 14
figure 14

Same as in Fig. 13 for the axial and vector couplings

technique [64, 65]. In this method, the PDF uncertainties are included in the \(\chi ^{2}\) using nuisance parameters \(s_j\) which are further constrained according to the tolerance criterion of the fit. The number of nuisance parameters \(s_j\) corresponds to the number of eigenvectors for each PDF set. While \(s_j\) are free parameters when minimizing the \(\chi ^{2}\), they do not change the number of degrees of freedom, because for each \(s_j\) the corresponding prior is added, represented by the last term in Eq. (11), which acts as a data point. As the tolerance criterion we use \(\Delta \chi ^{2} =1\). In this approach, one assumes that the new data are compatible with the theoretical predictions using the existing PDF set. No further theoretical uncertainties beyond the PDF uncertainties are considered when calculating the \(\chi ^{2}\). The uncertainties on the fitted parameters are obtained using the HESSE method which computes numerically the second derivatives of the \(\chi ^{2}\) with respect to the fitted parameters [48, 50]. This assumes that the dependence of the theoretical predictions on the parameters of interest is linear near the minimum of the \(\chi ^{2}\). We cross-checked the uncertainties using the MINOS algorithm which uses the profile likelihood method to compute asymmetric confidence intervals, as well as the MNCONT algorithm [63] which explicitly finds 2D contours where the \(\chi ^{2}\) is minimal, and found a good agreement with the hessian uncertainties.

The results of the fit are presented in Figs. 11, 12, 13, 14, 15, 16, 17, 18 and 19 as allowed regions for different pairs of corrections to the Z couplings to u- and d-type quarks at confidence level (CL) of 68%.Footnote 4 Note that we fit four couplings at a time, while for presentation purposes we show 2D projection plots with different pairs of couplings. The resulting uncertainties on the couplings to the d-type quarks are roughly a factor of two larger than the corresponding uncertainties for the u-type quarks. A similar impact of the \(A_{\text {FB}}\) measurements on the PDFs was found in Ref. [9]: it was shown that these measurements are most sensitive to the weighted sum \((2/3)u_v+(1/3)d_v\) of the valence u- and d-quarks, and we report the same finding in the present study. This is related to the valence quark content of the proton and also to the fact that the \(d\bar{d}\) initiated process gets more suppressed at high rapidity (where the sensitivity to the couplings is largest) in comparison to the \(u\bar{u}\) initiated ones [8]. The exact details of this effect depend on the quark PDF behaviour at high values of the partonic momentum fraction x [17, 24, 25, 46, 66]. We find a strong correlation (up to 0.95) between the different couplings, as illustrated in Figs. 11, 12, 13, 14, 15, 16, 17, 18 and 19, while the correlation coefficients between the couplings and the quark PDFs are moderate (about 0.5). Furthermore, the latter exhibit significant variability across different values of x, reflecting the complex correlation of the different parton distributions in the original PDF sets.

In Figs. 11 and 12 the allowed regions are shown as obtained using different PDF sets.Footnote 5 These results are presented using either couplings to right- and left-handed quarks or axial-vector couplings. Both the size and the shape of the allowed regions (i.e. the correlations between the fitted parameters) are similar, independent of the PDF set. To illustrate the impact of the PDF uncertainties on the results, in Figs. 13 and 14 we show the allowed regions obtained with and without including the PDF uncertainties into the fit. Furthermore, in Fig. 15 we show the average size of the uncertainties on the fitted couplings as obtained using different PDF sets. The impact of the PDF uncertainties is sizable. Namely, after including them into the fit the uncertainties on the couplings increase by a factor of \(\sim 3\), i.e. the resulting uncertainties on the couplings are dominated by the PDF uncertainties. When using different PDF sets, the size of the average uncertainty does not change by more than a factor of 1.5 indicating a reasonable consistency between the size of the PDF uncertainties for the modern PDF sets. However, it is worth mentioning also that the PDF uncertainties are constrained by the pseudodata when using the profiling technique. The behaviour of quark and antiquark densities at large Bjorken-x varies significantly between different PDF sets [7, 8, 17] strongly affecting the theoretical predictions for the \(A_{\text {FB}}\) in the high rapidity bins. However the discrepancies between different PDF sets will foreseeably reduce after the inclusion of additional experimental data covering the large-x regions [66]. The \(A_{\text {FB}}\) observable has proved particularly sensitive to quark and anti-quark densities in the large Bjorken-x region [9, 11, 17], which can be accessed through measurements of the asymmetry at high rapidities [9, 11]. In order to provide consistent results, in the future such analyses should comprise a simultaneous fit of the proton PDFs and the couplings.

In Figs. 16 and 17 we compare our results obtained for the HL-LHCFootnote 6 with the other analyses of existing data from LEP, Tevatron, HERA and LHC. The results are presented using either couplings to right- and left-handed quarks or axial and vector couplings. Namely, we compare with the analysis of the H1 Collaboration at HERA [67], the LEP+SLD combination [5], the analysis of D0 Collaboration at Tevatron [68] and the analysis of LEP, ATLAS and D0 data from Ref. [43].Footnote 7 In addition to the HL-LHC results, we present our results of analyzing all available 10 bins from the ATLAS measurement of \(A_{\text {FB}}\)  [3], while only 4 bins at the Z peak were used in the analysis of Ref. [43]. For the analysis of ATLAS data, we set the central data points to the theoretical predictions obtained at LO, while we use the data statistical and systematic uncertainties, as well as the PDF uncertainties. The analysis of ATLAS data follows the procedure used for the HL-LHC pseudodata. Such a procedure provides credible uncertainties on the Z couplings, while in order to get meaningful central values one would need to use theoretical calculations with higher-order QCD corrections. Given that the experimental uncertainties of the ATLAS data are larger than the uncertainties of the HL-LHC pseudodata, the PDF uncertainties play a moderate role in this case. The level of precision expected at the HL-LHC outperforms any existing data sets [5, 43, 67, 68, 71, 72]. Also, we note that using all the 10 data bins from the ATLAS measurement [3] provides more information on some of the couplings (e.g. \(\delta g_R^{Zd}\)) compared to the 4 bins at the Z peak together with the LEP and D0 data which were used in the analysis of Ref. [43].

Fig. 15
figure 15

The average size of the uncertainties on the fitted corrections to the Z couplings to u- and d-type quarks obtained using different PDF sets

In Figs. 18 and 19 we compare the results obtained for the HL-LHC with the results expected at the future colliders currently under discussion, LHeC [69, 73] and FCC-eh [70, 74]. For the LHeC, two electron beam energies of 50 or 60 GeV are considered, and two assumptions on the uncertainties. The results are provided for the so-called aggressive uncertainty scenario for \(E_e = 60\) GeV, and the conservative one for \(E_e = 50\) GeV (further details can be found in Ref. [69]). Furthermore, in Fig. 20 the average size of the uncertainties which can be obtained using current and future data sets are compared. A sub-percent level of precision is expected at the LHeC, FCC-eh and HL-LHC, which is one order of magnitude better than what can be obtained using existing data sets from LEP, Tevatron, HERA and LHC. More precisely, the average size of the uncertainties which are expected at the FCC-eh are a factor of 6 better than the one from our HL-LHC expectation, while the LHeC uncertainties are only 1.7–4 times smaller (depending on the scenario) than our HL-LHC results. Note that uncertainties may be also reduced at FCC-hh collider, due to increased cross section and luminosity compared to the HL-LHC. This could be studied using similar methods as those used for the HL-LHC in this paper. However, it requires a dedicated investigation of the detector acceptance and is beyond the scope of the current analysis.

Fig. 16
figure 16

Allowed regions for all pairs of corrections to the Z couplings to u- and d-type quarks obtained using HL-LHC pseudodata as well as different existing data sets

Fig. 17
figure 17

Same as in Fig. 16 for the axial and vector couplings

Fig. 18
figure 18

Allowed regions for all pairs of corrections to the Z couplings to u- and d-type quarks obtained using HL-LHC pseudodata compared to the ones for different future experiments

Fig. 19
figure 19

Same as in Fig. 18 for the axial and vector couplings

Fig. 20
figure 20

The average size of the uncertainties on the fitted corrections to the Z couplings to u- and d-type quarks for different future experiments using the linear (left) or logarithmic (right) scale

5 Conclusions

We have studied the possibility to improve constraints on the Z couplings to the u- and d-type quarks using the future measurements of \(A_{\text {FB}}\) at the HL-LHC with 3000 \(\hbox {fb}^{-1}\). We have investigated in detail the dependence of the \(A_{\text {FB}}\) on the various couplings in different regions of the invariant mass of the lepton pair, and we have shown that a wide range of \(45<M(ll)<145\) GeV is profitable to constrain the couplings. Furthermore, a measurement as function of the rapidity of the lepton pair provides a significant added value. Thus, we suggest double-differential measurements of the \(A_{\text {FB}}\) as a function of both the invariant mass and rapidity of the lepton pairs, as done e.g. in Refs. [3, 43]. Our quantitative analysis of the impact of the binning scheme on the precision of the extracted couplings suggests the choice of a specific binning scheme which ensures a substantial sensitivity to the couplings in such a measurement.

The resulting uncertainties on the couplings to the d-type quarks are found to be approximately a factor of two larger than the corresponding uncertainties for the u-type quarks. Since the \(A_{\text {FB}}\) observable is strongly sensitive to the proton PDFs, we find a significant dependence of the results on the PDF set used for such study, as was checked using the ABMP16, CT18, HERAPDF2.0, MSHT20 and NNPDF4.0 PDF sets. Preferably, in the future such analyses should comprise a simultaneous fit of the proton PDFs and the couplings.

The results were compared with the existing analyses of the LEP, HERA, Tevatron and LHC data, as well as with the results which are expected at the future colliders LHeC and FCC-eh. The uncertainties on the Z couplings to the u- and d-type quarks at the HL-LHC are expected at percent level, thus outperforming by an order of magnitude any determinations of these couplings using existing data sets. This level of precision is similar, but a little inferior to the one which is expected at the LHeC and FCC-eh.