1 Introduction

Several of the most important measurements at LHC experiments, such as Higgs boson searches [1, 2], precision measurements of the W boson mass at LHC and Tevatron [36] or measurements of Drell–Yan (DY) processes [7] and electroweak (EW) boson pair production [5, 8] rely on a precise reconstruction of momenta for the final-state leptons [9, 10]. A substantial effort of the experimental community was devoted to optimize detector design and understand detector responses. Precision of 0.1 % (even 0.01 % for lepton directions) is of no exception. For more details, see e.g. Refs. [912].

The QED effects of the final-state radiation (FSR) play an important role in such experimental studies. FSR is included in all simulation chains and indeed should be studied together with the detector response to leptons. It can not be separated, because of infrared singularities of QED. Different approaches based on various theoretical simplifications are in use at present. At the level of the collinear approximation, expressions for higher order FSR corrections are not only well defined but are in fact process independent. In general, QED calculations are process dependent, but methods for obtaining results with \(\mathcal{O}(\alpha^{2})\) corrections and resummation of higher order effects are well established as well as techniques for evaluating theoretical errors. There is no need for introducing effective models. Situation is different if intermediate or outgoing particles, such as Z or W, are unstable, an effort like documented in [13] is then needed. In case of Z and W decays if the narrow width approximation is used, the theoretical framework for QED FSR effects is unambiguous. In case of the Z decay (in fact for hard processes mediated by s-channel Z/γ exchange) the framework in which the QED final-state radiation is separated from other contributions is defined unambiguously as well. This was explored in e + e collisions at LEP and high precision solutions were proposed.

For LEP I experiments Monte Carlo simulation programs based on exclusive exponentiation and featuring second order matrix elements were developed [14, 15]. Considerable theoretical effort was invested in the reordering of the perturbative expansion. This opened a scheme for proper resummation of vacuum polarization diagrams into terms contributing to Z width and remaining vacuum polarization corrections [16]. The definition of final-state bremsstrahlung was affected in a minimal way. The by-product of this effort was separation of the electroweak corrections into QED parts, corrections to the Z boson propagator (which have to be resummed to all orders) and remaining weak corrections which, with the proper choice of the calculation schemes, were small.

In the case of LEP II this theoretical approach had to be reconsidered, because of the new W + W pair production processes. The gauge cancellation between diagrams of electroweak bosons exchanged in s and t channels, had to be carefully respected. The resummation of the dominant contribution of the vacuum polarization had to be restricted to the case of the constant width. It was not necessary however to reopen discussion on details of scheme for final-state bremsstrahlung calculations because of the relatively small available statistics of W-pair samples [17, 18] and thus limited interest in high precision calculations for the QED bremsstrahlung.Footnote 1

It is important to stress that the QED final-state effects in processes where leptons are produced through decays of W or Z/γ can be calculated and simulated with an essentially arbitrary high precision. They form a separate class of Feynman diagrams and already developed techniques should be explored at LHC, especially as the attractiveness of such an approach was confirmed [21] in the context of W mass measurement of CDF and D0 collaborations. With the ever increasing precision, effects beyond photonic final-state bremsstrahlung have to be of course considered as part of FSR effects as well, in particular those due to emission of extra lepton pairs. Also interference effects, such as initial-final-state bremsstrahlung interference, has to be taken into account.

The interference becomes an issue for separating QED FSR from the rest of the electroweak corrections, at the level of cross section. That is why, the interference of photon emission from the initial state quarks and the final-state lepton has to be discussed carefully. It is of practical importance to verify if the suppression of the interference due to boson’s lifetime survives experimental cuts. If it is the case, then separation of effects due to the QED final-state bremsstrahlung and remaining parts of the electroweak and strong interaction corrections can be conveniently explored in the experimental studies.

For this paper we assume, that in practical applications, all other corrections than QED final-state interactions, that is remaining electroweak, and initial state hadronic interactions are expected to be measured with the properly defined observables. Thanks to this approach we will be able to achieve very competitive precision of theoretical predictions on the class of corrections directly affecting lepton momenta measurements. This represents a complementary approach to the one used in [22]. There, emphasis was put on the use of electroweak calculations together with all hadronic initial state interactions necessary for complete predictions for observables.

Our paper is organized as follows. In Sect. 2 we describe two programs PHOTOS and SANC, and their theoretical base. Section 3 is devoted to tests of the first order QED calculations. This is of importance in itself but also represent consistency checks of definitions of this part of electroweak corrections which will be considered at the amplitude level as QED final-state radiation. Definition of observables and calculation schemes used all over the paper are also given in this section. Section 4 is devoted to discussion of results for multiphoton emission and theoretical uncertainties in \(\phi^{*}_{\eta}\) measurements. Section 5 is devoted to discussion for other than photonic bremsstrahlung effects which contribute to the theoretical error of QED FSR simulated using PHOTOS or SANC. In particular, comments on emissions of pairs and interferences between photon emission from the final state and other sources are given here. Section 6, Summary, closes the paper.

Our discussion of theoretical error is limited to systematic error of QED FSR only, but it is performed in the context of full event generation. Other effects, like orientation of spin state for the intermediate W or Z bosons, affecting input for calculation of QED FSR spin amplitudes are addressed in Sect. 5.3. Precision tag for the QED FSR calculation implemented in PHOTOS or SANC is finally given for a broad class of observables: first for the photonic bremsstrahlung and then for complete FSR corrections.Footnote 2

In our work we concentrate on the applications for LHC experiments of techniques, calculations and programs developed earlier: PHOTOS [2329], SANC [3039] and KKMC [14, 15] including systematic error estimation for these new applications. Substantial and well known effort documented in these references and which form theoretical basis of our present work, will not be recalled. This would require substantial increase of our paper size with the repetition of published material. Our naming conventions originate from those papers as well. Aim of our present work is establishing physics precision of the new results of importance for LHC applications. Whenever necessary, (essentially for validating compatibility of our calculation schemes), we will address issues of technical correctness as well.

2 Description of the programs

Usually when discussing phenomenological processes at LHC in the context of the detector response one concentrates on description of the hard process and initial state QCD effects embedded e.g. in general purpose Monte Carlo generators featuring parton shower, models of underlying events, and finally QCD NLO or NNLO corrections to the hard scattering process itself are taken into account, see e.g. [40] for a review.

In this paper we concentrate on the QED FSR effects. This determines the choice of programs which will be used by us for simulations. The final-state interactions consist of QED bremsstrahlung in decays of electroweak bosons (including cases of substantial virtualities) and to some degree on the other parts of weak corrections as well. Let us recall the massive effort of years 1980–2000 for establishing definition of calculation scheme at LEP [4143] where theoretically sophisticated and numerically essential separation from electroweak corrections of the initial-state, final-state, vacuum polarization (including definition of the width) and interferences was established.

In this context discussion of theoretical errors of all parts being necessary for predictions is important, since it is not straightforward to disentangle effects of new physics and their genuine weak backgrounds. This represents however further separate work on weak corrections. In the present paper we concentrate on final-state radiation, especially on the final-state photonic one.

For the purpose of these studies, two programs featuring QED final-state radiation for LHC applications will be first described and later used. We will start with presentation of the SANC system [30], because it features complete electroweak corrections as well. Description of PHOTOS [2325, 29] implementing the QED final-state photonic bremsstrahlung only, will follow.

It might be useful to note that abbreviations LO and NLO in SANC and PHOTOS have somewhat different meaning. In the case of SANC, “LO” (the Leading Order) means just the tree-level Born cross section, while an exclusive-exponentiation-like notation is adopted in PHOTOS, where “LO” supposes exponentiation/resummation of the terms responsible for leading logarithmic terms of photonic bremsstrahlung. Full coverage of multiphoton phase-space is assured. The same concerns “NLO”: in SANC it means the one-loop approximation (Born plus \({\mathcal{O}}(\alpha)\) EW corrections), while in PHOTOS exponentiation of the \({\mathcal{O}}(\alpha)\) result is assumed.

2.1 SANC

SANC is a computer system for Support of Analytic and Numeric calculations for experiments at Colliders [30]. It can be accessed through the Internet at http://sanc.jinr.ru/ and at http://pcphsanc.cern.ch/. The SANC system is suited for calculations of one-loop QED, EW, and QCD radiative corrections (RC) to various Standard Model processes. Automatized analytic calculations in SANC provide FORM and FORTRAN modules [35], which can be used as building blocks in computer codes for particular applications.

For Drell–Yan-like processes within the SANC project there are implemented:

  • complete one-loop EW RC in the charged current [31] and neutral current [32] processes;

  • photon induced Drell–Yan processes [33];

  • higher order photonic FSR in the collinear leading logarithmic approximation;

  • higher order photonic and pair FSR in the QED leading logarithmic approximation [38, 44, 45];

  • complete next-to-leading QCD corrections [34, 36];

  • Monte Carlo integrators [37] and event generators;

  • interface to parton showers in PYTHIA and HERWIG based on the standard Les Houches Accord format.

Tuned comparisons with results of HORACE [46] and Z(W)GRADE [47] for EW RC to charged current (CC) and neutral current (NC) Drell–Yan (DY) were performed within the scope of Les Houches ’05 [48], ’07 [49] and TEV4LHC ’06 [50] workshops. A good agreement achieved in these comparisons confirms correctness of the implementation of the complete one-loop EW corrections in all these programs.

An important feature of the SANC approach is the possibility to control and directly access different contributions to the observables being under consideration. In particular, SANC code allows to separate effects due to the final-state radiation, the interference of initial and final radiation, the so called pure weak contributions, etc.

Separation of the FSR contribution in the case of neutral current DY processes is straightforward, it naturally appears at the level of Feynman diagrams. But the corresponding separation in the case of the charged current DY processes is not so trivial. In general it is even not gauge invariant. For the sake of tuned comparison with PHOTOS, a special prescription for this separationFootnote 3 was introduced into SANC.

Let us consider a formal separation of the pure weak (PW) and QED contributions δ PW and δ QED to the total \(W^{+} \to u + \bar{d} \) decay width

$$ \varGamma_W^{\mathrm{PW}+\mathrm{QED}} = \varGamma^{\mathrm{LO}}\bigl(\delta^{\mathrm{PW}} + \delta^{\mathrm{QED}}\bigr). $$
(1)

This process at the O(α) is described by 6 QED-like diagrams with virtual photon line and 3 other ones with real photon emission which together lead to the formula

$$\begin{aligned} \delta^{\mathrm{QED}} =& \frac{\alpha}{\pi} \biggl[ Q_W^2 \biggl( \frac{11}{6} - \frac{\pi^2}{3} \biggr) \\ &{}+ \bigl(Q_u^2 + Q_d^2 \bigr) \biggl( \frac{11}{8} - \frac{3}{4}\log{\frac{M_W^2}{\mu^2_{\mathrm{PW}}}} \biggr) \biggr], \end{aligned}$$

where parameter μ PW is the ’t Hooft scale introduced for separation of QED and PW contributions. In order to separate the FSR QED contribution, we choose \(\mu_{\mathrm{PW}} = M_{W} \exp(-\frac{11}{12}) \). This value of the ’t Hooft scale makes the total QED contribution to the W boson decay being equal to zero. This is in agreement with the corresponding treatment in PHOTOS, where by construction the effect of FSR to a process does not change the normalization of the cross section.

2.2 PHOTOS

Already in the era of data analysis of LEP experiments simulation of bremsstrahlung in decays of resonances and particles required specialized tools. In parallel, two programs oriented toward highest possible overall precision for the whole processes in e + e collisions such as KKMC [14] or KORALZ [51], programs dealing with decays only, gradually became of a broad use. The PHOTOS Monte Carlo was one of such applications [23, 24]. Naturally comparisons with these high-precision generators became parts of test-beds for PHOTOS package.

The principle of PHOTOS algorithm is to replace, on the basis of well defined rules, the decay vertex embedded in the event record such as HEPEVT [52] or HepMC [53] with the new one, where additional photons are added. Such solution, initially not aimed for high precision simulations, turned out to be very effective and precise as well. Phase space parameterization was carefully documented in [26]. Gradually for selected decays [2628, 54], also exact matrix elements were implemented and could be activated in place of universal kernels.Footnote 4 Originally [23], only single photon radiation was possible and approximations in the universal kernel were present even in the soft photon region. With time, multiphoton radiation was introduced [25] and then installation of exact first order matrix elements in W and Z decays became available with C++ implementation of PHOTOS [29]. The algorithm of PHOTOS is constructed in such a way, that the same function, but with different input kinematical variables, is used if the single photon emission or full multiphoton emission is requested. Such an arrangement enables tests in a rigorous first order emission environment. For multiphoton emission, the same kernel is used iteratively, thanks to the factorization properties. Technical checks are thus spared. Optimal solution for the iteration was chosen and verified with alternative calculations [55, 56] based on the second order matrix element. It was later extended to the multiphoton case for Z decays in Ref. [27]. Numerical tests of that paper, for distributions of generic kinematical observables pointed to the theoretical precision for the simulation of photon bremsstrahlung of the 0.1 % level.

When presenting numerical results from PHOTOS we always refer to its C++ version [29] with matrix elements for W and Z decays switched on. If the LO level is explicitly mentioned, matrix elements are replaced by universal kernels and algorithm as of FORTRAN version 2.14 [25], or higher is used. At present, version [29] represents the up-to-date version of PHOTOS; it is available, with few technical updates, from LCG library [57] as well.

3 Definitions and results of the first order calculations

As a first step of our tests we have compared numerical results obtained from PHOTOS and from SANC programs in case of the single photon emission. These tests cross-check conventions used and numerical stability of the two calculations. They also verify the proper choice of parameters in PYTHIA8 generator, which is used to produce electroweak Born level events, on which PHOTOS is activated. We have monitored distributions for the following observables: pseudorapidity η of , transverse mass M T of \(\ell^{-}\bar{\nu}_{\ell}\) pair, and transverse momentum p T of in the case of charged current; and pseudorapidity η of , invariant mass M of + pair, and transverse momentum p T of in the case of neutral current.

The following set of input parameters was used:

$$\begin{aligned} &{G_\mu = 1.6637 \times 10^{-5}~\mathrm{GeV}^{-2},} \\ &{M_W = 80.403~\mathrm{GeV}, \qquad \varGamma_W = 2.091~\mathrm{GeV},} \\ &{M_Z = 91.1876~\mathrm{GeV}, \qquad \varGamma_Z = 2.4952~\mathrm{GeV},} \\ &{V_{ud} = 0.9738, \qquad V_{us} = 0.2272,} \\ &{V_{cd} = 0.2271, \qquad V_{cs} = 0.9730,} \\ &{m_e = 0.511~\mathrm{MeV}, \qquad m_\mu = 0.10566~\mathrm{GeV},} \end{aligned}$$
(2)

and the following experiment motivated cuts were applied on momenta of the final-state leptons:

$$\begin{aligned} &{\mathrm{N}\mathrm{C}{:}\quad \big|\eta\bigl(\ell^{+}\bigr)\big| < 10, \qquad \big|\eta\bigl(\ell^{-}\bigr)\big| < 10,} \\ &{\phantom{\mathrm{N}\mathrm{C}{:}\quad }p_\perp\bigl( \ell^{+}\bigr) > 15~\mathrm{GeV},\qquad p_\perp\bigl(\ell^{-}\bigr) > 15~\mathrm{GeV},} \\ &{\phantom{\mathrm{N}\mathrm{C}{:}\quad }70 < M\bigl(\ell^{+}\ell^{-}\bigr) < 110~\mathrm{GeV};} \\ &{\mathrm{C}\mathrm{C}{:}\quad \big|\eta\bigl(\ell^{-}\bigr)\big| < 10, \qquad p_\perp\bigl(\ell^{-}\bigr) > 0.1~\mathrm{GeV},} \\ &{\phantom{\mathrm{C}\mathrm{C}{:}\quad }p_\perp(\bar{\nu}_\ell) > 0.1~\mathrm{GeV}.} \end{aligned}$$
(3)

We work in the running width scheme for W and Z boson propagators and fix value of the weak mixing angle: cosθ W =M W /M Z , sin2 θ W =1−cos2 θ W . The value of the electromagnetic coupling α is evaluated in the G μ -scheme using the Fermi constant G μ : the effective coupling is defined by \(\alpha_{G_{\mu}} = \sqrt{2}G_{\mu}M_{W}^{2}\sin^{2}{\theta_{W}}/\pi\).

To compute the hadronic cross section we have used CTEQ6L1 set of parton distribution functions with running factorization scale \(\mu_{r}^{2} = \hat{s}\), where \(\hat{s}\) is squared total energy of the colliding partons in their center-of-mass system.

Comparison is performed for muon and electron final states. For the electron final states and for each observable the bare and calo results are provided. In the calo case the four-momenta of the final electron and photon are combined into effective four-momentum of the electron when the separation

$$\Delta R = \sqrt{(\Delta\eta(e,\gamma))^2 + (\Delta\phi(e,\gamma))^2} < 0.1. $$

To check that the normalization of Born cross sections is properly adjusted between simulations using SANC and PHOTOS, we have completed tests at a sub-permille precision level. The corresponding results for electrons and ratios of differential cross sections are shown in Fig. 1 for CC and in Fig. 2 for NC. Errorbars in this plot and in all that follow represent the statistical fluctuations of the corresponding Monte Carlo integration.

Fig. 1
figure 1

Ratios for Born-level distributions in W decay

Fig. 2
figure 2

Ratios for Born-level distributions in Zee decay

For SANC and PYTHIA+PHOTOS cases, the results for \(\mathcal{O}(\alpha)\) corrections which are defined by \(\delta = (\sigma^{\mathcal{O}(\alpha)}-\sigma^{\mathrm{Born}})/\sigma^{\mathrm{Born}}\) are presented on Figs. 38. As we can see from these figures agreement at the one-loop level was found to be excellent, at the level of 0.01 % for both Z and W decays, once biases due to the technical parameter separating hard and soft photon emission were properly tuned between the two calculations.

Fig. 3
figure 3

\(\mathcal{O}(\alpha)\) corrections for basic kinematical distributions from PYTHIA+PHOTOS and SANC in W decay

Fig. 4
figure 4

\(\mathcal{O}(\alpha)\) corrections for basic kinematical distributions from PYTHIA+PHOTOS and SANC in Wμν decay

Fig. 5
figure 5

\(\mathcal{O}(\alpha)\) corrections for basic kinematical distributions from PYTHIA+PHOTOS and SANC in W decay (calo electrons)

3.1 Dependence on technical parameters

While performing tests we had to address well known technical problem of the so called “k 0 bias”. In case of fixed order correction implemented into Monte Carlo algorithm, a threshold on energy for emitted photon, typically in the rest-frame of the decaying particle, has to be introduced. It regularizes infrared singularity. Below this threshold photons are simply integrated out and resulting contribution is summed up with virtual corrections to cancel the infrared singularity. Unfortunately k 0→0 limit can not be reached, unless one accept working with negative event weights. The dependence on the technical parameter k 0 is however small for inclusive quantities such as the total cross section. The effect becomes enhanced on differential distributions near the resonance peaks. Such an effect can be observed in Fig. 9 for invariant mass m μμ in Zμμ(γ) decay. Generally this dependence is nowadays of no interest, since in practical application options of the programs with multiple photon emission should be used. This example however may be instructive for studies of ambiguities in implementation like in [58] where PHOTOS is used for soft photon emission only, while the hard photon phase space is populated with the help of genuine POWHEG simulation for the final and initial state emissions simultaneously.

4 Multiple photon emissions

Let us now turn to the same observables, but calculated with the multiple photon emission option of SANC and PHOTOS, suitable for the actual comparisons with the data. One can see from Figs. 10, 11, 12 and 13, that agreement is a bit worse than for the single photon case, but well within the expected theoretical precision. The relative contribution of higher order corrections is defined as \(\delta = (\sigma^{\mathcal{O}(\alpha^{2}+\mathrm{higher})}-\sigma^{\mathcal{O}(\alpha)})/\sigma^{\mathrm{Born}}\), so that it can be directly summed with the first order effect considered in the previous section. One should stress that these results represent quantitative comparisons of different approximations used in the SANC and PHOTOS as well. In fact, the approximations used, contrary to the single photon case, are not identical. Results of PHOTOS represent the NLO calculations with exponentiation and resummation of the collinear terms of the first order photon emission matrix element. Results of SANC use the collinear leading logarithmic approximation which introduces by construction numerical dependence on the corresponding QED factorization scale μ 2. Differences due to non-optimal choice of the scale μ 2 in SANC are below several permille for differential distributions and below 0.1 % otherwise, see Figs. 10, 11, 12 and 13. Additional effort and care in estimation of the size of seemingly minor effects may be required, if further improvements on theoretical precision, beyond 0.1 % are needed.

One should keep in mind that in precision measurements of LHC experiments unfolding procedure is applied (see e.g. Table 1 and discussion in Sect. 7 of Ref. [65]) to obtain the so called dressed leptons (charged leptons and accompanying them collinear photons are recombined into a single effective lepton). In this way the conditions of the Kinoshita–Lee–Nauenberg theorem [59, 60] are fulfilled and the large corrections proportional to logarithms of the lepton mass cancel out. Numerically this effect can be seen from a comparison of Figs. 6 and 8 for the first order RC. Such a reduction happens with the higher order photonic corrections as well.

Fig. 6
figure 6

\(\mathcal{O}(\alpha)\) corrections for basic kinematical distributions from PYTHIA+PHOTOS and SANC in Zee decay

Fig. 7
figure 7

\(\mathcal{O}(\alpha)\) corrections for basic kinematical distributions from PYTHIA+PHOTOS and SANC in Zμμ decay

Fig. 8
figure 8

\(\mathcal{O}(\alpha)\) corrections for basic kinematical distributions from PYTHIA+PHOTOS and SANC in Zee decay (calo electrons)

Fig. 9
figure 9

Ratio of invariant mass distributions from SANC and PHOTOS in Z decay as function of \(k_{0}=\epsilon=2E_{\gamma,min}/\sqrt{s}\)

Fig. 10
figure 10

Higher order corrections for basic kinematical distributions from PYTHIA+PHOTOS and SANC in W decay

Fig. 11
figure 11

Higher order corrections for basic kinematical distributions from PYTHIA+PHOTOS and SANC in Wμν decay

Fig. 12
figure 12

Higher order corrections for basic kinematical distributions from PYTHIA+PHOTOS and SANC in Zee decay

Fig. 13
figure 13

Higher order corrections for basic kinematical distributions from PYTHIA+PHOTOS and SANC in Zμμ decay

4.1 Comparisons of PHOTOS with KKMC

In Ref. [27] we have demonstrated physics reasons behind very good, 0.1 % level, agreement between PHOTOS and KKMC results for final-state photonic bremsstrahlung. Still until now, only in case of Z/γ intermediate state rigorous tests with second order QED matrix element Monte Carlo are available, and only in restricted condition of no hadronic activities in the initial state. For the reference simulations of \(q \bar{q} \to Z/\gamma^{*} \to l^{+} l^{-} (n \gamma)\) processes KKMC Monte Carlo [14] was used.Footnote 5 Also the Born level \(q \bar{q} \to Z/\gamma^{*} \to l^{+} l^{-} \) events from KKMC were used for simulations with PHOTOS. For both cases the monochromatic series of events with fixed virtuality of intermediate Z/γ state have been generated. This provides source of particularly valuable benchmarks as KKMC is the only program which features exclusive exponentiation combined with spin amplitudes for double photon emissions. As the numerical results of such quite extensive tests, more than 1000 figures are collected on our web page [62]: plots for Zμ + μ () and Ze + e () are presented there.

For all plots of collection [62], before selection cuts are applied, events are boosted to the laboratory frame assuming fixed 4-vector of Z defined by its \(p_{T}^{Z}\) and η Z; the two dimensional grid in \((p_{T}^{Z},\eta^{Z})\) is constructed, with 7 bins in \(p_{T}^{Z}\) spanning region 0–50 GeV and 4 bins in η Z spanning from 0–2, for each point in the grid 40M events are simulated. Kinematical selection is applied on lepton and photon 4-momenta and kinematical distributions are constructed from accepted events. It is required that each lepton has \(p_{T}^{l}>20\) GeV and |η ±|<2.47, with the gap 1.37<|η ±|<1.52 excluded. The gap in η ± corresponds to transition region between central and end-cap calorimeter in ATLAS detector, and is used here to somewhat arbitrary enhance possible effect of exclusive selection. Then, straightforward comparison between electron and muon cases is available; the only difference originating from leptons’ masses. For photon |η γ|<2.37 and region 1.37<|η γ|<1.52 is excluded again, \(p_{T}^{\gamma}>15\) GeV is required. In example shown in Fig. 14, angle between photon and a closer lepton is shown for \(p_{T}^{Z}=9\) GeV and η Z=2. Results from KKMC and PHOTOS NLO are compared. In both cases multiphoton emission is generated in Z/γ rest frame for \(u \bar{u} \to Z/\gamma^{*} \to \mu^{+}\mu^{-} (n\gamma)\) production process, with no initial state activity of any sort and virtuality of intermediate state M Z +6 GeV=97.187 GeV, where M Z is the Z-boson mass.

Fig. 14
figure 14

The distribution of the angle in the laboratory frame between photon and the closer muon, as generated from KKMC (thick line) and PHOTOS (thin line); selection cuts are applied, see the text. The intermediate state of virtuality M Z +6 GeV=97.187 GeV with the transverse momentum p T =9 GeV and pseudorapidity η=2 was created in the \(u\bar{u}\) annihilation. The samples of 40M events were simulated. Rates of events with photon passing selection cuts defined in the text, differ for the two programs by a factor 0.9991. Relatively large (0.1 %) difference to unity is typical for the larger values of pseudorapidity. For LO results (see web page [62] for extended results) the ratio for the surfaces, is of the same order, for electrons it is closer to 1. The distribution would feature plateau if the angle and cuts were calculated in the rest frame of Z/γ state

For each choice of \(p_{T}^{Z}\) and η Z used to define Z/γ state 4-momentum a set of three observables: angle between photon and closer lepton, directions of leptons and an overall acceptance rate is monitored in [62]. A general agreement of 0.1 % can be concluded for all distributions. The results for LO restricted PHOTOS are collected there as well.

For all cases one has to bear in mind that overall normalization correction factors for cross section, like \((1+ \frac{3}{4}\frac{\alpha}{\pi})\) in the case of Z decay, have to be included when using PHOTOS package.

4.2 The \(\phi^{*}_{\eta}\) observable

Motivated by [63, 64] and by recent ATLAS publication [65] on precise observable \(\phi^{*}_{\eta}\), representing important improvement for the measurement of Z boson transverse momentum (\(p^{Z}_{T}\)) at LHC, we have decided to devote section of this paper to discussion on the respective QED FSR corrections.

The measurement of the Z boson transverse momentum (\(p^{Z}_{T}\) or \(\phi^{*}_{\eta}\)) offers a very sensitive way for studying dynamical effects of the strong interaction, complementary to the measurements of the associated production of bosons with jets. The knowledge of the \(p^{Z}_{T}\) distribution is crucial also to improve the modeling of the W boson production needed for a precise measurement of the W mass [6], in particular in the low \(p^{Z}_{T}\) region which dominates the cross section. The study of the low \(p^{Z}_{T}\) spectrum (\(p^{Z}_{T}< M_{Z}\)), has also an important implication for the understanding of the Higgs signatures [1] as well as for the New Physics searches at the LHC [66].

The precision of the direct measurement of the spectrum at low \(p^{Z}_{T}\) at the LHC and at the Tevatron using the Z leptonic decay is limited by systematic uncertainties related to the knowledge and unfolding of the experiments resolution, in particular lepton energy scale [67, 68].

In recent years, additional observables with better experimental resolution and less sensitive to experimental systematic uncertainties have been investigated [6972]. The optimal experimental observable to probe the low \(p^{Z}_{T}\) domain of Z/γ production at hadron colliders was found to be \(\phi^{*}_{\eta}\) which is defined as

$$ \phi^*_\eta= \tan(\phi_{\mathrm{acop}}/2) \cdot \sin\bigl( \varTheta^*_\eta\bigr), $$
(4)

where ϕ acop is an azimuthal opening angle between the two leptons and the angle \(\varTheta^{*}_{\eta}\) is the scattering angle of the leptons with respect to the proton beam direction in the rest frame of the dilepton system. The \(\varTheta^{*}_{\eta}\) angle is defined as

$$ \cos\bigl(\varTheta^*_\eta\bigr) = \tanh \biggl(\frac{\eta^{-} - \eta^{+}}{2} \biggr), $$
(5)

where η and η + are the pseudorapiditities of the negatively and positively charged lepton, respectively. The variable \(\phi^{*}_{\eta}\) is correlated to the quantity \(p^{Z}_{T}/m_{ll}\), where m ll is the invariant mass of the lepton pair. It therefore probes the same physics as the transverse momentum \(p^{Z}_{T}\) and can be approximately related to it by \(p^{Z}_{T} \simeq M_{Z} \cdot \phi^{*}_{\eta}\). From the experimental point of view the variable \(\phi^{*}_{\eta}\) relies entirely on the angle reconstruction of the leptons in pair production, therefore on the tracking devices of high precision.

We have studied the theoretical error on \(\phi^{*}_{\eta}\) distributions due to photonic bremsstrahlung effects. As in Sect. 4.1 comparison of results from PHOTOS and KKMC generators was performed and stored on web page [62]. Let us list the appropriate cuts and give sample results.

We have requested that for both leptons \(p_{T}^{l}>\) 20 GeV and |η ±|<2.4. Distributions of \(\frac{dN(Z\to l^{+}l^{-})}{d\phi^{*}_{\eta}}\) from KKMC and PHOTOS were monitored. As before, the monochromatic samples of \(q \bar{q} \to Z/\gamma^{*} \to l^{+}l^{-} (n\gamma)\) were generated for two virtualities: M Z +6 GeV (97.187 GeV) and M Z −4 GeV (87.187 GeV), for incoming up and down quarks. Generated events were boosted to the laboratory frame assuming fixed \(p_{T}^{Z}\) and η Z of intermediate Z/γ state. Again the grid of 7 bins in \(p_{T}^{Z}\) spanning region 0–50 GeV and 4 bins in η Z spanning from 0–2 was used.

In Fig. 15 an example for one point of \((p_{T}^{Z},\eta^{Z})\) grid is chosen and the \(\phi^{*}_{\eta}\) distribution is shown. An agreement significantly better than 0.1 % between KKMC and PHOTOS is observed. For other choices of flavor of incoming quarks, Z/γ momentum and virtuality agreement is equally good [62].

Fig. 15
figure 15

The distribution of the \(\phi^{*}_{\eta}\) in Z/γ μ + μ () as generated from KKMC and PHOTOS; selection cuts are applied, see the text. The intermediate state of virtuality 97.187 GeV with the transverse momentum p T =9 GeV and pseudorapidity η=2 decaying into muon pair was created in the \(u\bar{u}\) annihilation. The samples of 40M events were used, ratio of the surfaces under distributions is 0.9991. Relatively large (0.1 %) difference to unity is typical for the larger value of pseudorapidity. For LO results (see web page [62] for extended results) the ratio for the surfaces, is of the same order. For electrons, both in LO and NLO cases, this ratio is closer to 1

4.3 Case of universal kernel

So far, in all numerical results PHOTOS with first order matrix elements as available in C++ version were used both in Z and W decays. If only the universal kernel was used, as available in public FORTRAN PHOTOS version 2.14 or higher, the loss of precision would be noticeable, but the uncertainty calculated with respect to the total rate would remain at 0.2 % level, for photonic bremsstrahlung corrections to the shapes of distributions [62]. As in the NLO case, overall normalization factor has to be corrected separately.

5 Non-photonic final-state bremsstrahlung

In this section we concentrate on those effects which go beyond multiple photon emissions. They can be divided into three groups. Emission of additional pairs, the effect which certainly belongs to final-state emissions, effect of interference of initial-final-state QED effects and finally all effects which are not directly related to final-state radiation, but nonetheless may affect their matrix element calculations.

5.1 Emission of pairs

Emission of light fermion pairs should be included starting from the second order of QED, i.e. from the \({\mathcal{O}}(\alpha^{2})\) corrections. There are two classes of diagrams which need to be taken into account. Emission of real pairs (Fig. 16) and the corresponding correction to the vertex (Fig. 17). These two effects cancel each other to a large degree due to the Kinoshita–Lee–Nauenberg theorem. The generic size of the effect can be expected to be of the order of higher order photonic bremsstrahlung corrections discussed so far. Moreover, it is well known from direct calculations in particular cases that the pair corrections are typically several times smaller than the photonic bremsstrahlung ones in the same order in α. Let us recall that careful studies of pair radiation effects were performed at LEP times [39, 73].

Fig. 16
figure 16

A typical example of real pair correction

Fig. 17
figure 17

A typical example of virtual pair correction

In the context of this paper we will estimate theoretical uncertainty due to neglecting of pair effects in the Drell–Yan observables.

The SANC integrators allow to perform a quick calculation of the pair corrections within the leading logarithmic approximation. We apply here the formalism of electron structure (fragmentation) functions [38, 44, 45] which describe radiation in the approximation of collinear kinematics.

The leading logarithmic approximation (LLA) was applied to take into account the corrections of the orders \({\mathcal{O}}(\alpha^{n}L^{n})\), n=2,3. Let us remind that in the first order \({\mathcal{O}}(\alpha)\) SANC has the complete calculation. The large logarithm \(L=\ln(\mu^{2}/m_{l}^{2})\) depends on the lepton mass m l and on the factorization scale μ. The latter is taken to be equal to the c.m.s. incoming parton energy (other choices are also possible).

The pure photonic contribution to the non-singlet electron fragmentation function in the collinear leading logarithmic approximation reads:

$$\begin{aligned} {\mathcal{D}}^{\gamma}_{ee}(y,L) =& \delta(1-y) + \frac{\alpha}{2\pi}(L-1)P^{(1)}(y) \\ &{}+ \frac{1}{2} \biggl(\frac{\alpha}{2\pi}(L-1) \biggr)^2P^{(2)}(y) \\ &{}+ \frac{1}{6} \biggl(\frac{\alpha}{2\pi}(L-1) \biggr)^3P^{(3)}(y) + {\mathcal{O}}\bigl(\alpha^4L^4\bigr). \end{aligned}$$
(6)

Analytic expressions for the relevant higher order splitting functions can be found in Refs. [38, 45].

For numerical evaluations of the splitting functions regularized by the plus prescription, we applied the phase space slicing as follows:

$$\begin{aligned} P^{(1)}(y) =& \biggl[\frac{1+y^2}{1-y} \biggr]_+ \\ =& \lim _{\Delta\to 0} \biggl\{ \delta(1-y) \biggl(2\ln\Delta+\frac{3}{2} \biggr) \\ &{} + \varTheta(1-y-\Delta) \frac{1+y^2}{1-y} \biggr\} , \end{aligned}$$
(7)

where an auxiliary small parameter Δ is introduced. In actual computations we used Δ=10−4 and verified the independence of the numerical results from the variations of this parameter, definition of P i(y), i>1 is given in Refs. [38, 45].

The effect due to emission of real and virtual electron-positron pairs can be estimated using the non-singlet and singlet pair contributions to the LLA electron fragmentation functions:

$$\begin{aligned} {\mathcal{D}}^{\mathrm{pair}}_{ee}(y,L) =& \biggl(\frac{\alpha}{2\pi}(L-1) \biggr)^2 \biggl[\frac{1}{3}P^{(1)}(y) + \frac{1}{2}R^s(y) \biggr] \\ &{}+ \biggl(\frac{\alpha}{2\pi}(L-1) \biggr)^3 \\ &{} \times \biggl[ \frac{1}{3}P^{(2)}(y) + \frac{4}{27}P^{(1)}(y) \\ &{}+ \frac{1}{3}R^{s}\otimes P^{(1)}(y) - \frac{1}{9}R^s(y) \biggr] + {\mathcal{O}}\bigl( \alpha^4L^4\bigr). \end{aligned}$$
(8)

Expressions for the relevant singlet splitting functions R s and R sP (1) can be found in Ref. [45].

The differential cross section of the neutral current Drell–Yan process with FSR leading logarithmic corrections takes the factorized form

$$\begin{aligned} &{d^5\sigma^{\mathrm{FSR}}_{\mathrm{LLA}}\bigl(p+p\to X+ l^+ \bigl(y_1p^+\bigr)+l^-\bigl(y_2p^-\bigr)\bigr)} \\ &{\quad=d^3\sigma^{\mathrm{Born}}\bigl(p+p\to X+ l^+\bigl(p^+\bigr)+l^- \bigl(p^-\bigr)\bigr)} \\ &{\qquad {}\times d y_1 \bigl({\mathcal{D}}^{\gamma}(y_1,L) +{\mathcal{D}}^{\mathrm{pair}}(y_1,L) \bigr)} \\ &{\qquad {}\times d y_2 \bigl({\mathcal{D}}^{\gamma}(y_2,L) +{ \mathcal{D}}^{\mathrm{pair}}(y_2,L) \bigr).} \end{aligned}$$
(9)

For consistency we expand the product of the fragmentation functions in α and take in only terms of the order \({\mathcal{O}}(\alpha^{2}L^{2})\) and \({\mathcal{O}}(\alpha^{3}L^{3})\). Terms coming from the product of the photonic and pair parts of the fragmentation functions are treated as a part of pair corrections.

Numerical results for the contribution of pair corrections are presented in Figs. 18, 19, for W decays.

Fig. 18
figure 18

Higher order photonic and pair corrections (δ in %) for basic distributions from PYTHIA+PHOTOS and SANC in \(W^{-}\to \mu^{-} \bar{\nu}\) decay

Fig. 19
figure 19

Higher order photonic and pair corrections (δ in %) for basic distributions from PYTHIA+PHOTOS and SANC in \(W^{-} \to e^{-} \bar{\nu}\) decay

If further improvement in precision would be required, pair emission can be implemented e.g. into C++ PHOTOS generator [29].

We assume that in the experimental analysis, there is no specific implicit cut rejecting \(Z \to l^{+}l^{-} f \bar{f} \) events. Our cut on \(p_{T}^{l}\) may reject some \(Z \to l^{+}l^{-} f \bar{f} \) events, since with the real \(f \bar{f} \) pair emission, leptons l will have a somewhat softer energy spectrum. But this require relatively high energy to be carried out by the \(f \bar{f} \) pair, the dominant triple logarithmic term will thus cancel between contributions from diagrams of Figs. 16 and 17. Let us stress that special care is needed in case of selection cuts sensitive to the soft, small virtuality and collinear to primary lepton pairs. Only this region of phase space may contribute significantly to pair corrections. Cuts affecting additional leptons of larger energies are thus of no concern.

The direction of the leptons originating from Z decays (and therefore our \(\phi^{*}_{\eta}\) observable) may be affected by partial reconstruction of two leptons (or two hadronized quarks) of the extra pair. Resulting phenomena may be important only if the pair \(f \bar{f} \) is not collinear to any of the primary leptons l. Again this represents a non dominant effect, thus substantially below required precision goal of a permille level.

Summarizing, details of pair corrections are still not important for precision tag at the level of 0.1–0.2 %.

5.2 Initial-final-state bremsstrahlung interference

The effect of QED-type interference between spin amplitudes for emission from the initial and final states represents important, even if numerically small, class of corrections. Even if separation of the initial and final-state bremsstrahlung at the spin amplitude level is clear, large interference effects may make this separation of limited practical convenience. However in certain approaches interferences can be combined with QED final-state effects in a convenient way.

Before we will present numerical results, let us recall first some details of the discussion from LEP times, see Ref. [73]. The discussion was devoted to energies higher than resonance peak, that is why interference cancellations as explained e.g. in [74, 75] did not apply. In the case of LHC observables discussed in this paper, oriented on production of W and Z resonances, such suppression is nonetheless expected. The reason is of a physical origin; time separation between boson production and decay. However, as a consequence of the uncertainty principle this suppression can be broken with strong event selection cuts if these cuts would constrain the final-state energies.

In the studies of Drell–Yan processes at LHC one can restrict discussion of the interference to the first order only. On the technical level control of the QED \({\mathcal{O}}(\alpha)\) interference contribution is realized in the SANC Monte Carlo integrator rather simply. The corresponding effect is computed by switching a respective flag in the code.

For the W decays similar arguments related to intermediate state life-time apply. In this case however, size of corrections is calculation scheme dependent, i.e. depends on the way how diagrams of photon emission off W line are treated. Studies with SANC demonstrated that the interference is below 0.1 % for LHC applications. They were completed not only for W but for Z as well, see e.g. [76].

Let us show in Fig. 20, as an example, corrections from interference to the \(\phi^{*}_{\eta}\) observable discussed previously. The effect is small, below 0.1 % for \(\phi^{*}_{\eta}<0.15\) and rises to 0.5 % for \(\phi^{*}_{\eta}\simeq 0.3\). After combining the interference and FSR corrections with QCD parton showers (for this purpose we used SANC Monte Carlo generator interfaced to Pythia8 program) the effect becomes significantly smaller (see Fig. 21). This is due to the fact that parton showers strongly affect the distribution of \(\phi^{*}_{\eta}\) variable, increasing the population of bins with \(\phi^{*}_{\eta}> 0.004 \) by a few orders of magnitude.

Fig. 20
figure 20

IFI/FSR ratio in Z decay for ϕ distribution. For \(\phi^{\ast}_{\eta}>0.2\) interference effects become sizable

Fig. 21
figure 21

IFI/FSR ratio in Z decay for ϕ distribution combined with parton showers

Recent measurement from ATLAS collaboration [65] of \(\phi^{*}_{\eta}\) observable extends to 1.3 but with large so far, statistical errors in range \(\phi^{*}_{\eta}>0.2 \). If precision requirements would become more demanding, the effect should be included together with the FSR corrections. At present, size of the interference effect can be used to estimate the size of the corresponding theoretical uncertainty due to its omission.

We can conclude that the initial-final-state interference does not represent a problem for separating final-state photonic bremsstrahlung from the remaining electroweak corrections in processes of W and Z production at LHC. This conclusion is justified for the precision of 0.1 %, but it will have to be studied in more detail for more exclusive configurations, like e.g. larger values of \(\phi^{\ast}_{\eta}\) distribution or for observables defined for off resonance peak regions of lepton pair invariant masses.

5.3 Relations with other electroweak and hadronic interactions

In calculation of final-state radiation matrix element, dependence on the direction of incoming quarks is present. However one can see from [27] that such effects due to e.g. initial state interactions, affect numerical results for final-state bremsstrahlung in a minor way, through the term which is in itself at the permille level, thus well below present precision requirement of 0.1 %.

For decays of narrow width states or when gauge symmetry can be used to separate phenomena from other parts of the interactions, there are no major difficulties to identify QED effects at the spin amplitude level. In the general case, QED FSR can be defined and its systematic error can be discussed as well. However, if e.g. contribution of diagrams featuring t-channel exchange of bosons complicate the separation, discussion of systematic errors of other parts of calculations may become scheme dependent. Our discussion on QED FSR corrections only will nonetheless be still useful for experimental applications.

6 Summary

We have addressed question of theoretical error for predictions of QED final-state bremsstrahlung in decays of W and Z bosons, used in precision measurements at LHC and Tevatron experiments.

Tests and comparisons of PHOTOS versus SANC programs for final-state photonic bremsstrahlung were performed in realistic conditions. Hard processes and initial state hadronic interactions were simulated with the help of PYTHIA8 [77] Monte Carlo program, for results with PHOTOS. For SANC its own set-up was used. Related differences in electroweak Born-level processes required careful tuning until agreement was established.

We have started our discussion with technical tests and results obtained at the first order. Separation into QED FSR and remaining electroweak corrections have been studied and verified at the spin amplitude level. We have checked that in PHOTOS and SANC numerically compatible, down to 0.01 % precision level, schemes of such separation were defined. This agreement confirms proper installation of matrix elements and numerical stability of SANC and PHOTOS as well. Then the comparison was repeated after allowing multi-photon emission in both programs. An agreement necessary to estimate systematic error in implementation of the QED final-state photonic bremsstrahlung for the precision level of 0.1 % was found. This conclusion holds for the decays of intermediate states, produced from annihilation of light quarks, predominantly close to the W and Z resonance peaks but with tails of the distributions taken into account. The conclusion holds if NLO kernel is active in PHOTOS and for SANC multiphoton option. For PHOTOS with the LO kernel theoretical precision is estimated to be 0.2 %. This conclusion is limited to effects resulting from shapes of distributions and for the selection cuts discussed in the paper. In principle, whenever new type of cuts is applied such comparison needs to be repeated. Effects on normalization have to be taken into account independently, either as part of genuine electroweak corrections (thanks to the proper choice of μ PW in W decays), or as an simple overall factor, like \((1+ \frac{3}{4}\frac{\alpha}{\pi})\) in case of the Z decay.

We have estimated the size of the higher orders QED photonic bremsstrahlung corrections using other programs. The KKMC [14, 15] Monte Carlo program of LEP era, featuring exclusive exponentiation and second order matrix element for final-state photonic bremsstrahlung was used for reference results. With the help of this program monochromatic intermediate Z/γ states of fixed virtuality were produced from annihilation of light quarks. This provided interesting test while grid of predefined values of (p T ,η) was populated, in particular for \(\phi^{*}_{\eta}\) observable. The differences versus NLO PHOTOS was found below 0.1 % (0.2 % for PHOTOS kernel restricted to LO only).

Contributions due to interference of the initial and final state QED radiation were found to be below 0.1 % for selected W and Z observables, as expected from the physics arguments. Separation of the final-state radiation from the remaining electroweak effects is of a practical importance as it facilitate phenomenological work. Our calculation schemes are convenient from that point of view. Interference effect was found to be below required precision level.

We estimate precision level of photonic final-state corrections at 0.1 %. With such precision tag separation of QED FSR from the rest of the process can be used for the sake of detector studies on final-state leptons. Such detector studies represent also a well defined segment in comparison of theoretical predictions with the measured data. One exception is \(\phi^{*}_{\eta}\) distribution in region of large \(\phi^{*}_{\eta}> 0.15\) if initial state parton shower is not taken into account. Already at \(\phi^{*}_{\eta}\simeq 0.3\) interference reach 0.5 %. In this region of phase space spin amplitudes for bremsstrahlung in the initial and final states become gradually of comparable size. Emission of additional pairs was discussed as well and a size of effect was estimated at 0.1 % level.

We estimate an overall systematic error for FSR implementation in PHOTOS and SANC at 0.2 % (0.3 % for PHOTOS with LO kernel). Further improvement of precision is possible, but requires more detailed discussion. Details of experimental acceptance have to be taken into account.

At a margin of the discussion we entered investigation of dependence on scheme specific parameters such as electromagnetic factorization scale μ 2 or photon energy threshold k 0=ϵ used in fixed order simulations. This may be of some interest for further studies of uncertainties resulting from some choices of matching of FSR with hard process and/or initial state interactions and/or hard emission matrix elements.