Production of Z-bosons in the parton branching method

TL;DR

This paper combines TMD parton distributions from the Parton Branching method with NLO calculations using MCatNLO to predict Z-boson production spectra at the LHC, analyzing uncertainties and nonperturbative effects.

Contribution

It introduces a novel combination of PB TMDs with NLO MCatNLO calculations for Z-boson production, providing detailed uncertainty analysis and insights into low transverse momentum regions.

Findings

Predictions agree well with LHC measurements within uncertainties.

Renormalization and factorization scale uncertainties dominate.

Moderate impact of TMD uncertainties on predictions.

Abstract

Transverse Momentum Dependent (TMD) parton distributions obtained from the Parton Branching (PB) method are combined with next-to-leading-order (NLO) calculations of Drell-Yan (DY) production. We apply the MCatNLO method for the hard process calculation and matching with the PB TMDs. We compute predictions for the transverse momentum, rapidity and $\phi^*$ spectra of Z-bosons. We find that the theoretical uncertainties of the predictions are dominated by the renormalization and factorization scale dependence, while the impact of TMD uncertainties is moderate. The theoretical predictions agree well, within uncertainties, with measurements at the Large Hadron Collider (LHC). In particular, we study the region of lowest transverse momenta at the LHC, and comment on its sensitivity to nonperturbative TMD contributions.

Equations

$$x{\mathcal{A}}_a(x,k_t^2,{\mu }^2)$$

$${\mathcal{A}}_{0,b}(x,k_{t,0}^2,{\mu }_0^2)$$

Full text

1 Introduction

The production of $\mathup{{{Z}}}$-bosons in a Drell-Yan (DY) process [1] in $\mathup{{{p}}}\mathup{{{p}}}$ collisions is one of the most precisely measured processes at high energies at the LHC [2, 3, 4, 5]. At large transverse momenta $p_{\rm T}$ of the $\mathup{{{Z}}}$-bosons (for $p_{\rm T}>{\cal O}(m_{\mathup{{{Z}}}}))$, higher-order calculations in perturbative QCD including several additional jets are necessary to describe the measurements, while at small transverse momenta ($p_{\rm T}<{\cal O}(m_{\mathup{{{Z}}}}))$ soft-gluon perturbative resummations and nonperturbative contributions are needed [6, 7, 8, 9, 10]. An accurate description of the low-$p_{\rm T}$  spectrum of the vector bosons is important for precision measurements of the $\mathup{{{W}}}$-boson mass $m_{\mathup{{{W}}}}$, since it influences the simulation of the transverse mass or the transverse momentum of the decay leptons from which $m_{\mathup{{{W}}}}$ is extracted.

Soft-gluon resummation in the low-$p_{\rm T}$  spectrum can be achieved by analytic resummation methods [11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23] or by parton shower in multi-purpose Monte Carlo (MC) event generators [24, 25, 26, 27] matched with higher-order matrix elements [28, 29, 30, 31, 32, 33]. Nonperturbative effects can be incorporated via transverse momentum dependent parton distribution functions (TMD PDFs, or TMDs) [34] at low transverse momenta.

In Refs. [35, 36] a Parton Branching (PB) approach has been proposed which, similarly to parton shower event generators, is based on the unitarity picture [37] of parton evolution, but, unlike these event generators, uses this picture to define and evaluate TMD distribution functions. The PB method [35, 36] incorporates soft-gluon angular ordering in the parton evolution and running coupling. This enables it to achieve leading-logarithmic (LL) and next-to-leading-logarithmic (NLL) accuracy in the soft-gluon resummation, consistently with the formulation [38, 39] of coherent branching. In this respect, the results of the PB method can be compared with results based on the resummation method [10]. The main feature of the PB approach, compared to parton shower event generators, is that TMDs can be obtained and fitted to experimental data, so that the non-perturbative parameters can be fixed, and predictions can then be constructed with no further free parameters.

In this article we apply the PB approach to $\mathup{{{Z}}}$-boson DY production at the LHC. We use the TMDs obtained within the PB method and fitted [40] to inclusive deep-inelastic scattering (DIS) precision data from HERA together with a next-to-leading-order (NLO) calculation of the DY process. All parameters are kept as obtained from the fits to HERA DIS data, and no further adjustment is performed. The article is organized as follows. In Sec. 2 we describe the TMDs obtained from the PB method. In Sec. 3 we describe the DY NLO calculation. In Sec. 4 we present results for the $\mathup{{{Z}}}$-boson spectra in transverse momentum, rapidity and $\phi^{*}$ variable. We compare our results with measurements from the LHC Run-I at $\sqrt{s}=$ 8 TeV, and present predictions for $\sqrt{s}=$ 13 TeV. We give conclusions in Sec. 5.

2 Parton distributions from the PB method

The PB method allows evolution equations for collinear [41, 42, 43, 44] and TMD parton distributions to be solved numerically in an iterative procedure, by making use of the concept of resolvable and non-resolvable branchings and by applying Sudakov form factors to describe the evolution from one scale to another without resolvable branching. The PB method is described in detail in Refs. [40, 36, 35].

As discussed in detail in [40], the TMD parton density distributions are obtained from convoluting the perturbative evolution kernel ${\cal K}$ with the non-perturbative starting distribution ${\cal A}_{0,b}(x^{\prime},k_{t,0}^{2},\mu_{0}^{2})$:

[TABLE]

In general, the starting distribution ${\cal A}_{0}$ at scale $\mu_{0}$, where $\mu_{0}\approx{\cal O}(1\,{\rm{GeV}})$, can have flavor-dependent and $x$-dependent $k_{t,0}$ distributions. However, for maximal simplicity we use here a factorized form,

[TABLE]

in which the intrinsic $k_{t,0}$ distribution is given by a Gauss distribution with $\sigma^{2}=q_{0}^{2}/2$ for all parton flavors and all $x$, with a constant value $q_{0}=0.5$ .

Collinear and TMD parton distributions were obtained in [40] from fits of the parameters of the starting distribution $f_{0,b}(x,\mu_{0}^{2})$ to the inclusive-DIS precision measurements from HERA, after QCD evolution and convolution with the coefficient functions at NLO. Two different sets of parton distributions were obtained: Set 1, which corresponds at collinear level to HERAPDF2.0NLO [45], and Set 2, which differs by the choice of the scale used in the running coupling $\alpha_{\mathrm{s}}$, namely, it uses the transverse momentum (instead of the evolution scale), corresponding to the angular-ordering approach. An additional parameter $q_{cut}=1~{}$ is introduced in $\alpha_{\mathrm{s}}(\max(q^{2}_{cut},|{\bf q}_{t,}^{2}|))$ and a model dependence was estimated in Ref.  [40] with a variation around the default choice.

In Fig. 1 the collinear parton densities are shown for up-quark, strange-quark and gluon distributions at evolution scales of $\mu=10,100$ . At $\mu=10$   differences especially for the gluon are observed between Set 1 and Set 2. At scales relevant for $\mathup{{{Z}}}$ production ($\mu=100$ ) the differences between the two sets are small. The collinear parton densities are available in a format compatible with the one employed in LHAPDF [46], and can be used in the calculation of NLO processes.

In Fig. 2 we show the transverse momentum distributions for up-quark, strange-quark and gluon partons at $x=0.01$ (typical for $\mathup{{{Z}}}$ production at $\sqrt{s}=8$ ) and $\mu=100$ . The lower panels show the uncertainty bands obtained from experimental and model uncertainties.

In Fig. 3 we show the transverse momentum distributions for up-quarks and gluons at $x=0.01$ and $\mu=100$ . The band shows the uncertainty coming from changing the width of the Gauss distribution $q_{0}$ in Eq. (2) by a factor of 2 up and down in the fit as described in Ref. [40]. This variation will be used later to estimate the uncertainty of the $\mathup{{{Z}}}$-boson $p_{\rm T}$-spectrum coming from the intrinsic $k_{t}$-distribution.

We conclude this section with some comments on Eqs. (1),(2). The structure of the TMD distribution is not the same for quarks and gluons in the PB method: both the evolution kernels ${\cal K}$ and the intrinsic distributions ${\cal A}_{0}$ in Eq. (1) are in general different for different parton species. Taking a simple flavor-independent (and $x$-independent) Gaussian in Eq. (2) is not a feature of the PB method, but rather it is motivated by our finding in Ref. [40] that the precision DIS data from HERA used for the fits are not sensitive to the flavor structure of the intrinsic distribution.

A similar remark applies to nonperturbative contributions to evolution kernels. At present, the kernel ${\cal K}$ in Eq. (1) does not include any nonperturbative components. But in principle these could be introduced in the PB framework as nonperturbative contributions to the Sudakov form factors, and parameterized in terms of nonperturbative functions to be determined from fits to experimental data. Similarly to the comment above, this is not done at the moment mainly because HERA and LHC have little sensitivity to these long-distance effects.

3 PB-TMDs and NLO calculation of Drell-Yan production

In Ref. [40] leading order (LO) matrix elements for the calculation of $\mathup{{{Z}}}$ production in pp collisions at the LHC are used. A rather good description was obtained for the small-$p_{\rm T}$  region. Predictions applying PB-TMDs for the calculation of $\mathup{{{Z}}}$-boson production in proton-lead collisions have been recently reported in Ref. [47]. In the following we will describe how to use the PB TMD-parton distributions together with higher order calculations. We make use of MadGraph5_aMC@NLO (version 2.6.4, in the following labelled MC@NLO) [32] framework and apply the NLO PB parton distributions with $\alpha_{\mathrm{s}}(M_{\mathup{{{Z}}}})=0.118$ for the NLO calculations of inclusive Drell-Yan production.

In MC@NLO subtraction terms, corresponding to the parton shower used, are calculated and subtracted from the NLO cross section. Since the PB-method follows angular ordering, with the same choice of the argument in $\alpha_{\mathrm{s}}$ and in the kinematics as used in the parton shower generator Herwig++[26] and Herwig6 [48, 49], the calculation of the hard process with the subtraction terms of Herwig  is used. NLO splitting functions in the PB-method are used (consistent with NLO collinear parton densities in the MC@NLO calculation), while the subtraction terms include only LO splitting functions. However, these differences appear at an order, beyond the present next-to-leading order.

The NLO event generator MC@NLO generates events (with weights) which are stored in a format which can be read by parton shower event generators (LHE format) [50]. Instead of a parton shower event generator, we apply the PB-TMDs and modify the kinematics of the initial state partons (and as a consequence the final state partons and particles) according to the transverse momentum distributions given by the PB-TMDs. Since adding transverse momenta requires changes also in the longitudinal momenta, we require that the invariant mass of the partonic system $\hat{s}$ as well as the rapidity of the system is not changed (a method which was applied in Herwig6 and Herwig++ [26]).

The transverse momentum spectrum of $\mathup{{{Z}}}$-bosons obtained from the calculation of MC@NLO at a purely partonic level (LHE level) is shown in Fig. 4 (left). Shown are the distributions obtained with Herwig6 and Herwig++ subtraction terms. In Fig. 4 (right) the distribution is shown after transverse momenta are added according to the PB-TMDs. The differences between using calculations with the different subtractions terms are very small, as seen in the ratios in the lower panels. The PB-TMD contribute to the $p_{\rm T}$ spectrum of the $\mathup{{{Z}}}$-boson up to the scale of the hard process, not only in the non-perturbative region (since the TMDs extend to large $k_{t}$, as can be already seen from Fig. 2, 3).

In order to avoid double counting between the contribution of the real emission treated by the matrix element calculation and the contribution from the PB-TMD (or parton shower), a matching scale $\mu_{m}$ needs to be defined. This scale is determined by the NLO calculation and is transferred to the user via the parameter SCALUP (included in the LHE file).

The PB-TMDs depend, as indicated in Eq. (1), on the longitudinal momentum fraction $x$, the factorization scale $\mu$ and the transverse momentum $k_{t}$. The factorization scale $\mu$ and the longitudinal momentum fraction $x$ are used in the calculation of the hard process for collinear kinematics and a scale corresponding to $\mu=\frac{1}{2}\sum_{i}\sqrt{m^{2}_{i}+p^{2}_{t,i}}$ is chosen by default, with the sum running over all final state particles, that is, in the case of Z production, over the decay products and the final jet.

In order to allow the full phase space to be covered for the transverse momentum in the PB-TMD, the factorisation scale $\mu$ is set to the invariant mass of the hard process $\mu=\sqrt{\hat{s}}$ for the underlying Born configuration. For the real emission configuration the scale is changed to $\mu=\frac{1}{2}\sum_{i}\sqrt{m^{2}_{i}+p^{2}_{t,i}}$, as in the MC@NLO calculation. Another possible choice is to set $\mu$ to the maximum transverse momentum of the most forward and backward parton to ensure a proper matching with the angular ordering evolution of the PB-TMDs; the results are very close to the results obtained with $\mu=\frac{1}{2}\sum_{i}\sqrt{m^{2}_{i}+p^{2}_{t,i}}$. Finally, the transverse momentum is constrained to be smaller than the matching scale $\mu_{m}=$SCALUP. The calculation are performed with the Cascade3 package  [51, 52] (version 3.0.X), which allows to read LHE files and to produce output files to be analysed with Rivet [53].

The resulting transverse momentum spectrum of $\mathup{{{Z}}}$-bosons obtained from the calculation of MC@NLO after inclusion of PB-TMDs for different subtraction terms is shown in Fig. 4 (right). It is interesting to note that, for the $p_{\rm T}$ spectrum of $\mathup{{{Z}}}$-bosons calculated at inclusive level, no sensitivity to the subtraction terms is observed after the PB-TMDs are included, because the difference in Fig. 4 (right) translates into an effect of at most a percent in the total rate in Fig. 4 (left). It has been checked explicitly that a similar result is obtained when using Pythia8 [24] parton showers.

3.1 PB-method and parton showers

The basic principles of the PB-method and parton showers are very similar: both methods rely on the definition of resolvable branchings and treat every branching separately, where the kinematics can be reconstructed.

The PB-method is a method to determine the parton density, and has been applied to precision DIS measurements from HERA to determine the free parameters of the starting distributions [40]. The parton densities (collinear as well as transverse momentum dependent ones) have been obtained with NLO DGLAP splitting functions and two-loop $\alpha_{\mathrm{s}}$ with $\alpha_{\mathrm{s}}=0.118$. The evolution scale and the definition of resolvable branching is fixed. The intrinsic $k_{t}$-distribution is not really constrained from inclusive DIS measurements, and a variation of the mean of the intrinsic Gauss distribution leads to the same fits of collinear parton densities. An uncertainty is assigned coming from a variation of the mean of the Gauss distribution by a factor of two, as well as uncertainties coming from the experimental uncertainties of the data points and model uncertainties as discussed in Ref. [40]. Apart from the scale uncertainties of the hard process calculation and uncertainties of the matching scale, all parameters are fixed.

In a traditional parton shower approach (e.g. [24, 25, 26, 27]) the evolution is not constrained by the parton density, and $p_{\rm T}$- ordered or angular ordered evolution are used with any collinear parton density, together with specific choices on the definition of resolvable branchings. The free parameters are usually constrained by fits to measurements, the Monte Carlo (MC) event generator tunes, provided by the MC authors directly or by experiments (e.g. [54, 55, 56, 57]). The value of $\alpha_{\mathrm{s}}$ is usually different from the one used in the parton density, as discussed for example in Ref. [33], and special tunes are needed for LO or NLO calculations.

While the principles of PB-method and parton shower are similar, the major advantage of the PB-TMDs lies in a consistent use of parton evolution, the concept of resolvable branching and the determination of the transverse momentum during the evolution in a way which is constrained by fits to inclusive data, as the usual collinear parton densities. The central values of the predictions might be very similar in both approaches, however the uncertainties obtained in the PB-TMD come directly from the fit to inclusive DIS data, while the ones coming from a traditional parton shower are obtained separately and are not constrained by fits to parton densities. Numerically the uncertainties for the $p_{\rm T}$-spectrum of the $\mathup{{{Z}}}$-boson in the PB-method are of the order of a percent or smaller (see discussion around Fig. 5), while in a parton shower approach the uncertainties range from 50% at LO to ca. 5% for NLO splitting functions [58].

4 $\mathsfup{{{Z}}}$-boson production at the LHC

4.1 Transverse momentum and $\phi^{*}$ spectra at 8

In Fig. 5 we show the prediction for the transverse momentum spectrum of $\mathup{{{Z}}}$-bosons at $\sqrt{s}=8~{}$ obtained with a calculation using MC@NLO together with the PB-TMD [40], and compare it with the measurement by ATLAS [3]. The calculation is performed with Cascade [51, 52] (version 3.0.X) and Rivet [53]. Fig. 5 (left) shows the prediction using PB-2018-Set1 and Set2 parton distributions for both the collinear and the TMD calculations. The prediction based on PB-2018-Set1 overshoots the measurement at small transverse momentum of the $\mathup{{{Z}}}$-boson, while the prediction based on PB-2018-Set2 agrees very well with the measurement. This difference in the $p_{\rm T}$-spectrum is a direct consequence of the differences in the $k_{t}$-spectrum of the TMDs, shown in Fig. 3: in Set 2 the value of $\alpha_{s}$ at small $k_{t}$  is larger, leading to a higher probability for radiation and therefore a depletion of the distribution at small $k_{t}$. In the following we will use PB-2018-Set2 only. In Fig. 5 (left) the uncertainties of the TMD PDF are shown which come from experimental and model uncertainties (as in a fit to collinear parton densities), they are very small ($\sim 2\%$). In addition we show the effect of varying the mean of intrinsic $k_{t}$-distribution by a factor of two in Fig. 5 (left). Only in the lowest $p_{\rm T}$ bin an effect of ca. 5% can be observed. The uncertainties coming from the TMD densities and the uncertainties coming from the scale variation ($\mu_{F}$ and $\mu_{r}$ are varied by a factor of 2 up and down independently) in MC@NLO are shown in Fig. 5 (right). The uncertainties from the TMD determination are small compared to the uncertainties coming from scale variation.

The bump in the $p_{\rm T}$ distribution at $\sim 30~{}$ is an effect of the scale choice in MC@NLO, and it has been explicitly verified that a similar structure is observed when using Herwig6 or Pythia instead of the PB-TMD. The deviation of the prediction at higher transverse momenta comes from missing higher order contributions in the matrix element calculations, as only ${\cal O}(\alpha_{\mathrm{s}})$ corrections are included, and the restriction of transverse momentum of the TMD by $\mu_{m}$.

In Fig. 5 (right) the contribution from DY+1 jets (with $p_{\rm T}^{jet}>10~{}$) at NLO obtained with PB-2018 Set 2 parton distribution is shown in addition, with the band showing the uncertainty coming from the scale variation. At larger $p_{\rm T}$ the higher order contribution plays an important role and improves the description of the measurements. At this stage, we do not attempt to merge DY and DY+1 jet calculations. This will require additional studies in the use of the TMDs.

In Table 1 we show a quantitative comparison of our predictions with the measured $p_{\rm T}$ distribution. We calculate $\chi^{2}/n$ between the $n$ measurement points and the prediction taking into account the experimental uncertainties (statistical, systematic and luminosity uncertainty, added in quadrature) and the theoretical uncertainties (uncertainties from the TMD determination and scale uncertainties, as shown in Fig. 5, added in quadrature). The agreement between the measurement and the prediction using Set 2 is very good for $p_{T}^{cut}<50~{}$, although no parameters are fitted. The better $\chi^{2}$ obtained with Set 2 compared to Set 1 supports the use of transverse momentum (instead of the evolution scale), corresponding to the angular-ordering approach, as the argument of $\alpha_{s}$.

In Fig. 6 we show a comparison of the calculation with the $\phi^{*}$ distribution as measured in [3].

4.2 Predictions for $\mathsfup{{{Z}}}$-boson production at 13

In Fig. 7 we show predictions for the transverse momentum and rapidity spectra of $\mathup{{{Z}}}$-bosons at $\sqrt{s}=13~{}$ obtained, as in the previous subsection, with a calculation using MC@NLO together with the PB-TMD Set 2. The $\mathup{{{Z}}}$-bosons are selected from decay leptons with $p_{\rm T}>25~{}$ and $|\eta|<2.4$ and $|m_{ll}-m_{\mathup{{{Z}}}}|<15~{}$, following closely the selection in [59]. The uncertainties coming from the TMD parton density are shown as the blue band, while the uncertainties from a variation of the factorisation and renormalization scales are shown as the red band. In addition are shown predictions coming from a variation of the mean of the intrinsic $k_{t}$-distribution by a factor of two up and down.

In Fig. 8 the prediction for $\phi^{*}$ is shown, including TMD and scale uncertainties.

4.3 The region of small transverse momenta

The transverse momentum spectrum of the $\mathup{{{Z}}}$-bosons at small $p_{\rm T}$ is a direct measure of the intrinsic motion of the partons inside the protons as well as an important probe of the perturbative soft-gluon resummation, either in terms of TMDs or in terms of parton showers. In Fig. 9 we show predictions for the $p_{\rm T}$ spectrum with a fine binning obtained with the PB - method and standard parton showers with recent tunes: Pythia8 [24] with tune CUETP8M1 [55], Herwig++ [26] and Herwig6 [48] with default parameter settings (version 6.5.21 in MadGraph5_aMC@NLO (version 2.6.4)). All predictions make use of MC@NLO calculations of the hard process with the same collinear parton density (PB-NLO-2018-Set2) but with the appropriate subtraction terms included in the MC@NLO calculation. As expected, all calculations agree at larger $p_{\rm T}$ , while differences of up to 20 % are observed at small $p_{\rm T}<5~{}$. The prediction using Herwig6 uses parameter settings which were not tuned to recent measurements, and serves as an illustration of the sensitivity of MC tunes. With dedicated measurements in the region of $\mathup{{{Z}}}$-boson $p_{\rm T}<5-10~{}$ with fine enough binning, differences in the resummations and parton showers can be distinguished.

5 Conclusion

The PB-TMDs, which have been determined from NLO fits to inclusive DIS data, have been used to predict the $\mathup{{{Z}}}$-boson transverse momentum spectrum in \mathup{{{p}}}$$\mathup{{{p}}} collisions at the LHC for $\sqrt{s}=8$ and $13~{}$, using NLO collinear matrix element calculations for inclusive $\mathup{{{Z}}}$-production. The MC@NLO framework as implemented in MadGraph5_aMC@NLO with the Herwig6 subtraction terms is used. The PB-TMDs are used to generate transverse momenta of the incoming partons. The matching of PB-TMDs with the NLO calculation is performed by defining the factorisation scale $\mu$ for the TMD, different for the Born and real emission contributions.

The principle on which the PB method is based is similar to that of parton showers, but the difference is that in the PB method TMD densities are defined and determined from fits to experimental data, which places constraints on fixed-scale inputs to evolution.

The prediction for $\mathup{{{Z}}}$-boson production obtained with PB-TMD together with MC@NLO has been compared to measurements of ATLAS at $\sqrt{s}=8~{}$, and very good agreement with the measurement at small $p_{\rm T}$ is found using the PB-NLO-2018 Set2. The uncertainties coming from the determination of the PB-TMDs are quite small, and only in the lowest $p_{\rm T}$-region a sensitivity to the intrinsic transverse momentum spectrum is found, of the order of 2–3 %. The scale dependence of the matrix element calculation dominates the overall uncertainty.

The PB-TMD combined with MC@NLO has been used to predict the transverse momentum spectrum of $\mathup{{{Z}}}$-bosons at 13 . Measurements of the transverse momentum spectrum with a very fine binning at small $p_{\rm T}$ would allow one to separate details of the resummation procedure and the intrinsic transverse momentum distributions.

The PB-TMD have been combined for the first time with NLO collinear matrix element calculations and very good agreement with measurements for the transverse momentum and rapidity cross sections of $\mathup{{{Z}}}$-boson production is observed, without further adjusting any parameters, in contrast to what is needed in traditional parton shower approaches.

Acknowledgments. We are grateful to R. Frederix for many discussions and advice on aMCatNLO, and to A. Apyan and G. Gomez Ceballos Retuerto for very useful comments on the manuscript. FH acknowledges the support and hospitality of DESY, Hamburg while part of this work was being done. HJU thanks the Polish Science and Humboldt Foundations for the Humboldt Research fellowship during which part of this work was completed and K. Kutak for many valuable and instructive discussions. STM thanks the Humboldt Foundation for the Georg Forster research fellowship and gratefully acknowledges support from IPM.