Separating electroweak and strong interactions in Drell–Yan processes at LHC: leptons angular distributions and reference frames

Among the physics goals of LHC experiments, precision tests of the Standard Model in the Strong and Electroweak sectors play an important role. Because of nature of the proton-proton processes, observables based on the measurement of the direction and energy of leptons provide the most precise signatures. In the present paper, we concentrate on the angular distribution of Drell–Yan process leptons, in the lepton-pair rest-frame. The vector nature of the intermediate state imposes that distributions are to a good precision described by spherical polynomials of at most second order. We show that with the proper choice of the coordinate frames, only one coefficient in this polynomial decomposition remains sizable, even in the presence of one or two high pT\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p_T$$\end{document} jets. The necessary stochastic choice of the frames relies on probabilities independent from any coupling constants. This remains true when one or two partons accompany the lepton pairs. In this way electroweak effects can be better separated from strong interaction ones for the benefit of the interpretation of the measurements. Our study exploits properties of single gluon emission matrix elements which are clearly visible if a conveniently chosen form of their representation is used. We rely also on distributions obtained from matrix element based Monte Carlo generated samples of events with two leptons and up to two additional partons in test samples. Incoming colliding protons’ partons are distributed accordingly to PDFs and are strictly collinear to the corresponding beams.


Introduction
The main purpose of the LHC experiments [1,2] is to search for the effects of New Physics. This programme continues after the breakthrough discovery of the Higgs boson [3,4] a e-mail: Z.Was@cern.ch and measurement of its main properties [5]. In parallel to searches of New Physics, see e.g. [6][7][8], a program of precision measurements in the domain of Electroweak (EW) and Strong (QCD) interactions is on-going. This is the keystone of the programme of establishing the Standard Model as a fundamental theory. It is focused around two main directions: searches (setting upper limits) for anomalous couplings and precision measurements of the Standard Model parameters. Precision measurements of the production and decay of intermediate Z and W bosons represent the primary group of measurements of the second domain, see e.g. [9][10][11][12].
In the following, let us concentrate on the phenomenology of distributions for leptons originating from Drell-Yan processes, and in particular on possible separation of electroweak and strong phenomena. There is a multitude of publications devoted to the theoretical preparation for such measurements, e.g. [13][14][15]. They are mostly discussing evaluation of higher order effects of strong interactions. The commonly used scheme for controlling the angular distribution of leptons is a Collins-Soper frame [16] motivated by the factorisation theorem [17].
However, it is long known that the distribution of leptons produced from a vector intermediate state such as e + e − → μ + μ − can be described as a sum of contributions consisting of Born level distributions with factorized terms independent of the relative direction of the leptons. The analytical form for these distributions is a result of an explicit first order calculation in O(α Q E D ), given for example in [18], formula (3.4). This expression is directly applicable to initial state emission of gluons in case of qq → + − hard process in proton-proton collisions as well. It is valid all over the phase space. Let us not go into details of the mass effects, that is the terms of the order of m 2 s . This complicates details, but does not change the general principles. In [19] it is argued that such result may hold in more general case as well. In particular, that they break only at the level of complete calculations for two hard emitted partons. Such considerations were essential in establishing the PHOTOS Monte Carlo [20] algorithm for bremsstrahlung in decays of resonances and particles from its very beginning. It is of interest if considerations of that sort can be useful for the phenomenology of initial state parton emissions as well.
In [21] it was shown, that the distributions of leptons in the rest frame of the lepton pairs, for the processes of pp collisions are described by spherical harmonics of at most second order. In the formulation of the decomposition into a sum of trigonometrical polynomials, up to 8 non-zero coefficients are expected, showing dependence on the lepton pair invariant mass, transverse momenta and rapidity. In contrast, at the Born level, only one coefficient with functional dependence on electroweak parameters suffices.
The question is, whether the wealth of 8 coefficients brings some more information on the hard electroweak process and on QCD dynamics or is a consequence of geometrical alteration due to presence of accompanying jets only; thus bringing no additional physical information. On the other hand, they are based on precisely measured kinematics of leptons, usually more precise than that of jets.
By inspection of formula (3.4) from [18], we find that in case of gluon (photon) emission the distribution of the lepton angles in the lepton pair system consists of two Bornlike contributions weighted respectively by squares of the energies of the two incoming partons as seen in rest-frame of the lepton pair system. Such weight depends only on the 2 → 3 process kinematics (no dependence on couplings) and can be used over the entire range of energy and momentum of the emitted parton. That offers easy application. The scattering angle is then the angle (also in the same frame) between the direction of outgoing leptons and respectively first/second incoming parton. The azimuthal angle can be chosen as the same one for the two Born's as of formula (3.4) from [18], as the angle between two oriented half-planes built from: (i) the direction of the beams and emitted parton and (ii) the direction of the outgoing leptons and emitted parton. It does not matter if such an angle is defined in the lepton pair frame or reaction frame effects of secondary decays, out of scope of the present publication, there is a freedom in definition of reference direction for azimuthal angle. Matrix element does not depend on it.
A similar solution was used for the PHOTOS Monte Carlo where results from [22,23] were instrumental.
Let us stress that the above strategy was developed in the first place for initial state QED bremsstrahlung. It is thus of interest to which degree it can be extended to the case of pp collisions where instead of bremsstrahlung photons, jets are emitted: not always gluonic ones, but sometimes of quark partons as well. For that purpose we study samples generated with the help of Monte Carlo programs featuring explicit matrix elements.
Our paper is organized as follows. Section 2 is devoted to the presentation of strong and electroweak interaction features which are important for the explanation of our results. In Sect. 3 we list properties and definitions important for the analysis, using an expansion of the angular distribution of decay products of Drell-Yan process in the rest-frame of the lepton pair. We elaborate on possible choices of the coordinate frame orientation. The choice of orientation can be optimized with the help of matrix element results or leading logarithm calculations; the variants are presented and discussed.
In Sect. 4 we collect numerical results for the distributions in the case of pp → l + l − j and pp → l + l − j j, for invariant mass of the lepton pair around the Z -boson peak. The aim is to demonstrate that with the help of probabilistic splitting of each event into configurations defined with respect to two distinct frames, the results of formula (3.4) from [18] are reproduced and indeed only one non-zero coefficient in the decomposition of the angular distribution is needed. We use series of pp → l + l − j events generated with MadGraph [24]. We first concentrate on a sample where the outgoing parton is a gluon, which nonetheless carries substantial transverse momentum, p T . In the consecutive results and subsections we release this restriction of outgoing gluon only. Later, we turn to pp → l + l − j j events. Again, a sample of events generated by MadGraph is studied. Finally, we provide also results for Powheg+MiNLO [25,26] generated sample of pp → Z j events at QCD NLO.
In the following Sect. 5 we discuss results obtained for different sub-groups of events, where different production mechanisms dominate, and how the consequences of choices for reference frames point to the general pattern. Finally, in Sect. 6 we conclude the paper.

QCD and EW
The predictive power of QCD is based on the factorisation theorem [17]. It provides a framework for separating out long-distance effects in hadronic collisions. In consequence it allows for systematic prescriptions and tools to calculate the short-distance dynamics perturbatively, the same time allowing for identifying the leading nonperturbative longdistance effects which can be extracted from experimental measurements or by numerical calculations in Lattice QCD.
Let us list some references to illustrate the developments in this area without ambitions of their completeness. Theoretical investigations of Drell-Yan pair production [27] have a long history. It is one of the few processes in hadron-hadron collisions where the collinear QCD factorization has been rigorously proven [28][29][30][31]. Within this framework, the NLO pQCD calculations of inclusive cross sections have been performed [32,33] and later it was done on up to NNLO accuracy [34,35]. In particular, fully exclusive NNLO pQCD calculations became available, including the leptonic decay of the intermediate Z boson [36][37][38].
The question of the input from Electroweak sector of the Standard Model is important, especially for distributions of leptons originating from intermediate Z /γ * state. We have addressed numerical consequences of this point recently in [39] in the context of τ lepton polarization in Drell-Yan processes at the LHC. A wealth of publications was devoted during last years to this point, see e.g. [15,40,41]. Let us point however to constraints from limitations to separating interactions into Electroweak and QCD ones, which are well known, since more than 15 years, see e.g. [42].
In the following, we will concentrate on numerical analysis of processes at tree level for parton-parton collisions into lepton pair and accompanying jets. Even though approach is limited to leading, tree-level predictions, it provides input for general discussions. Such configurations constitute parts of the higher order corrections, or can be seen as the lowest order terms but for observables of tagged high p T jets.

Angular decay distribution
The measurement of the angular distribution of leptons from the decay of a gauge boson V → where V = W, Z or γ * , produced in hadronic collisions via a Drell-Yan-type process h 1 + h 2 → V + X provides a detailed test of the production mechanism. If the gauge boson is produced at high transverse momentum, one can define an event plane spanned by the beam and the gauge boson directions. This reference plane can be used for study of lepton-hadron current correlation effects.
In fact, these correlations are described by the set of nine hadronic structure functions which can be calculated within the context of the parton model using perturbative QCD. Such measurements are already completed for Tevatron [43] and LHC [44,45].
Following the conventions and notations of [21,46], let us briefly recall that the lepton-hadron correlations are described by the contraction of the lepton tensor L μν with the hadron tensor at the parton level H μν , where L μν acts as an analyser of the structure of H μν which carries the effective information on the polarisation of the gauge boson produced in the interaction. The angular dependence can be extracted introducing helicity cross-sections corresponding to the nonzero combinations of the polarisation density matrix elements where m, m = +1, 0, −1 and are the polarisation vectors for the gauge boson defined with respect to the chosen lepton-pair (gauge boson) rest-frame. The angular dependence of the differential cross-section can be written as dσ dp 2 T dY d * = 9 where the g α (θ, φ) represent harmonic polynomials of the second order, multiplied by normalisation constants and dσ α denote helicity cross-sections, corresponding to nine helicity matrix elements. The angle θ and φ in d * = d cos θ dφ are the polar and azimuthal decay angles of the leptons in the lepton-pair rest-frame. We can conveniently rewrite Eq. (3), explicitly defining polynomials and corresponding coefficients dσ dp 2 T dY d cos θ dφ where dσ U +L denotes the unpolarised differential crosssection (a convention used in several papers of the 80's). The coefficients A i ( p T , Y ) are related to ratios of corresponding cross-sections of intermediate state helicity configurations [21].
The dynamics of the production process is hidden in the angular coefficients A i ( p T , Y ). This allows us to treat the problem in a model independent manner. In particular, as we will see, all the hadronic physics is described implicitly by the angular coefficients and it decouples (the A 4 is an exception) from the well understood leptonic and intermediate boson physics. As we have already pointed out, the spherical angles θ and φ define the orientation of the lepton direction in the lepton-pair rest-frame. Let us stress that the actual choice of the orientation of coordinate frames one uses represents an important topic; we will return to it later.
The lepton angular distribution in the lepton-pair restframe is determined by the gauge boson polarisation and by the range of the invariant mass of the lepton pairs, which determines indirectly the admixture of the Z /γ * interference to the cross-section. The general structure of the angular distribution is given by nine helicity cross-sections corresponding to the nine spin density matrix elements for the gauge boson [21]. In the LO parton subprocesses (qq → V → ) only the transverse polarisation gauge boson contributes to the helicity density matrix. At higher order, in general, all three polarisations of the gauge boson contribute. At O(α S ) in perturbative QCD the angular distribution is described by six helicity cross-sections, which, for fixed range of invariant mass of lepton pairs, are functions of the transverse momentum and rapidity of the gauge boson. At the O(α 2 S ), contributions from three additional helicity cross-sections are non-zero. It is important that higher order corrections do not break the above picture in numerically important terms.
As we will see, for the choices of the frames discussed here, in the limit of zero transverse momenta all coefficients, except A 4 which is defined by electroweak couplings, vanish. The ( p T , Y ) dependence of the A i coefficients vary strongly with the choice of the reference frame. Let us present now the variants of the reference frame definition we are going to use.

Collins-Soper frame
The well known and broadly used [16] Collins-Soper reference frame is defined as a rest-frame of the lepton-pair (Z -boson), with the polar and azimuthal angles constructed using proton directions in that frame. Since the Z -boson has a transverse momentum, the directions of initial protons are not collinear in the lepton-pair (Z -boson) rest frame. The polar axis (z-axis) is defined in the lepton-pair rest-frame such that it is bisecting the angle between the momentum of one of the proton and inverse of the momentum of the second one. The sign of the z-axis is defined by the sign of the lepton-pair momentum with respect to z-axis in the laboratory frame. To complete the coordinate system the y-axis is defined as the normal vector to the plane spanned by the two incoming proton momenta and the x-axis is chosen to set a right-handed Cartesian coordinate system with the other two axes. Polar and azimuthal angles are calculated with respect to the outgoing lepton and are labeled θ and φ respectively. In the case of zero transverse momentum of the lepton-pair, the direction of the y-axis is arbitrary. Note, that there is an ambiguity in the definition of the φ angle in the Collins-Soper frame. The orientation of the x axis here follows convention of e.g. [21,47,48].
The formula for cos θ can be expressed using directly momenta of the outgoing leptons in the laboratory frame [49]. with where E i and p z,i are respectively the energy and longitudinal momentum of the lepton (i = 1) and anti-lepton (i = 2) and p z ( + − ) denotes the longitudinal momentum of the lepton system, m( + − ) its invariant mass. The φ angle is calculated as an angle of the lepton in the plane of the x and y axes in the Collins-Soper frame.
Only the four-momenta of outgoing leptons and incoming proton directions are used. That is why the frame is very convenient for experimental purposes.

Mustraal frame
The Mustraal reference frame is also defined as a rest frame of the lepton pair. It have been proposed and used for the first time in Mustraal Monte Carlo program [18] for the parametrization of the phase space. The Mustraal Monte Carlo program was constructed for muon pair production at LEP. Since then, it has been used successfully, over LEP time, in other Monte Carlo programs as well. Resulting optimal frame was used to minimise higher order corrections from initial state radiation to the e + e − → Z /γ * → ff for the algorithm implementing genuine weak effects in the LEP era Monte Carlo program KORALZ [50]. A slightly different variant was successfully used in the Photos Monte Carlo program [51] for simulating QED radiation in decays of particles and resonances. The parametrization was useful not only for compact representation of single photon emissions but for multi-emission configurations as well.
In this paper, implementation of the Mustraal frame has been extended to the case of pp collision and one or two partons in the final state accompanying Drell-Yan production of the lepton pairs. Let us describe in more detail the implementation of this phase space parametrization. As we will use it for previously generated events, we will calculate angles and energies from the outgoing particles four-momenta, rather than the opposite way around.
The parametrization is constructed in consecutive steps, each following one, using information obtained from the previous.
• We start from the following information, which turns out to be sufficient: (i) The 4-momenta and charges of outgoing leptons τ 1 , τ 2 . (ii) The sum of 4-momenta of all outgoing partons. • The orientation of incoming beams b 1 , b 2 is fixed as follows: b 1 is chosen to be always along positive z-axis of the laboratory frame and b 2 is anti-parallel to z axis. The information on incoming partons of p 1 , p 2 is not taken from the event record. It is recalculated from kinematics of outgoing particles and knowledge of the center of mass energy of colliding protons. In this convention the energy fractions x 1 and x 2 of p 1 , p 2 carried by colliding partons, define also the 3-momenta which are along b 1 , b 2 respectively.
• The flavour of incoming partons (quark or antiquark) is attributed as follows: incoming parton of larger x 1 (x 2 ) is assumed to be the quark. This is equivalent to choice that the quark follow direction of the outgoing system, similarly as it is defined for the Collins-Soper frame. 1 This choice is necessary to fix sign of cos θ 1.2 defined later.
• The 4-vectors of incoming partons and outgoing leptons are boosted into lepton-pair rest frame. • To fix orientation of the event we use versorx lab of the laboratory reference frame. It is boosted into leptonpair rest frame as well. It will be used in definition of azimuthal angle φ, which has to extend over the range (0, 2π). • We first calculate cos θ 1 (and cos θ 2 ) of the angle between the outgoing lepton and incoming quark (outgoing antilepton and incoming anti-quark) directions.
• The azimuthal angles φ 1 and φ 2 corresponding to θ 1 and θ 2 are defined as follows. We first define e y 1,2 versors and with their help later φ 1,2 as: and similarly for φ 2 : • Each event contributes with two Born-like kinematics configurations θ 1 φ 1 , (θ 2 φ 2 ), respectively with wt 1 (and wt 2 ) weights; wt 1 + wt 2 = 1 where . (9) In the calculation of the weight, incoming partons energies E p1 , E p2 in the rest frame of lepton pair are used, but not their couplings or flavours. That is also why, instead of σ B (s, cos θ) the simplification (1 + cos 2 θ) is used in Eq. (9). Dependence on the sign of cos θ drops out. 2 • The above definitions of φ 1,2 are not used directly. The φ 1 and φ 2 are retrieved as acos(cos φ 1 ) and acos(cos φ 2 ) as explained above. However, to extend to the full range (0, 2π), the directions of the incoming partons p 1 , p 2 and outgoing parton (sum of outgoing partons) p 3 of the laboratory frame and τ 1 , τ 2 of the τ -pair rest frame are used to calculate only sign of the following For sign < 0, azimuthal angle φ = 2π − acos(cos φ) is taken. • Let us point out once again, that our weights wt 1,2 are independent from coupling constants of matrix element. The four-momenta of incoming partons are not needed; necessary information is obtained from final states. That is also why the solution is stable in cases where further sources of p T beyond one or two reconstructed jets would be present.
The algorithmic logic of the above construction, the extended Mustraal frame, is quite similar to the one developed for the PHOTOS Monte Carlo [20] or for the inclusion of genuine weak corrections (and final state bremsstrahlung) in KORALZ [50]. Experience from unfinished projects [54,55] was of importance as well.
Unlike the case of the Collins-Soper frame, the Mustraal frame requires not only information on 4momenta of outgoing leptons but also on outgoing jets (partons). The information on jets (partons), is used to approximate the directions and energies of incoming partons for 2 Stronger simplification is explored in [52]. As a consequence they predict e.g. different dominant dependence on electroweak couplings of A 7 than in [21]. In our work we evaluate an effect of approximation with the help of numerical studies of results from Monte Carlo programs based on explicit matrix element calculations. Our choice of frame does not depend on couplings present in the matrix elements. However it depends on some general geometrical properties of the emission amplitudes. That is why our observations are not in conflict with statements of reference [53]: we take into account dominant parts of intrinsic properties of the production mechanism because they are of a geometrical nature. calculation of wt 1,2 and θ , φ angles. This does not have to be very precise but can introduce additional experimental systematics, and thus requires attention.

Numerical results
Let us now present numerical results. We will use samples of events generated with the MadGraph matrix element Monte Carlo for Drell-Yan production of τ -lepton pairs, with m τ τ = 80 − 100 GeV and 13 TeV pp collisions. Lowest order spin amplitudes are used in this program for the parton level process. For the EW scheme we have used default initialisation of the MadGraph with on-shell definition of sin 2 θ W = 1 − m 2 W /m 2 Z = 0.2222, which determines value of the axial coupling for leptons and quarks to the Z-boson. The incoming partons are distributed accordingly to PDFs (using CTEQ6L1 PDFs [56] linked through LHAPDF v6 interface) and remain precisely collinear to the beams. Two samples of 2 and 20 M events were prepared, 3 the first one features τ τ j final states and the other one τ τ j j. At this level, jet (j) denotes outgoing parton of unspecified flavour. In the generated sample information on incoming and outgoing partons flavours are stored and we will use this information to define subsamples. Figures 1 and 2 show distributions of cos θ and φ calculated respectively in the Collins-Soper and Mustraal frames for pp → τ τ j process in case of qq annihilation or qG scattering. Clearly the asymmetric shape of cos θ is preserved despite an increasing threshold on p τ τ T transverse momenta and the φ modulation is almost completely removed with the choice of Mustraal frame. One can conclude from those two figures that in the Mustraal frame angular correlations remain robust against QCD initial state effects, preserving information on the V-A structure of electroweak couplings in the Z-boson vertex, not only in case of outgoing gluonic jet but also of outgoing quark. No deformation of the shape in cos θ distribution, which remains as at the Born level, is observed. In the case of the Collins-Soper frame, deformations are present. It was expected that choice of Mustraal frame preserves the Born shape in case of gluonic emission: in the definition of this frame properties of qq → τ τ j matrix elements were explored. The case of an outgoing quark, where a similar pattern is observed, is beyond what we would expect from the matrix element representation which was used for the Mustraal Monte Carlo construction. We will return to this point in Sect. 4.2. 3 The analysed sample of events show a rather large instability of the numerical predictions for the A i coefficients which by far exceeds statistical fluctuations. This may be the consequence of limitations due to pretabulations. That is why we have not pursued generation of higher statistic samples. Nonetheless, the above comparison is inspiring for more checks, namely if such properties will hold in more general cases, like when more partons are present in the final state, or when higher order corrections (QCD NLO) are taken into account. Let us continue with preparation of the phenomenological picture for such studies before continuing with further presentation of numerical results.

Angular coefficients in the Monte Carlo samples
The angular coefficients A i ( p T , Y ) are not explicitly input to the theoretical calculations nor to the MC event generators. They can however be extracted from the shapes of the angular distributions with the method proposed in [46], owing to the orthogonality of the spherical polynomials. The weighted average of the angular distributions with respect to any specific polynomial isolates an average reference value or moment of its corresponding coefficient. Reference [46] argues that only the coefficients of spherical harmonics of the second order will contribute, unless complete effects of α 2 s are taken into account. The moment of a polynomial P(cos θ, φ) over a specific range of p T , Y is defined as follows: As a consequence of Eq. (4) we obtain: Thanks to the discussion in [46], we expect that terms beyond Eq. (4), i.e. of higher order polynomials, should be of order 0.01 or smaller. Let us point out that two jets of large p T constitute contributions to corrections to the Drell-Yan process of such an order, so the effect at some corners of the phase space may be noticeable.
From now on, we will monitor the A i coefficients rather than the θ , φ distributions themselves. We will not discuss higher than second order coefficients. We have checked, that they are indeed small and remain within expected range. The size of those terms can be deciphered from the presented plots as well, that is why we do not devote to these point much attention. Our cross checks of those aspects are left unpublished, as at present they are of limited interest and are also easy to obtain.

Numerical results for configurations beyond single gluon emission
In the three Figs. 3, 4 and 5 we collect results for angular coefficients A i in the processes of one jet in the final state: first when it is from single gluon emission, later when it can be initiated by quark or anti-quark, finally when both subprocesses are combined. In the above and in the following plots, we show sets of five angular coefficients A 0 − A 4 only; the remaining ones A 5 − A 7 are close to zero over the full p τ τ T range for both definition of frames, Collins-Soper and Mustraal. We present also A 02 = A 0 − A 2 , it is known to be sensitive to the QCD sector and should be precisely zero at the LO QCD and the Collins-Soper frame, because of the Lam-Tung relation [57].
As the next step, we turn our attention to processes with two jets in final state, first of unspecified type, Fig. 6, and later for processes of quark anti-quark annihilation only, Fig. 7. Then the final state jets originate dominantly from a pair of gluons or from an intermediate gluon decaying to a quarkantiquark pair. These cases are similar from the point of view of spin structure to the single gluon emission. Finally, the last two figures are for the processes where one, Fig. 8, or two gluons, Fig. 9, are present in the initial state.
The separation between processes is based on the flavour of incoming partons (quark, anti-quark or gluon), as stored for each event in the generated sample. This allows splitting the sample of two jets events into groups of different spin nature of the intermediate states. Such separation cannot be universal, or applied to data of course, where we work with multiple production mechanisms. It is convenient for illustration purposes; it could be explored with more formal explanations following reasoning as in Refs. [59,60], see Sect. 5.2.

Results with NLO simulation
So far, we have discussed results for samples of fixed order tree level matrix elements of single or double parton (jet) emission. These results represent a particular sector of Drell-Yan processes. In general, configurations with a variable number of jets and effects of loop corrections and parton shower of initial state should be used to complete our studies. We have performed this task partially only, with the help of 40M weighted Z + j events generated with Powheg+MiNLO Monte Carlo, again for pp collisions at 13 TeV, invariant mass of lepton pair restricted to range 80-100 GeV, and the effective EW scheme 4 with sin 2 θ W = 0.23147 but otherwise of default initialization. The PowhegBox v2 generator [25,26], augmented with MiNLO method for choices of 4 Note that the axial coupling and as a consequence A 4 coefficient will be different than in case of MadGraph initialisation even in the zero transverse momenta limit. scales [61] and inclusion of Sudakov form factors [62], by construction achieves NLO accuracy for distributions involving finite non-zero transverse momenta of the lepton system.
In Fig. 10 results are completed for such a sample, following the same approach as in our previous plots. Comparisons of results using Mustraal and Collins-Soper frames feature again the usual pattern. This confirms the robustness of our conclusions. Further studies should include not only parton shower effects in the initial and final state but detector reconstruction as well. This is however beyond the scope of the present paper.

Main features of results
In all collected plots, we have compared results for the cases when Collins-Soper and Mustraal frames were used for the cross-section decomposition into angular coefficients and spherical polynomials. We can conclude that the first choice was better suitable to measure the effects of QCD dynamics, whereas in the second choice electroweak effects were far more dominating the picture, even in case of high p T jets. This may indicate which choice is better for the particular use. It is important that in both cases we rely on precisely measured leptons.
To confirm the origin of differences between Mustraal frame and Collins-Soper frame observed for A 0 , A 1 , A 2 and A 3 coefficients, we have analysed sample of pp → τ τ j events generated with MadGraph configured in the EW sector with sin 2 θ W = 0.25. In this way the A 3 and A 4 are set close zero over full p τ τ T range because the dominant contributions from the Z -boson peak was minimised, as they are proportional to a v ∼ (1/4−sin 2 θ W ) = 0 (in case of Mustraal A 3 was close to zero even for sin 2 θ W = 0.22222). We have observed that the differences in A 0 , A 1 and A 2 coefficients between Collins-Soper and Mustraal frames were quite similar as in the case of sin 2 θ W = 1 − m 2 W /m 2 Z = 0.22222, which has been used for numerical results presented for MadGraph. This observation suggests that differences in A 0 , A 1 and A 2 of two frames are due to QCD effects, and do not depend substantially on the EW sector. Therefore differences in A i can be used to provide a precision reference point for discussing the dynamics of QCD in the Drell-Yan processes, similarly as the φ * observable [9,63].
As stated already in [21], the theoretical uncertainties due to the choice of the factorisation and renormalization scales are very small for the cross section ratios A i . Most of the uncertainties of the structure functions and of the choice of the factorisation scheme cancels in the ratios too.

Angular coefficients in Collins-Soper frame
The coefficients A 0 and A 2 are increasing functions of p V T (where V denotes W or Z boson) and the deviations from lowest order expectation (A 0 = A 2 = 0) are quite large even at modest value of p V T , i.e. for p V T = 20 − 50 GeV. These coefficients are exactly equal at LO, thanks to Lam-Tung relation [57]. This is no longer true at NLO but predicted corrections are of less than 10 % [46]. They were recently measured [45], and it seems that in fact larger corrections are needed. This topic was recently covered in [58] at second order in perturbation theory. The deviation of A 1 and A 3 from zero is much smaller even at large p V T . At the Z -boson peak, the A 4 coefficient is proportional to the product of vector and axial couplings of the intermediate vector boson to leptons and quarks, (v /a v q /a q ), and is therefore sensitive to the Weinberg angle θ W .

Variations of angular coefficients if Mustraal frame is used
In general, for this choice of rest-frame, only A 4 is substantially different from zero and is fairly consistent with the same coefficient of Collins-Soper frame. The electroweak sector of the interactions dominates in the decomposition to spherical polynomials. Sensitivity to sin 2 θ W remains, but in contrast sensitivity to QCD effects is minimized. One may wonder, why the property of the matrix element of single gluon emission extends to other cases, when final states feature quark or there are two outgoing partons, gluonic or quarks. In principle, to understand this, one would need to study the properties of all spin amplitudes, e.g. following approach of [59,60]. For the purpose of the present paper, such an ambitious task is not necessary. We restrict ourselves to inspection of numerical results only. We are encouraged that more formal arguments may be provided if needed, up to the precision limits as explained in [19].
The QCD impact on angular coefficients in the Collins-Soper choice of frames is dominated by the jet(s) geometry. The observation can be made that the effects is due to the choice of frames orientation for which angles θ and φ are defined.
On the first look, one could observe that at LO QCD the Lam-Tung relation breaks in Mustraal case for qG scattering subprocess and holds in the case of Collins-Soper frame. For now we don't have explanation of this effect. However, at NLO QCD as we can observe from Fig. 10, the breaking of the rule itself is comparable for the two choices of the frames, especially at larger values of p τ τ T . One should not forget that the Mustraal choice is not, and should not be expected to set A 0 and A 2 to zero, for the parton level processes of incoming gluons, the two coefficients are nonetheless sizably smaller in this case.

Summary
The interest in decomposition of results for measurement of final states in Drell-Yan processes at the LHC into coef-ficients of second order spherical polynomials for angular distributions of leptons in the lepton-pair rest-frame, was recently confirmed by experimental publications [43][44][45]. Inspired by those measurements, we have investigated possible variants of choice for quantization frames orientation. We have exploited geometrical properties of matrix elements for vector quanta emissions known since early 80's.
We have found that the choice based on those properties helps to separate effects due to the Electroweak couplings and due to strong interactions: emissions of jets. Most of the coefficients can be eliminated or their size can be substantially reduced. The approach of stochastic attribution to each event orientation of reference frame, based on the properties of matrix element for single gluon emission was further simplified, to eliminate dependence on coupling constants and to become of a purely geometrical nature.
We have checked that such an approach retains sufficiently its nature and that it can be generalized to the case of events featuring more than one jet and/or of quark jets too. For that purpose, we have confronted numerical results of angular coefficients for all parton level processes with one or two outgoing jets (not necessarily gluonic one), for our Mustraal orientation of frame, with the one of Collins-Soper. We have found that with our choice, most of the coefficients in the expansion vanish or are reduced substantially, except for the one existing already at partonic Born level.
This may help to separate the interpretation of experimental results into quantities sensitive to electroweak and strong interaction effects. Because directions of leptons are measured with very high precision, measurement of angular coefficients can be more effective than direct measurement of hadronic jets, even if in principle, no more physics input is obtained in such a way. The gain of precision may be substantial, as in case of φ * η observable [9,63], but it requires more experimentally oriented studies. The approach of stochastic attribution to each event orientation of reference frame may be useful for the optimalization of TauSpinner methods [39,64], when (2 → 2) Born level amplitudes are used.
In this paper, we do not discuss the kinematic effects of real photon emissions and corrections due to higher order EW corrections. These effects, depending on the calculation scheme and in particular how photons close to the directions of incoming partons and outgoing leptons are treated [49,65], may be sizable or at 1 % level only. Appropriate solution need to be chosen, for the simulation purposes and observables definition. The NLO programs like [40] can be useful; to obtain predictions, and for discussion of systematic errors.