Initial state radiation in the Herwig 7 angular-ordered parton shower

We study the simulation of initial-state radiation in angular-ordered parton showers in order to investigate how different interpretations of the ordering variable affect the logarithmic accuracy of such showers. This also enables us to implement a recoil scheme which is consistent between final-state and initial-state radiation. We present optimal values of the strong coupling and intrinsic transverse momentum to be used in each version of the parton shower, tuned using Z0-boson production at the LHC at 7 TeV. With these tuned showers, we perform a phenomenological study of the Drell-Yan process at several centre-of-mass energies.


Introduction
The complex hadronic final states of high energy particle collisions require the use of General Purpose Monte Carlo (GPMC) event generators to turn the fundamental theory into predictions of the observable final states. These event generators take a calculation of a hard scattering process and evolve it down in scale using a parton shower (PS) algorithm to a low energy scale where perturbation theory breaks down and a hadronization model is used to convert the coloured partons into hadrons, which can then also decay using a non-perturbative model. A model of Multiple Parton Interactions (MPI) is also required to describe hadron-hadron collisions.

JHEP01(2022)026
Research in this area over the last decade has focused on matching the parton shower with fixed-order hard processes, but there have also been recent attempts to improve the parton shower algorithms themselves. It was one such work [1] on the problems dipole showers 1 face even at leading-logarithmic accuracy that prompted us to make similar analyses of the Herwig 7 angular-ordered parton shower [4].
Logarithmic accuracy can be defined in terms of how well the PS reproduces the singular limits of soft-collinear matrix elements of QCD. In the limit where gluon emission is widely separated in angle and energy then for n emissions there will be 2n large logarithms, L, and the emission probability is proportional to α n s L 2n . These are the Leading Logarithmic (LL) contributions. For emissions which are collinear or occur at commensurate energies the emission probability is proportional to α n s L 2n−1 and these are the next-to-leading logarithmic contributions. Parton shower algorithms are formally only LL accurate, but by using quasi-collinear splitting functions, the transverse momentum as the scale of the strong coupling and the Catani-Marchesini-Webber (CMW) [5] scheme for the two-loop coupling, they are also able to account for LL (leading logarithmic) and NLL (next-toleading logarithmic) contributions with the exception of soft wide-angle gluon emission.
In ref. [4] we examined how the choice of recoil scheme in the shower affects the logarithmic accuracy for final-state radiation (FSR). We found that of the three schemes presented, the "p ⊥ preserving" and the "dot-product preserving" schemes both satisfied conditions to be NLL accurate while the "q 2 preserving" scheme, which was the default of Herwig 7 prior to that publication, was not NLL accurate.
In addition to FSR, there is also initial state radiation (ISR) to consider. The modelling of ISR can have a large impact on the transverse momentum distributions of colour-singlet systems produced in hadron collisions. The Drell-Yan (DY) process, in which a pair of leptons is produced in a hadron collision, is an example of a "standard-candle process", as the lepton momenta can be measured precisely to validate our understanding of QCD. Thus, it is of paramount importance to provide accurate theoretical predictions for this process. Issues arising from transverse momentum recoil due to ISR in dipole showers were addressed in ref. [6]. These were mainly related to the fact that a final state singlet, such as in Drell-Yan production, receives a non-vanishing transverse momentum contribution only from the very first shower emission. The authors proposed to allow the incoming partons to take the transverse momentum recoil, and then realign them with the beam directions at the end of the showering phase. Although this solution mitigates the problem, it is not sufficient to reach NLL accuracy, because the shower will still face the same recoil-scheme related problems present for final-state radiation [1].
In angular-ordered showers these issues are absent, both for the initial-and final-state showers, but spurious NLL terms may arise depending on the interpretation of the ordering variable in the presence of multiple emissions, as already found in ref. [4] in the context of FSR. For this reason, in this work we present the findings of our study into the effects of recoil scheme on ISR in the angular-ordered shower, focusing on the schemes which are JHEP01(2022)026 NLL accurate for FSR. This also enables us to formulate a recoil scheme which is consistent between FSR and ISR, and perform phenomenological studies of Drell-Yan processes. In this work we considered several prescriptions to improve the description of the hardest emission, and we provided dedicated tunes for all the considered options.
In section 2 we describe the implementation of the kinematic mapping for multiple emissions for different interpretations of the ordering variable. In section 3 we discuss the radiation pattern for the case of two soft emissions well separated in rapidity. The global recoil necessary to ensure full momentum conservation for the production of a colour singlet in hadron-hadron collisions is discussed in section 4 (more details and the implementation for more generic processes is discussed in appendix A). Comparisons with experimental data for Z-and W -boson production are presented in section 5. We present our conclusions in section 6.

Kinematics
We consider the emission of initial-state radiation from a parton coming from a proton. Conventionally, we consider such a parton as shower progenitor and evolve the shower backwards from the scale of the hard emission to the cutoff. We adopt the Sudakov decomposition for particles such that where the reference vectors P and n are the momentum of the parent hadron and a lightlike vector which points in the opposite direction. 2 k ⊥l is orthogonal to both P and n.
For the parton that enters the hard scattering process α l = x, where x is the fraction of the hadron's momentum that enters the hard scattering process. The parameter β l can be found by imposing that the external particles are on their mass shell.
When we consider an ISR splitting ij → i, j we denote with i the space-like child and with j the time-like one. Space-like partons are always considered massless, while the time-like parton can have a mass m j = 0. This means that when a g → bb splitting is considered, one quark will be considered as massless, while the final-state one will be treated as massive. The parent parton is on-shell and massless and it acquires an increasingly negative virtuality as it emits towards the hard scattering. The momenta obey so that the transverse momenta are defined relative to the direction of P and n. This also means that k ⊥l is space-like. 2 In the hadron-hadron frame for colour-singlet processes such as Drell-Yan, or if the colour partner is in the initial state. For processes with an outgoing colour partner, for example DIS, the direction of the colour partner in the Breit frame is used.

One emission
We consider the case of one emission. If we use the usual Sudakov decomposition the momenta are where z is the light-cone momentum fraction The coefficient β can be determined by requiring that p j is on its mass-shell The virtuality of the space-like parton is and the ordering variable can be defined as (2.7)

Multiple emissions
When multiple emissions are considered (i.e. j acquires a positive virtuality p 2 j and/or ij undergoes a further initial-state splitting), we cannot preserve simultaneously p 2 i , p j · p ij and p 2 ⊥ . It can be shown that in this case the virtuality of i is By comparing eq. (2.8) with eq. (2.6), we notice that the on-shell mass m 2 j has been replaced with the virtuality p 2 j and we have also a term proportional to p 2 ij , which is zero when only one emission is concerned. In order to complete the kinematic reconstruction we thus need to find an expression of |p ⊥ | 2 as a function ofq 2 , p 2 j and p 2 ij .

p ⊥ -preserving scheme
If we preserve the transverse momentum of each emission, we havẽ which immediately yields (2.10) JHEP01(2022)026

q 2 -preserving scheme
If we preserve the virtuality of the emission, we obtaiñ By inverting eq. (2.8) and using the above expression for p 2 i we get Since p 2 ij < 0 and p 2 j > 0, during the evolution |p ⊥ | 2 will decrease, and can eventually become negative, as already found in the context of FSR [4].

Dot-product preserving scheme
The final choice we examine is the dot-product preserving scheme: (2.14) By inverting eq. (2.8) and using p 2 In this case we see that during the evolution of the time-like parton j, |p ⊥ | is reduced, while during the evolution of the space-like parton ij it is increased. We need to ensure that there is still a physical solution, i.e. |p ⊥ | 2 ≥ 0, following subsequent initial-and final-state radiation. As further ISR can only increase |p ⊥ | 2 there is no problem, however subsequent FSR can reduce |p ⊥ | 2 and so we must ensure that further final-state radiation satisfies As shown in ref. [4], the angular-ordering condition for final-state radiation,q j < (1 − z)q, will ensure that for further time-like radiation 3 and therefore the angular-ordering condition ensures that

Summary
Remembering that we always consider the case in which m i = m ij = 0, we can summarize all the |p ⊥ | 2 expressions using where in the p ⊥ -preserving scheme, in the q 2 preserving scheme, and in the dot-product preserving scheme. It can be shown that eqs. (2.19) and (2.8) are valid also in case of ISR from a resonance. The derivation is identical to the pure ISR case, with the exception that m i = m ij = 0 and that the ordering variable for one emission is (2.23)

Double gluon emission
To discuss the impact on the logarithmic accuracy of the recoil scheme, we focus now on the case of double soft gluon emission from an incoming quark line, as shown in figure 1. Using the Sudakov decomposition described in the previous section: where the β i are determined by the requirement that k 2 i = 0 and α = x/(z 1 z 2 ). In a spacelike shower backwards evolution is used, so the partons are labelled in ascending order the further they are from the hard process.
All recoil schemes are equivalent for the last emission, therefore where we have introduced 2 = 1 − z 2 → 0 in the soft limit. The values of p 2 1 and |p ⊥1 | depend on the interpretation of the ordering variable.

p ⊥ -preserving scheme
For the p ⊥ -preserving scheme, we have Thus, from eq. (2.8) As for final-state radiation, this implies that the largest contribution to the virtuality can also come from subsequent emissions if for example 1 2 . However, since the transverse momentum of previously emitted gluons is unchanged, this scheme will reproduce the pattern of multiple independent soft emissions widely separated in angle.

q 2 -preserving scheme
For the q 2 -preserving scheme, we have thus inverting eq. (2.8) which means that |p ⊥,1 | 2 decreases and is not guaranteed to be positive. In the soft limit Clearly if 2 1 , the kinematics of the first emission is significantly modified. This, as we saw in ref. [4], causes NLL issues, so we will not to consider this scheme further.

Dot-product preserving scheme
When preserving the dot-product From eq. (2.23) we have (3.10) Thus in the soft limit It is interesting to notice that in the case of ISR, subsequent emissions tend to increase |p ⊥1 | 2 , while in the case of FSR the opposite behaviour takes place (see ref. [4]). In any case, when 1 and 2 are both small, regardless of which one is smaller, this scheme reproduces the p ⊥ -preserving one, thus it is capable of describing the matrix element for multiple independent soft gluon emissions.

Global recoil strategy for Z and Z+jet production
In a parton shower algorithm in which each parton showers independently, it is necessary to 're-assemble' the individual partons to produce the full event. Since the partons' virtualities are shifted by the showering, it is not possible to preserve all momentum components simultaneously, and some momentum must be shuffled between partons. We call the algorithm by which this is done the global recoil strategy.
In this section we summarize the global recoil strategy applied to the cases of Z and Z+jet production. More details and the generalization to an arbitrary colour structure are discussed in appendix A. We stress that in the presence of only soft and/or collinear emissions, the global recoil strategy amounts to power suppressed changes to the rapidity and the transverse momentum of partons emitted from initial state radiation and thus does not alter the discussion presented in the previous section.
Let us consider a Drell-Yan process, which at LO is described by the annihilation of a qq pair into a massive gauge boson. Each incoming parton is identified as a shower progenitor by the Herwig angular-ordered shower, and thus showered independently. The two initial partons form a colour-connected neutral system so the transverse momentum of each new emission is defined with respect to the qq direction. After the showering phase, the original incoming partons have acquired a negative virtuality and a transverse momentum. The total transverse momentum imbalance is reabsorbed by the Z boson, but there is some freedom in how to impose longitudinal momentum conservation. Three options are possible: 1. we fix the Z-boson rapidity; 2. we fix the longitudinal momentum of the Z boson; 3. we preserve the new off-shell momentum of the shower progenitor that does not contain the hardest emission.

JHEP01(2022)026
The latter option is the default behaviour, as it allows for a simpler matching with higher order matrix elements. Indeed for one emission it exactly reproduces the kinematics of the Catani-Seymour initial-initial dipole [7]. All these options ensure that, in the case of multiple soft emissions well separated in rapidity, the transverse momentum of the Z is given by the vector sum of all the emissions' transverse momenta and that the Z rapidity can receive only suppressed power corrections O p 2 ⊥ /m 2 z . In section 5 we will present our results for Z and Z+jet matched predictions. In the latter case, also the final state quark can be identified as a shower progenitor, which is colour connected to an initial parton. The global recoil strategy applied in this case is a hybrid between the one we just discussed for the production of a colour singlet, and the one developed for Deep Inelastic Scattering (DIS) processes. To be concrete, let us consider the qe − → qe − DIS process, which proceeds through a colour neutral t-channel exchange. The kinematic reconstruction and the transverse momentum are defined in the Breit frame, where the two quarks are back-to-back. Such a mapping must leave the momenta of the two electrons unchanged, thus the final-state quark and its children, i.e. the partons produced during its parton-shower evolution, need to absorb the transverse momentum imbalance due to initial state radiation. A longitudinal boost is also applied to the incoming (outgoing) parton and its children to ensure that the t-channel propagator is preserved.
For Z+jet production, the shower progenitor which leads to the hardest emission is reconstructed first together with its colour evolution partner. If they form an initial-initial dipole, the kinematic reconstruction devised for the Drell-Yan case is adopted, while if they form an initial-final dipole the one for DIS is implemented. Then one proceeds to the reconstruction of the remaining shower progenitor momentum, which is colour connected to a gluon jet which has already been reconstructed. In this case, the gluon jet will be boosted again to absorb the recoil of its second colour partner.

Drell-Yan production
To investigate the performance of the new recoil scheme for the angular-ordered shower, we compared predictions for the DY process at √ s = 7, 8 and 13 TeV with ATLAS and CMS data.
We have considered several options for the treatment of the hardest emission. We simulated leading order (LO) predictions, with matrix element corrections (MEC) to improve the description of the hardest emission and to populate the dead zone, i.e. the region of phase space not populated by the angular-ordered shower. We also simulated next-toleading order (NLO) predictions, obtained from Matchbox [8] machinery, which allows the inclusion of next-to-leading order corrections in either the MC@NLO [9] or POWHEG [10] matching schemes.
For LO+MEC predictions, by default Herwig uses leading-order parton distribution functions (PDFs). However, the usage of NLO PDFs is also possible as this introduces JHEP01(2022)026 higher-order differences beyond the level of accuracy of the calculation. 4 Thus, for illustrative purposes, in the next sections we compare LO results obtained with both LO and NLO PDFs, finding small differences as expected.
When using the POWHEG method, it is possible to separate the real emission contribution into a singular and a non-singular part, and exponentiate in the Sudakov form factor only the former contribution, while the latter is generated as a Born-like event with a higher particle multiplicity. This separation is somewhat arbitrary and in Herwig is controlled by the hard scale profile [11], which we turn off in our simulation (i.e. we exponentiate all the real corrections).
When the hardest emission generated by the POWHEG algorithm is inside the phase space region accessible to the parton shower algorithm (i.e. it is not in the dead zone), the kinematics are always reconstructed as a Z event. However, in the case that the hardest emission is in the dead zone, we consider two options. The first is to reconstruct the event as a Z+jet hard process when the hardest emission is in the dead zone. We label this behaviour as "POWHEG" in the plots below. But we also consider the option of always reconstructing the event as a shower emission from a Z event, which we label as "POWHEG (aS)". We consider the first to be more self-consistent, but the second is actually the current Herwig default. This choice has a negligible impact on the description of the small Z-boson transverse momentum, but leads to sizeable differences in the high-p ⊥ tails of distributions.

Tuning
Before comparing Herwig predictions with the experimental data, we need to tune the parameters sensitive to ISR: the intrinsic p ⊥ , i.e. the non-perturbative intrinsic transverse momentum for the partons inside the incoming hadron, and the strong coupling α s (which is given in the CMW [5] scheme). To this end, we have introduced the possibility of setting the value of α s for ISR independently from that which was previously tuned for FSR. We do not consider an independent variation of the minimum transverse momentum used as a shower cutoff, p min ⊥ , because it is strongly anti-correlated with α s . Therefore it is sufficient only to tune α s and set p min ⊥ to its FSR value of ≈ 1 GeV [4]. 5 We consider only the dotproduct and transverse-momentum preserving schemes, as the virtuality preserving scheme was found to have NLL issues in ref. [4].
The parametric dependence of the distributions is obtained from the Professor [12] program, which also finds the parameter values that minimise the χ 2 distribution for the data we are using and parameters we are fitting. To avoid being dominated only by points with a small experimental uncertainty, during the minimization procedure we set the minimum experimental uncertainty to 1%. We compare the events analysed with 4 Here we refer to the PDF employed to evaluate the matrix element and to perform the ISR evolution.
The underlying event (UE) is not included, as it has no impact on the transverse momentum of the Z boson, which is only affected by ISR, in the Herwig 7 model. For the UE, the usage of LO PDFs is recommended as the combination of leading order matrix elements and PDFs is a better approximation to the NLO result than LO matrix elements and NLO PDFs due to the different behaviours of the gluon distribution at small x. 5 For the dot-product preserving scheme p min ⊥ = 0.958 GeV, while for the p ⊥ preserving one p min ⊥ = 0.900 GeV.

JHEP01(2022)026
Rivet [13] to the ATLAS Drell-Yan Z-boson production data at 7 TeV [14][15][16]. To reduce correlations, we consider the transverse momentum of the Z boson reconstructed from muons, and the angular correlations between e + e − pairs produced from the Z decay.
Since in our simulation we use only the matrix elements for Z plus 0 or 1 jet, we are unable to perform a realistic description of the high p ⊥ spectrum, which would require the matching with higher multiplicity matrix elements as the parton shower approximation is not valid in this region. For this reason, when tuning we consider only bins in which the Z transverse momentum is smaller than 50 GeV and φ * η < 0.8. 6 The results of our tuning procedure for the several predictions described in section 5 are shown in table 1. We notice that the behaviour of the dot-product and p ⊥ preserving schemes is similar, except when the POWHEG scheme is adopted and the hardest emission inside the dead zone is treated as a shower emission (aS). For all the considered cases, the new dot-product preserving scheme always yields a better chi-squared. We also find that the value of α s for ISR is considerably larger than the one for FSR obtained in ref. [4]. In particular, at LO (+MEC) the tuned value of α s is very close to 0.1256 in the CMW scheme, which corresponds to the well known value α s = 0.118 in the MS scheme. The usage of NLO PDFs for LO predictions does not yield a significant difference with respect to the default choice of using LO PDFs. On the other hand, when adopting the MC@NLO matching the value of the strong coupling is always smaller and the χ 2 is slightly worse. The predictions obtained with the POWHEG scheme are those with the largest chi-squared, and the tuned value of α s is always larger than the expected value of 0.1256. When we always treat the hardest emission as a shower emission (aS), the dot-product preserving scheme yields the best chi-squared, with α s = 0.1255. This is in contrast with our expectations, as the treatment of emissions in the dead zone as part of the hard process is better motivated. Furthermore, as we have already said, this is the only case where we find a significant discrepancy (both in the chi-squared and in the tuned value of the strong coupling) between the p ⊥ and the dot-product preserving schemes.
It is not a surprise that the value of the coupling obtained by tuning the dot-product preserving scheme predictions is smaller than the one from the p ⊥ scheme, as subsequent emissions tend to increase the transverse momentum of previously emitted partons (see section 3.3). This behaviour is opposite to the FSR case, where instead the dot-product preserving scheme yields an α s value larger than in the p ⊥ scheme [4].
In POWHEG we find a larger value of α s since for the hardest emission, which is completely handled by Matchbox, the CMW prescription is not included in the strong coupling evaluation. We have re-run the fits with the CMW prescription implemented and obtain α s ≈ 0.125 with distributions and χ 2 that are nearly identical. 6 The variable φ * η is introduced in ref. [15] and is defined as with ∆φ being the azimuthal opening angle between the two leptons and θ * the scattering angle of the leptons with respect to the proton beam direction in the rest frame of the dilepton system.

Results
In this section we present the results of our simulations of vector boson production at the LHC and the impact of the recoil scheme and matching procedure on the accuracy of these simulations.
Z production at 7 TeV. We begin by illustrating the distributions of Z boson transverse momentum and φ * η at 7 TeV. We expect to find good agreement between all of our predictions and the data for small values of both transverse momentum and φ * η as we have tuned our parameters using the low-p ⊥ range of these distributions.
In figure 2 we show the CMS measurement [17] of the transverse momentum of the Z boson in the low p ⊥ region. We can observe that the dot-product and p ⊥ preserving schemes are almost identical (left pane), and very small differences are found including higher order corrections (right pane). Figure 3 illustrates the distribution of the φ * η parameter of the Z boson measured by the ATLAS collaboration [15]   Herwig results, obtained using the POWHEG matching prescription, are compared to CMS [17] and ATLAS [15] data at 7 TeV. φ * η all the theoretical predictions are very similar and agree well with data. We observe deviations between the recoil schemes (left panel) and among the several matching prescriptions (right panel) only for larger values of φ * η . It is interesting to notice that the LO predictions obtained using a NLO PDF (green curve in the left plot) are very similar to those obtained with the default LO PDF (red curve). We also notice that, at LO, the dot-product preserving scheme yields good agreement with data up to φ * η ∼ 0.8, which is the upper value we used in our tuning procedure, while 5% differences with respect to JHEP01(2022)026  the data arise at φ * η ∼ 0.8 in the p ⊥ -preserving scheme (blue curve in the left plot). Both the MC@NLO and POWHEG NLO predictions (right plot) give a poorer description of the data, with small discrepancies between φ * η ∼ 0.02 and φ * η ∼ 0.1, becoming very significant for φ * η > 0.2 (POWHEG) or 0.5 (MC@NLO).
In figure 4 we compare the POWHEG predictions obtained treating hardest emissions in the dead zone as part of the hard process (i.e. as a genuine "real emission"), with those obtained treating the hardest emission always as a shower emission (aS). We notice that for small values of the Z-boson p ⊥ (left panel), all the predictions are in good agreement. However, large differences arise when looking at the hard tail of the φ * η distribution (right panel), where all the Herwig predictions underestimate ATLAS data. The "real-emission" scheme, which is the best theoretically motivated, leads to the largest discrepancies with respect to the data. The "as shower" treatment yields significantly different predictions between the dot-product and p ⊥ preserving schemes. As already seen in table 1, POWHEG dot-product "as shower" (aS) results (green curve) are the ones that yield the best description of the experimental data.
W production at 7 TeV. As we have tuned using only the Z-boson transverse momentum distribution, we can use the impact of our tune on the transverse momentum of the W boson to assess the universality of our tuning procedure. In figure 5 we compare the AT-LAS 7 TeV measurement [18] with our predictions in the dot-product (left) and p ⊥ (right) preserving schemes. For W -boson transverse momentum values smaller that 50 GeV, all the predictions agree fairly well with the data, which are however plagued by large uncertainties. For larger values we see that all the theoretical predictions are systematically lower than the data, particularly those obtained using the NLO matched simulations.
Z production at 8 TeV. We now examine the Z boson distributions at 8 TeV, i.e. with a centre-of-mass energy slightly above the one we adopted for the tuning. The ATLAS measurements of ref. [19] allow us to investigate the behaviour below and above the Zboson mass peak region. Predictions in the peak region are similar to those described in the previous section and not shown here.
In figure 6 experimental data is compared to Herwig results in the dot-product preserving scheme. Above the Z-mass peak we get good agreement with the data for p ⊥ < 100 GeV and φ * η < 1. However, for masses below the peak the LO predictions overpopulate the lowp ⊥ (φ * η ) region and under populate the region with moderate p ⊥ values. This feature is present also at NLO, but it is milder and discrepancies are only of the order of a few percent. In all cases, the high-p ⊥ (φ * η ) tail is not properly described. This is not surprising, as the proper treatment of this region would require higher order matrix elements.
Z production at 13 TeV. We conclude our phenomenological study by comparing Herwig predictions to CMS measurements at 13 TeV [20]. Figure 7 shows the distribution of JHEP01(2022)026  Table 2. χ 2 results for Z Drell-Yan distributions at 13 TeV compared to CMS data [20]. For the Z-boson transverse momentum distributions, the χ 2 is computed only in the region p ⊥ < 50 GeV, while for φ * η we consider only φ * η < 0.8. The "p ⊥ , y bins" rows refer to the p ⊥ distributions broken up into 5 distributions based on y.
the Z boson transverse momentum for several rapidity ranges, with the upper-left plot showing the inclusive case. In table 2 we present the results of χ 2 calculations for the full distributions of Z boson production at 13 TeV: p ⊥ , φ * η and y (Z boson rapidity). The table also shows the total χ 2 of the p ⊥ distributions for the y sub-ranges. For this calculation we consider only p ⊥ < 50 and φ * η < 0.8.
We notice that all predictions, particularly those obtained using NLO matching prescriptions, track the data fairly well around the peak region, p Z ⊥ ≈ 10 GeV, however large differences are found both in the large and in the small p ⊥ regions. For very small p Z ⊥ values, the distribution is very sensitive to the modelling of the intrinsic p ⊥ , which we have tuned at 7 TeV. For an improved description of the data, a new tuning could be performed for each centre-of-mass energy. Alternatively, the model of [21] could be used to predict the amount of transverse momentum built up by non-perturbative smearing throughout the perturbative evolution, as a function of the partonic and hadronic centre-of-mass energies.
Furthermore, from table 2, we notice that given a matching procedure, the results obtained in the dot-product preserving scheme always yield a smaller χ 2 . The POWHEG predictions always lead to the worst descriptions of the data, apart from the rapidity distribution of the vector boson (y). This is not unexpected, since the parton shower is not allowed to change the rapidity of the Z boson, 7 thus the χ 2 value associated with the y distributions reported in table 2 depends only on the matching procedure, and not on the recoil scheme.

JHEP01(2022)026 6 Conclusions
In this paper we have generalized our previous study [4] of the logarithmic accuracy of angular-ordered parton showers to include initial-state radiation.
Similarly to the final-state case, we find that the q 2 -preserving scheme has problems in the strongly-ordered regime, with later emissions being able to modify the kinematics of an earlier soft emission and diminish its logarithmic accuracy. For this reason, we have not considered it further.
We have extended the definition of the dot-product preserving scheme introduced in ref. [4] to ISR. The implementation of this recoil option thus enables a consistent treatment of radiation which is emitted from the initial and final states. Like the p ⊥ -preserving scheme, this does not modify earlier kinematics significantly, so that multiple-gluon emission inherits the correct radiation pattern of single-gluon emission. One significant difference compared to the final-state case is that in the initial state, additional soft gluon emissions add to the transverse momentum generated by the first emission, while in the final-state case, they reduce it. This has an impact on the value of α s fitted to data, which we allow to be different for ISR and FSR.
Based on their logarithmic accuracy, we conclude that either the p ⊥ -preserving or dotproduct preserving scheme can be used for ISR. For the sake of consistency between the ISR case and the FSR case, where the dot-product preserving scheme is recommended, we recommend it also for ISR.
We have performed a dedicated tune of the non-perturbative parameters using Zproduction data at 7 TeV for each recoil scheme and matching procedure that we considered. Due to the importance of gauge boson transverse momentum distributions for the LHC physics program, we have performed a phenomenological study of the Drell-Yan process using these tuned showers. Since the description of data is better with the dotproduct preserving scheme than with the p ⊥ -preserving scheme (regardless of the matching prescription employed), we recommend the use of the dot-product preserving scheme for the additional reason that it provides a better modelling of the data.
It is impossible to consider the description of data without also considering the choice of hard process and higher order matching scheme. Of the NLO schemes, the POWHEG (aS) scheme with the dot-product preserving shower gives the best description of data. We consider this scheme to be somewhat inconsistent, because it treats real emission in the dead zone that the shower cannot populate, as if it was a parton shower emission. Nevertheless, because of its significantly better description of data, we recommend that it remain Herwig's default setting for now. The dot-product preserving reconstruction scheme will become the default also for ISR from the next Herwig release.

JHEP01(2022)026
ST/T001011/1, ST/T001038/1) and the European Union's Horizon 2020 research and innovation programme as part of the Marie Skłodowska-Curie Innovative Training Network MCnetITN3 (grant agreement no. 722104). S.F.R.'s work was also supported by the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement No. 788223, PanScales). G.B. thanks the U.K. Science and Technology Facilities Council for the award of a studentship.

A Global recoil
In this section we briefly describe the global recoil that is applied to hadron-collider events at the end of the showering phase to achieve momentum conservation. An exhaustive description of the recoil strategy of the previous Herwig++ generator, which shares many similarities with the current one (adopted since the 7.0 release), can be found in section 6.5.2 of ref. [22].
After the parton shower evolution, the space-like shower progenitors (i.e. the partons colliding in the hard process) and the time-like ones (i.e. the final-state partons arising from the hard scattering) are no longer on their mass shell, having acquired a negative or positive virtuality. Furthermore, the colliding particles have also acquired some transverse momentum that must be redistributed among the final-state progenitors and their daughter partons. We therefore need to perform some momentum reshuffling to ensure momentum conservation. How this is performed depends on whether the colour partner is an initialor final-state parton.
The details of the algorithm for final-final correlations can be found in section 6.4.2 of ref. [22], here we focus on the case where ISR is involved.

A.1 Drell-Yan: initial-initial correlations
When we consider the production of colour-singlet systems, such as electroweak gauge bosons in Drell-Yan processes, we only have an initial-initial dipole.
We use the hadronic beam momenta p ⊕ and p to define the Sudakov basis for the initial-shower algorithms. The suffix ⊕ denotes the particle incident from the +z direction, while from the −z direction. The momenta of the colliding partons q ⊕ and q can be written q ± = α ± p ± + β ± p ∓ + q ⊥± , (A.1) where q ⊥⊕ and q ⊥ are two space-like vectors orthogonal to the beam momenta. We denote by p cms the original final state momentum, i.e. the momentum of the colour singlet system. Prior to the inclusion of the shower but now the sum of the momenta of the incoming shower progenitors, q cms = q ⊕ + q is different from p cms . We can thus introduce two rescaling factors, k ± to define the shuffled momenta q q ± = α ± k ± p ± + β ± k ± p ∓ + q ⊥± , (A.3)

A.3 General case
For more complicated colour structures (like e.g. Z+jet production) we need a more general procedure. The default approach used by Herwig++ from version 2.3 and the one employed by FORTRAN HERWIG and Herwig++ versions prior to 2.3 are both detailed in ref. [22], here we want to present the new approach introduced in Herwig 7, which uses the information on the colour structure as much as possible.
The jet associated with the progenitor which leads to the hardest emission is reconstructed first. By default, its evolution partner is also reconstructed. 10 The procedure is then repeated with the next unreconstructed jet with the hardest emission. Since a gluon has two colour partners, it is shifted twice, once by the recoil from each of its partners.
For the case of DY production with matrix element corrections, the hardest emission is considered as part of the hard process (and thus treated as a shower progenitor and not as a shower emission) only when it is inside the dead zone. For POWHEG matched DY production, by default a profile function is employed to decide whether the first emission should be treated as a shower emission (and exponentiated in the Sudakov) or as part of the hard process. However, in this work, we switch off this profiling mechanism and we always exponentiate the hardest emission, but we treat it as a shower progenitor when it is inside the dead zone. For MC@NLO matched DY production, the first emitted parton is always interpreted as part of the hard process, i.e. as a shower progenitor.
In any of these cases, the global recoil is obtained by rescaling the momenta of the progenitors. However, since such rescalings are equal to 1 plus power-suppressed corrections, they do not interfere with the logarithmic structure of the result.
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.