Anomalous coupling, top-mass and parton-shower effects in ${W^+W^-}$ production

We calculate the process $pp\to W^+W^-\to e^+ \nu_e \mu^-\bar{\nu}_\mu$ at NLO QCD, including also effective field theory (EFT) operators mediating the $ggW^+W^-$ interaction, which first occur at dimension eight. We further combine the NLO and EFT matrix elements produced by GoSam with the HERWIG7/MATCHBOX framework, which offers the possibility to study the impact of a parton shower. We assess the effects of the anomalous couplings by comparing them to top-mass effects as well as uncertainties related to variations of the renormalisation, factorisation and hard shower scales.


Introduction
Among the most important goals for the next phase of LHC data taking are precision tests of the electroweak symmetry breaking sector and the search for signs of new physics. In this respect, the final state of W + W − plays a prominent role. For example, the continuum pp → W + W − → llνν is the dominant background in the measurement of H → W + W − → llνν. The process pp → W + W − (+jets) also can be a major background for new-physics processes involving missing energy. Therefore it is very important to have good theoretical control on the pp → W + W − cross section.
Final states with two massive vector bosons recently have attracted additional interest due to the fact that a slight excess at about 2 TeV in the search for di-boson resonances has been reported by both ATLAS [1,2] and CMS [3,4], most pronounced in the hadronic decay channel, which however does not seem to persist in Run II. Further, both the ATLAS and CMS measurements for the W + W − total inclusive cross sections using the leptonic decay channels, at 7 TeV [5,6] and at 8 TeV [7,8], are about 10-20% higher than the NLO predictions obtained from MCFM [9,10], which include the gg-initiated sub-process [11]. However, the latter discrepancy has been largely reduced by the NNLO predictions which became available recently [12][13][14]. In addition, it has been noticed that resummation of large logarithms arising from the jet veto condition needs to be taken into account carefully [15][16][17][18], and that the discrepancy for the fiducial cross section is only at the 1σ level, such that the way the extrapolation from the fiducial cross section is done should be revisited [17]. Considering all these recent developments, the need for precise phenomenological studies, also at the level of differential distributions and in view of possible BSM contributions, is evident.
Let us briefly review the history of higher-order calculations in the W + W − (+jets) channel: The process gg → W + W − has been calculated in continuously improving approximations in the literature: the calculation for on-shell W bosons has been performed in [19,20]. Leptonic decays of the W bosons were included in [21] for massless fermion loops and extended to include the masses of the top and bottom quarks in [22]. Analytic results, including the mass of the top quark, were presented in [10,11], together with a phenomenological study of interference effects with H → W + W − . Focusing on a Higgs-boson mass of about 125 GeV, an update of interference effects has been performed in [23,24] and in [25], where the latter includes higher-order corrections to the interference in a soft-collinear approximation up to NNLO. Very recently, the NNLO corrections to the process pp → W + W − were calculated in [12], removing the discrepancy to the data at 7 TeV, and decreasing the excess at 8 TeV to a level below 1σ. Electroweak corrections to the full 4-lepton final state, including also mass effects, have been calculated in [26]. For a phenomenological study of electroweak and QCD effects see also [27]. A study of combined electroweak [28,29] and QCD corrections (assuming that they factorise), including also matching to the angular-ordered parton shower, has been performed in [30] and is also available in Herwig7 [31,32], the successor of Herwig++ [33].
The NLO QCD corrections to the process qq → W + W − for on-shell W bosons have been calculated in [34,35]. The helicity amplitudes for the process including decays have been calculated in [36], followed by phenomenological studies in [9,37]. Matching with parton showers of these processes has been included in MC@NLO [38]. Weak-boson pair production with NLO QCD corrections, matched to a parton shower with the Powheg method [39,40], has been directly implemented in Herwig++ [33].
The process pp → H → W + W − also has attracted recent interest in view of measuring the Higgs width using information from off-shell production and decay, as proposed in [23,41,42] and further investigated in [43,44]. Such a measurement already has been performed based on the ZZ final state [45].
Calculations of the process pp → W + W − + jet, without including the gg initial state, have been performed in [46][47][48], and recently, including also NLO electroweak corrections, in [49]. The loop-induced process gg → W + W − + jet has been studied in [50]. A very detailed NLO study of 4-lepton plus 0,1-jet final states, including NLO matching to a parton shower and merged samples, H → W W * interference studies and squared quark-loop contributions, has been presented in [51].
In this paper, we calculate the process pp (→ W + W − ) → e + ν e µ −ν µ at NLO QCD, combining the hard matrix elements produced by GoSam [52,53] with the Herwig7/Matchbox [31,32,54] framework, which offers the possibility to study the impact of a parton shower. In addition, we particularly focus on the loop-induced process gg (→ W + W − ) → e + ν e µ −ν µ , where we investigate how new-physics effects which modify the effective gg W + W − coupling could affect various distributions. To this aim we include the most general effective field theory (EFT) operators mediating the gg W + W − interaction, which first occur at dimension eight, in our automated setup. This allows us to assess the impact of these operators in various effective coupling scenarios.
2 Details of the calculation The diagrams contributing to the process gg → e + ν e µ −ν µ (see Fig. 1) comprise of diagrams involving two W propagators ("doubly resonant") as well as diagrams involving only one W propagator ("singly resonant", see Fig. 1b). Note that the latter are important to maintain gauge invariance. Nonresonant diagrams, i.e. diagrams containing no W propagator, do not contribute to the e + ν e µ −ν µ final state, but they would contribute to the e + ν e e −ν e final state, which we are not considering here. We include massive top-and bottom-quark loops, all other quarks (and leptons) in the hard process are assumed to be massless. The photon-exchange graphs vanish due to Furry's theorem.
The Z-exchange diagrams are proportional to (m 2 u − m 2 d )(p 2 3 − p 2 4 ), when summed over up-and down-type contributions. Therefore these diagrams also vanish for massless quarks. For arbitrary invariant masses of the charged lepton -neutrino pairs, i.e. p 2 3 = p 2 4 , and the third generation, where we assume m t = 0 and m b = 0 in the loops, we find that the contributions from doubly-resonant (Fig. 1a) and singly-resonant (Fig. 1b) diagrams with internal Z-propagator cancel each other. The only triangle graphs that contribute are thus the Higgs-exchange diagrams (Fig. 1c) where the amplitude contains the QQH Yukawa couplings. The box diagrams do not involve these couplings and therefore form a gauge-invariant subset. Figure 1: Examples of diagrams contributing to the process gg → e + ν e µ −ν µ .

Operators parametrizing the gg W + W − coupling
The first operators that mediate a four-boson interaction between the two gluons and the two W bosons occur at dimension eight. The gluonic field-strength tensor is defined as the field-strength tensor of the W is defined as We can write the SU(2) fields W I in terms of the physical fields: where Λ denotes the scale below which the EFT description is valid.
Combining the SM Lagrangian with the one including the effective couplings, we have Combining higher-order QCD corrections from L SM with the part containing the higherdimensional operators means that we are performing a simultaneous expansion in α s /2π and in c i /Λ 4 . This requires a careful assessment of the relative importance of the various terms in such an expansion and of the range of validity of the effective theory. We will come back to this issue in Section 3.
As in the SM the gg W + W − coupling is loop induced, we will calculate the following contributions to the gg (→ W + W − ) → e + ν e µ −ν µ cross section, depicted schematically in Fig. 2: Note that the last term above is suppressed by 1/Λ 8 . In the following we list the Feynman rules corresponding to the dimension-eight operators. All momenta are considered to be incoming:

Loop-induced processes in GoSam
The virtual amplitudes for the process pp (→ W + W − ) → e + ν e µ −ν µ , as well as the spin-and colour-correlated tree amplitudes, have been generated with the program GoSam [52,53], which is an automated package to generate one-loop amplitudes. It is based on a Feynman-diagrammatic approach using Qgraf [55] and Form [56,57] for the diagram generation, and Spinney [58] and Form to produce optimised Fortran90 code. For the reduction of the one-loop amplitudes the user can choose between three different reduction libraries (or a combination thereof): Ninja [59][60][61], Golem95 [62][63][64] or Samurai [65,66]. The default setup uses Ninja in combination with Golem95 as a rescue system for numerically problematic phase-space points. The scalar basis integrals have been evaluated using Golem95 and OneLOop [67]. We use the complex mass scheme [68] throughout our calculation.
The calculation of loop-induced processes is straightforward with GoSam, as it is based on a Feynman-diagrammatic approach. However, we have improved the rescue system for loop-induced processes: in the standard case where a tree-level amplitude exists, the accuracy at a given phasespace point is assessed by comparing the value for the coefficient of the single infrared (IR) pole with the exact value (which is obtained from the universal IR behaviour of the subtraction terms). As this coefficient must be zero in the loop-induced case, we have implemented an accuracy check which tests the precision of this zero. In more detail, we have added the following options to the GoSam input card: • PSP_chk_li1: allows to set the desired precision of the pole part (which should be zero) in comparison to the finite part. If the pole part is at least PSP_chk_li1 orders smaller than the finite part, the point is accepted.
• PSP_chk_li2: for loop-induced processes, this option is used instead of PSP_chk_th2. It is the threshold to declare a phase-space point as "bad", based on the precision of the pole in comparison to the finite part. Points with precision less than this threshold are directly reprocessed with the rescue system (if switched on), or declared as unstable. According to the verbosity level set, such points are written to a file and not used when the code is interfaced to an external Monte Carlo program in accordance with the updated BLHA standard [69].
• Similarly, PSP_chk_li3 is used instead of PSP_chk_th3 as threshold for the rotation test in the loop-induced case.
• PSP_chk_li4 sets the minimum pole precision for the points which already have been reprocessed by the rescue system. If a rescued point gives a pole coefficient which is at least PSP_chk_li4 orders smaller than the finite part, the point is accepted.
For the Standard-Model case, we have validated our results by comparing to MCFM [11].

Extended BSM support in GoSam
We have implemented various new features in GoSam which facilitate the calculation of corrections beyond the Standard Model, as well as the interference between SM loop corrections and BSM effects described within an EFT framework. Among the new features are: • the import of BSM model files in UFO format [70] in combination with the updated BLHA standard [69] for the definition of new couplings has been implemented, • BSM-SM interference terms can be calculated by specifying the orders in the corresponding couplings, • sub-process specific settings in the GoSam input card are possible, • the renormalisation term for the Wilson coefficient in the effective ggH coupling has been added, • the coefficients of the effective couplings multiplying the higher-dimensional operators can be modified by the user in the GoSam input card.
For the process considered here we use a model file in UFO format [70] which we create by extending the SM Lagrangian of FeynRules [71] with the EFT operators outlined in Sec. 2.2. The SM parameters, which are provided per default with the UFO model file, are then overwritten in accordance to the updated BLHA standard [69] by the ones listed in Sec. 3. In order to be able to compute the pure SM NLO QCD contributions as well as the additional interference between the SM one-loop contribution and the effective coupling with only one model file, we restrict the one-loop contributions to allow only for SM couplings, with the help of GoSam's diagram filter facilities, while the tree-level contributions are allowed to include the effective coupling.

Interface to Herwig7 and computational setup
Herwig7 features the full simulation of particle collision events up to the particle level, i.e. perturbative as well as non-perturbative physics. The perturbative part provides the simulation of hard processes at NLO QCD (including several built-in LO and NLO matrix elements, LH event file input as well as the fully automated assembly of NLO QCD calculations for almost all Standard-Model processes, utilizing various interfaces to several external matrix-element providers), shower Monte Carlo algorithms, as well as the corresponding LO and NLO matching procedures (dedicated matrix element corrected shower plug-ins and built-in matched cross sections, as well as a fully automated matching machinery).
The capabilities of Herwig7 for integration, unweighting and sub-process parallelization, as well as the steering at the level of input files, are significantly improved. By virtue of the Matchbox framework new integrator modules were introduced, which provide for an efficient, automated multi-channel phase-space sampling: the one which we employ, for the process considered here, is based on the standard sampling algorithm contained in the ExSample library [72].
Based on the Matchbox framework, Herwig7 facilitates the automated setup of all subprocesses and ingredients necessary for a full NLO QCD calculation in the subtraction formalism: an implementation of the dipole subtraction method based on the approach by Catani and Seymour [73], interfaces to various external matrix-element providers, or plug-ins to various in-house calculations for the hard sub-processes or to built-in colour and helicity sub-amplitudes.
In the case where the necessary one-loop and tree-level parts are obtained from external matrixelement providers there exist several possibilities: either at the level of colour-ordered sub-amplitudes or at the level of squared matrix elements through the updated BLHA standard interface [69]. A more detailed description of the interface between Herwig++ and GoSam has been given in [74], for the example of Z+jet production. The interface between Herwig7/Matchbox [32,54] and GoSam-2.0 [53] is fully automated for one-loop QCD corrections, and extensions thereof to handle loop-induced processes and additional contributions from EFT operators are employed for the process considered here.
Fully automated matching algorithms are available, inspired by MC@NLO [38] and Powheg [75] (referred to as subtractive and multiplicative matching respectively), for the systematic and consistent combination of NLO QCD calculations with both shower variants in Herwig7 (facilitating two coherent shower algorithms -an angular-ordered parton shower [76] as well as a dipole shower [77], including the simulation of decays with full spin correlations). For the process studied here, we eventually combine our fixed-order NLO result (supplemented with loop-induced and EFT contributions) with the Herwig7 angular-ordered parton shower [76] through the subtractive (i.e. MC@NLO-like) matching algorithm based on [38,54].

Phenomenological studies
For our phenomenological studies of the process pp → e + ν e µ −ν µ + X we choose a center-of-mass energy of 13 TeV. We also study the behaviour of the BSM effects when going from 8 TeV to 13 TeV.
We use the MMHT2014nlo68c_nf4 [78] parton distribution functions (PDFs), with 4 massless quark flavours 1 in the initial state, and we set α s (M Z = 91.1876 GeV) = 0.12 in accordance with the PDF set we use.
Our default scale choice for the 13 TeV results is a dynamic scale, The mass of the top quark has been set to M t = 174.2 GeV, the Higgs mass to M H = 125.7 GeV. We further use a non-zero top width of Γ t = 1.4 GeV, and a Higgs width of Γ H = 4.11 MeV.
For the electroweak input parameters we follow the SLHA [79] scheme 2 for a set of SM low-scale input parameters to fix the electroweak sector: the electroweak input parameters we choose are M Z = 91.1876 GeV, α = α(M Z ) = 1/128.91 and G F = 1.16637 · 10 −5 GeV −2 , from which M W and sin 2 (θ w ) are subsequently derived. Furthermore we choose Γ Z = 2.4952 GeV and Γ W = 2.085 GeV.
The events are analysed using an in-house analysis for Rivet-2.4 [80] interfaced by Herwig7. The W bosons are directly reconstructed from their leptonic decay products (not via the built-in WFinder function of Rivet).
Our cuts on the analysis level are as follows. To mimic W -identification cuts, we employ a cut on the invariant mass of each same flavour lepton-neutrino pair of 60 GeV ≤ m lν l ≤ 100 GeV. We further require the net transverse momentum of the lepton-neutrino pair to be larger than 10 GeV, and a minimum p T of 25 GeV for each identified lepton and for the missing transverse energy of the event. The identified leptons are required to be in the rapidity range −3 ≤ y l ≤ 3.
For numerical stability we also employ cuts at the generator level, which are of course less restrictive than the cuts employed at the analysis level.
In the following three subsections we will first concentrate on the fixed-order results, and then discuss the impact of a parton shower in Sec. 3.4.

Gluon-induced contributions and effects of higher-dimensional operators
We start the discussion of the results with the gg-initiated processes, i.e. contributions which are either loop induced in the SM, or require dimension-eight operators to contribute at the tree level. In the following we will distinguish three different contributions, see Eq. (2.7). One is the pure Standard-Model loop-induced contribution, where the matrix element is the square of the gg-initiated one-loop amplitude. In the plots this contribution is denoted as gg_SM. The second contribution is the interference term between the SM one-loop amplitude and the EFT tree-level amplitude, which means it is a linear term in the higher-dimensional operators. This contribution is labeled as gg_Interf. Finally the third contribution stems from the square of the EFT tree-level amplitude and is therefore quadratic in the higher-dimensional operators. This contribution is labeled as gg_Eff2. This distinction allows us to separate the term linear in the higher-dimensional operators from the quadratic one and study their behaviour independently.
Consequently gg_SM+Interf and gg_Eff2+Interf denote the combination of gg_SM or gg_Eff2 with gg_Interf respectively, while gg_All finally denotes the combination of all contributions to the gg initial state (cf. Fig. 2 or Eq. (2.7)).
For the numerical results we have set  Let us first focus on the invariant-mass distribution of the W -boson pair, shown in Fig. 3a. The invariant mass is defined via the momenta of the decay products: The most striking feature is the fact that the interference term of the SM loop-induced gg → W + W − amplitude with the amplitude induced by the dimension-eight operators is negative. As this cannot be displayed in a logarithmic plot, we display it with the sign switched, denoted by gg_NegInterf. As expected, the term linear in the dimension-eight operators dominates over the pure EFT contribution (gg_Eff2) which is quadratic in these operators. This is an obvious behaviour as the quadratic terms receive an additional suppression of a factor c i /Λ 4 . However its contribution increases with the center-of-mass energy (here √ŝ = m W W ) and will eventually dominate over the linear term. The point where this happens depends on the setup and in particular on the chosen value for the anomalous coupling constants. In our example this happens at about 500 GeV (where the yellow and purple curves cross). As the linear term is negative, this is related to the point where the sum of the two higher-dimensional contributions (gg_Eff2+Interf) becomes positive, which happens a bit earlier at about 400 GeV (where the green and the dashed red curves cross). While the SM contribution drops rapidly as m W W is increasing, the dimension-eight contributions increase and start to dominate at around m W W ∼ 500 GeV. This is the expected behaviour as the contribution from a dimension-eight operator can increase maximally with s 2 /Λ 4 .  The scale-variation uncertainties (shown as shaded bands in the plots) are relatively large, due to the fact that the results for the gg-initiated subprocess are leading-order predictions.
The fact that the term linear in the dimension-eight operators is negative leaves us with a phenomenologically interesting constellation. In the region where the linear term is dominant we get a decrease in the cross section compared to the SM prediction, and with increasing invariant mass we observe a (partial) cancellation between the linear and the quadratic term. At the point where linear and quadratic term are of the same magnitude, we recover exactly the Standard-Model contribution. This means that putting experimental constraints on these types of couplings will be more difficult, and signs of new physics would be masked: in the low energy region the effects of the dimension-eight operators are anyway suppressed by a factor of 1/Λ 4 , and for larger energies we find (partial) cancellation between the linear and the quadratic term. However, we also emphasize that the sign of the dimension-eight operators is not necessarily fixed to be positive. A negative value would lead to a constructive interference between the linear and the quadratic term rather than to a cancellation. The bounds on negative values will therefore be much more restrictive than on positive values.
The region where the quadratic term becomes equally important and eventually dominates over the linear term has to be interpreted with care. The two terms being equally important means that the suppression of the quadratic terms by the additional factor c i /Λ 4 is compensated by a factor of s 2 . In other words s ∼ Λ 2 , which means that we are probing the scale of New Physics and which is the point where the EFT approach becomes invalid, as it is based on the assumption that it is a low-energy effective theory and that the scale of New Physics is much higher than the scale we are probing. It is this assumption that allows us to be confident that operators of lower dimensions are more important compared to higher-dimensional operators. If higher-dimensional operators were not sufficiently suppressed, it would not be justified to neglect dimension-ten operators, whose linear terms are actually less suppressed than the quadratic term of a dimension-eight operator. And even worse, in the case of s ∼ Λ 2 all higher-dimensional operators could contribute equally and there is no physically meaningful expansion anymore. This point simply denotes the breakdown of the EFT approach. A related issue is the possible violation of unitarity, which we will discuss in Sec. 3.1.1.
In Fig. 3b we display the observable ∆R W W = (y 1 − y 2 ) 2 + (φ 1 − φ 2 ) 2 , the separation in rapidity and azimuthal angle between the two W bosons. Here we see that the effects of the higher-dimensional operators lead to an enhancement of the distribution over the whole range. Note that the region below ∆R W W = π is not populated because we only show the fixed-order results in this subsection.
In Figs. 4a and 4b we show the same observables calculated at a center-of-mass energy of √ s = 8 TeV. We observe that the on-set of the BSM effects in the m W W distribution is around 550 GeV. However, the relative size of the BSM contributions with respect to the SM contribution is much larger at √ s = 13 TeV, as can be seen by comparing Figs. 3b and 4b. We have verified that the qualitative behaviour in comparing 13 to 8 TeV stays the same if we also use a fixed scale (M W ) for the √ s = 13 TeV case.  In Fig. 5 we compare the two scale choices µ r = µ F = M W (fixed) and µ r = µ F = m W W = (p e + + p µ − + p νe + pν µ ) 2 (dynamic), including the scale uncertainty band, obtained as usual by varying by a factor of two up and down from those central scale choices. The fixed scale M W , being relatively low, leads to a larger value of α s and therefore an increase in the cross section. Since the bands do not overlap, this also means that the scale variations by factors of two up and down are not sufficient to capture the uncertainties correctly.
Another interesting distribution is the relative azimuthal angle between the two charged leptons, ∆φ e + µ − , shown in Fig. 6a. The contributions from the higher-dimensional operators lead to more highly boosted W bosons, and therefore the associated leptons are more likely to be "back-to-back". A similar behaviour can be seen in the ∆R distribution of the leptons (see Fig. 6b).
We also show the various contributions to the transverse momentum of the positron from the W + decay in Fig. 7a and the invariant mass of the charged leptons in Fig. 7b. The BSM effects lead to a clear enhancement in the p e ⊥ spectrum, which becomes quite substantial already for p e ⊥ values as low as ∼ 60 -100 GeV. In the m e + µ − distribution, the effect of the higher-dimensional operators is also clearly visible for energies larger than about 150 GeV to 190 GeV already (taking scale-variation uncertainties into account). However, as this concerns only the gg-initiated contribution, the effect will be washed out once all sub-processes contributing to the pp initial state, plus higher-order  Figure 6: Distributions of (a) relative azimuthal angle ∆φ e + µ − (b) and ∆R e + µ − , for the gg-initiated contributions. Ratio plots are with respect to gg_All. The shaded bands show the scale-variation uncertainties.

Herwig7+GoSam
gg All gg SM gg Eff2 gg SM+Interf gg Eff2+Interf gg NegInterf corrections, are taken into account, as will be discussed in Sec. 3.3.
In order to investigate differences between the three dimension-eight operators, we will now consider them one at a time, always setting the coupling constant of the two others to zero respectively. In Figs. 8 and 9 the effects of the individual operators are shown for two observables, the angle between the decay planes of the W bosons, cos (Ψ eν,µν ), and their invariant mass, m W W . For these comparisons we have always set one of the c i coefficients to the value 0.3 and the other two to zero respectively. Looking at the decay planes of the W bosons in Fig. 8, we find that the first two operators, O 1 and O 2 , show the same angular dependence, whereas the angular dependence of the  third operator O 3 , which contains the dual field-strength tensorW I,µν , is different, which is seen to be more prominent for √ s = 13 TeV. To distinguish between the first and the second operator, the invariant-mass distribution is also a suitable observable, as can be seen in Fig. 9a. Here the first operator leads to a stronger decrease of the distribution in the region around m W W ∼ 500 GeV. Therefore, the combination of these two observables in principle allows for a distinction between the three operators. However, it should be noted that this can only be a qualitative discussion, as the impact of the dimension-eight operators strongly depends on the size (and on the sign) of the coupling constants c i . The overall size of the BSM effects for our default choice of the anomalous couplings is discussed in Sec. 3.3, where we combine all sub-processes contributing to the pp initial state.

Unitarity bounds
In the context of higher-dimensional operators it is also important to talk about unitarity. As the effects of these operators grow with increasing center-of-mass energy, they will eventually violate unitarity. For the case of stable W bosons, i.e. for 2 → 2 scattering, a unitarity bound on the total cross section can be derived along the lines of Ref. [81]. In more detail, we can start from Eq. (48) of Ref. [81] (but use total angular momentum J = 0 for the gg-initiated case), where the bound for an inelastic 2 → 2 scattering amplitude T in , summed over the final state helicities λ 3 , λ 4 , is given by To obtain the bound for the total cross section, we include the flux factor 1/(2ŝ), average over initial state colours and helicities, and sum over the colour and helicity configurations contributing to the cross section. This leads to where we have used λ1,λ2∈{+,−} = 4. The sum over the colour states contributing to the amplitude is given by δ ab δ ab = N 2 c − 1 (the trace of the identity matrix in the adjoint representation). Therefore, with N c = 3, we have σ ggW W ≤ π 2ŝ . (3.6) The derivation of the unitarity bound on the total cross section in Eq. (3.6) is based on a 2 → 2 scattering process. To obtain an estimate for the unitarity bound including the decay of the W bosons one could integrate the 2 → 2 process numerically and rescale the result with the branching ratios of the two W bosons decaying into leptons. This procedure, however, does not take into account the effect of the cuts on the leptons, and therefore we refrain from showing a unitarity bound in the plots for the distributions.
Unitarity arguments can also be employed to calculate an upper bound on the absolute value of the anomalous coupling constants. To do so we use the ansatz to require unitarity of the amplitude for a given set of helicities and project the amplitude onto partial waves [82]. Looking at the scattering a + b → c + d with the corresponding helicities λ a , ..., λ d , the partial wave decomposition reads with λ = λ a − λ b and µ = λ c − λ d , and where θφλ c λ d |T (E)|00λ a λ b denotes the transition matrix element. Its unitarity must hold for each partial wave independently, i.e.
Therefore we project the full amplitude onto single partial waves, where the strongest constraints typically come from the lowest order partial waves. In the case where λ = µ = 0, the d J functions reduce to the Legendre polynomials, i.e. d J 00 (θ) = P J (cos θ). Usually it is assumed that the strongest constraints stem from longitudinally polarized W bosons, as in the limit of large momentum k the longitudinal polarization vector behaves like Projecting onto the 0th partial wave we find (3.10) For the third operator, O 3 , the contribution vanishes for longitudinally polarized W bosons. It is also interesting to note that the contributions from the first two operators increase more mildly with energy than naively expected. For dimensional reasons the denominator in Eq. (3.10) could be ∼ s 2 , which in turn would mean that the amplitude itself could be ∼ s 2 . However we find only a behaviour which grows like ∼ s.
The fact that the third operator vanishes for longitudinal polarizations, and that we do not find the strongest possible increase with the center-of-mass energy, suggests to also consider the situation where the W bosons are transversely polarized. Projecting these amplitudes onto the 0th partial wave we find This means that the strongest constraints for energies above the weak scale come from transversely polarized W bosons and, in order to maintain unitarity, one can roughly assume (3.12)

Impact of heavy-quark loop contributions
We have taken both bottom-and top-quark masses into account for the quark loops mediating the SM gg → W + W − interaction. The effect of the bottom-quark mass is negligible (of the order of the Monte Carlo integration error), while top-quark mass effects have a considerable impact on the partonic cross section in the gg-initiated channel. Fig. 10 shows the effect of massive top-quark loops on the invariant-mass distribution of the charged leptons and on the R-separation between the charged leptons. We observe that top-mass effects decrease the m e + µ − distribution by more than 30% below values of m e + µ − ∼ 250 -300 GeV. The SM contribution with M t set to zero (gg_SM_massless) is of the same size as the SM+BSM result with masses taken into account (gg_All) at m e + µ − ∼ 200 GeV, which means that neglecting the top-quark mass in the SM calculation could potentially "fake" BSM effects.
It should be noted here that in the case of massless top quarks also the Yukawa coupling between the Higgs boson and the top quark vanishes. This result is contrasted to a calculation where the top-quark loops have been dropped altogether 3 , shown in Fig. 11. This has a considerable impact on the m e + µ − distribution beyond about 150 GeV, however, the effect is much less pronounced than in the case where top-quark loops are taken into account but the top-quark mass is neglected (cf. Fig. 10). Even though the mass effects are below the 10% level once the full pp initial state including the NLO corrections is taken into account (see Fig. 12), this study demonstrates that massive top-quark loops should be taken into account to describe measurements of highly boosted W bosons.  Figure 12: Impact of the top-quark loops on (a) the invariant mass of the charged leptons and (b) ∆R e + µ − , considering the full pp initial state. Ratio plots are with respect to pp+gg_SM_massless.

Combination of gluon-and quark-initiated channels
In this section we compare the gg-initiated contribution to the full process pp (→ W + W − ) → e + ν e µ −ν µ at NLO QCD, considering first the results at the fixed-order level. Shower effects will be discussed in Sec. 3.4. The important questions here are how visible the effects of the anomalous couplings are if the quark-initiated Standard-Model contributions are added, and to what extent the effects of the higher-dimensional operators are hidden in the theoretical uncertainties.
In Fig. 13 we show the invariant-mass distribution of the W -boson pair including all SM as well as EFT contributions. The effects of scale variations are plotted as well, where the scale uncertainty bands have been obtained by varying by a factor of two up and down from the dynamic scale choice µ R = µ F = m W W . We show the SM NLO contribution with and without the EFT contributions, and in comparison to that the effects of the higher-dimensional operators in the gg-initiated contributions alone. This allows to directly assess the impact of the anomalous couplings. The loop-induced, gg-initiated Standard-Model contribution leads to an O(10%) increase over the quark-or quark-gluon-initiated NLO result. We therefore observe that in the full result, combining all channels, a visible deviation from the SM prediction begins to show at larger m W W values, of about 700 GeV, while in the gg-initiated contribution, shown in Fig. 3a, the deviation already starts to be visible at about 500 GeV to 600 GeV (taking scale-variation uncertainties into account). While for values of m W W around 700 GeV the size of the BSM effects is comparable to the size of the scale uncertainties, shown in Fig. 13, for values of m W W of about 800 -900 GeV, the deviations from the Standard Model due to the higher-dimensional operators start to become clearly visible. On the other hand, the region beyond 1 TeV already probes energies where the EFT approach starts to become invalid.  Figure 13: Distributions for (a) invariant mass m W W and (b) separation ∆R W W of the reconstructed W -boson pair in pp (→ W + W − ) → e + ν e µ −ν µ for the sum of all partonic channels, including the effects of the higher-dimensional operators. pp+gg_SM includes all partonic channels, i.e. all the quark-initiated channels up to NLO QCD plus the loop-induced SM gg-initiated contribution gg_SM. pp+gg_BSM is the same but includes the loop-induced SM+BSM gg-initiated contribution gg_All instead of just gg_SM. In addition we show the SM and SM+BSM gg-initiated contributions separately. The shaded bands show the scale-variation uncertainties. Ratio plots are with respect to pp+gg_SM.  a parton shower. We therefore combine our fixed-order NLO results (supplemented with loopinduced and EFT contributions) with the Herwig7 angular-ordered parton shower [76] through the subtractive (i.e. MC@NLO-like) matching algorithm based on [38,54]. 4 Uncertainties in the shower are mainly quantified by varying the hard shower scale µ Q which provides a reliable estimate of missing logarithmic orders as well as the impact of large-angle, hard and thus unreliably modelled emissions. It also serves as a check to verify the improvements expected from NLO plus parton-shower matching. While, as expected, the invariant mass of the reconstructed W + W − system is not affected by the parton shower, a number of other infrared sensitive observables receive large contributions, both through the NLO real radiation and further subsequent parton-shower emissions. Typical infrared-sensitive distributions in this case are the R-separation of the two W bosons (where in the zero-jet limit ∆R is composed solely of a difference in rapidity, while ∆φ = π), as well as the transverse momentum of the reconstructed W + W − system, shown in Fig. 14. Both observables show the expected behaviour with respect to additional radiation; in the region ∆R < π, both the NLO real emission as well as shower emissions off the gg-induced channel contribute. The first contribution includes a small shower uncertainty, as this kinematic range has been improved by the NLO matching. Once the BSM contribution to the gg-channel becomes dominant, pure shower emissions off this sub-process become more important and hence yield a larger uncertainty. Ultimately, NLO QCD corrections, or at least a leading-order multi-jet merging are desirable in this case. Similar features are present in the transverse momentum of the reconstructed W + W − pair. Azimuthal and R-separations of the charged leptons are sensitive to the BSM contribution and very stable against QCD activity, as shown in Fig. 15.

Parton-shower effects
We finally discuss a few observables which are relevant to the experimental reconstruction of the W + W − final state, particularly lepton-jet separations and the distribution of missing transverse momentum, displayed in Fig. 16. While shower uncertainties at the level of 10% are observed, the lepton-jet separation is rather stable against QCD activity, and BSM contributions only affect the normalization in the small-∆R region; the experimentally required lepton-jet isolation is thus not introducing any bias. Larger impact is observed on the transverse momenta of the charged leptons (e.g. as shown for the positron in Fig. 16), which, however, turn out to be rather stable with respect to parton-shower scale variations.

Conclusions and outlook
The production of electroweak gauge-boson pairs is amongst the most important signatures at the LHC. These final states are important Higgs-boson decay channels, and they allow us to study the electroweak sector, with the aim to reveal the mechanism of electroweak symmetry breaking. We have studied the production of a pair of W bosons at NLO QCD, in the light of additional anomalous couplings. In particular, we have also included the gg-initiated (loop-induced) process gg (→ W + W − ) → e + ν e µ −ν µ , which is formally a contribution to the NNLO result, but is enhanced due to the large gluon luminosity at the LHC. In addition to the Standard-Model gg-initiated contribution, we have included gg-initiated contributions stemming from dimension-eight operators which induce a tree-level coupling between gluons and electroweak gauge bosons. This possibility has not been discussed in the literature so far. Their presence leads to an interference between the Standard-Model gg-induced one-loop amplitude and a tree-level amplitude mediated by dimension-eight operators.
We have discussed their effects on a variety of important observables. We have found that the presence of dimension-eight operators can lead to substantial effects in the high-energy tail of the distributions, which can be used by the LHC experiments to constrain the parameter space for the associated effective couplings.
Furthermore, we have investigated the importance of heavy (SM) quarks in the loop-induced process, leading to corrections of up to 10%, depending on cuts and center-of-mass energy.
Finally, by combining our fixed-order results with the Herwig7 angular-ordered parton shower, we have studied the effects of a parton shower, including variations of the hard shower scale, on the leptonic observables and on observables related to the reconstructed W + W − system.