The Effective Vector Boson Approximation in high-energy muon collisions

Due to the inclination for forward gauge radiation, lepton colliders beyond a few TeV are effectively electroweak (EW) boson colliders, suggesting the treatment of EW bosons as constituents of high-energy leptons. In the context of a muon collider, we revisit the validity of W and Z parton distribution functions (PDFs) at leading order in 2 → n process. We systematically investigate universal and quasi-universal power-law and logarithmic corrections that arise when deriving (polarized) weak boson PDFs in the collinear limit. We go on to survey a multitude of 2 → n processes at s\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \sqrt{s} $$\end{document} = 2–30 TeV via polarized and unpolarized EW boson fusion/scattering. To conduct this study, we report a public implementation of the Effective W/Z and Weizsäcker-Williams Approximations, which we collectively call the Effective Vector Boson Approximation, into the Monte Carlo event generator MadGraph5_aMC@NLO. This implementation lays the groundwork for developing matrix-element matching prescriptions involving EW parton showers and renormalized EW PDFs. To further with this agenda, we give recommendations on using W/Z PDFs.

theorems and Sudakov exponentiation in QED and pQCD, which at times rely more on the presence of multiple, well-separated (hierarchical) mass scales than on being unbroken gauge theories [62][63][64][65][66][67][68][69][70]. Clearly, being an Abelian/non-Abelian or weakly/strongly coupled theory is less crucial for sufficiently inclusive processes. At the same time, differences between collinear factorization in pQCD and the EW theory must exist since lepton and hadron beams are not composed of weak isospin-averaged states.
More specifically, the fact that muons carry EW quantum numbers implies that their collisions do not represent an inclusive summation over all initial-state weak isospin charges. (This would require µ − ν µ and ν µ − ν µ beams.) As a result, infrared logarithms beyond lowest order in perturbation theory do not fully cancel, leading to violations of the Block-Nordsieck Theorem [4,[71][72][73][74][75][76]. The analogy in pQCD is the violation of the Collinear Factorization Theorem at three-loops when applied to exclusive hadronic final states, e.g., dijet production [70,77,78]. However, despite this violation, application of the Collinear Factorization Theorem, which is presently only proved for a handful of processes [64,65,70], to arbitrary processes remains a quantitatively successful paradigm. Motivated by this success, we consider whether such a paradigm can also work for high-energy lepton collisions.
As a step to better understanding collinear factorization in the EW sector and to further explore the ability of the EWA to predict total and differential cross sections, we consider a framework that combines the EWA for helicity-polarized W and Z bosons with the Weizsäcker-Williams Approximation for helicity-polarized photons. We collectively label this the Effective Vector Boson Approximation (EVA). † In this framework and in the context of a multi-TeV µ + µ − collider, we investigate the impact of and validity of (helicity-polarized) γ/W/Z PDFs in 2 → n process. To focus on the role of partonic kinematics, we restrict ourselves to leading order (LO) matrix elements and bare, i.e., unrenormalized γ/W/Z PDFs, which are finite at LO. Processes that we consider include: associated and many-Higgs production, many-boson production, as well as associated and multi-top production. We extend recent studies [22,26,42,79,80] by investigating universal and quasi-universal corrections to weak boson PDFs that appear naturally in their derivations. Specifically, we study universal power corrections of the form (p 2 T /Q 2 ), which spoil the accuracy of collinear factorization, and quasi-universal power corrections of the form (M 2 V /Q 2 ), which spoil the accuracy of the Goldstone Equivalence Theorem [81,82]. Importantly, we also consider the role of universal and quasi-universal logarithmic corrections of the form δσ/σ ∼ O[log(µ 2 f /M 2 V )], by exploring scale variation and when the evolution variable in weak boson PDFs is defined in terms of transverse momentum (p T ) or virtuality (q). This is in addition to studying the role of helicity in both total and differential cross sections. We note that this study is complementary to extensive studies on uncertainties of the EWA [79,80,83].
We find strong sensitivity to power corrections when hard-scattering scales Q are below Q ∼ 1 TeV; for larger Q, we report agreement between full and approximated MEs when scale uncertainty bands, which can be large, are taken into account. More explicitly: we find that computations with the EWA can reproduce total and differential results within (large) scale uncertainties, so long as factorization-breaking power corrections are sufficiently suppressed. Even for asymptotically large energies, we find scale uncertainties remain large, demonstrating a need for renormalization group (RG) evolution in our factorization theorem for high-energy muon collisions. To strengthen the parallels with pQCD, we give a proof-of-principle demonstration of matrix element matching with transverse weak boson PDFs and full MEs. Given these criteria, we go on to survey nearly two dozen 2 → n VBF processes with the EVA in µ + µ − collisions at √ s = 2 − 30 TeV. Cross sections and their scale uncertainties are presented for both helicity-polarized and unpolarized initial states. To conduct this study, we report an implementation of the EVA into the generalpurpose Monte Carlo event generator MadGraph5_aMC@NLO (mg5amc). Notably, this public implementation lays the groundwork for developing QCD-like matching prescriptions with initial/final state EW boson radiation as well as (polarized and unpolarized) PDFs that are RG evolved via the EW theory and pQCD. To further with this agenda, we also give some recommendations on using weak boson PDFs in high-energy lepton collisions. The remainder of this work continues in the following order: In Sec. 2, we summarize the EVA formalism used throughout this work and present a formula for EW boson scattering in high-energy µ + µ − collisions. In Sec. 3 we document our computational setup and numerical values for SM inputs. Sec. 4 is the first of two principle sections and where we revisit the validity of the EWA. Sec. 5 is the second of two principle sections and where we give a survey of 2 → n VBF processes in the EVA. We conclude in Sec. 6. There, we give an extended discussion of our findings, reflecting particularly on the parallels we find with more subtle aspects of PDFs in QCD, e.g., the phenomenon of "precocious scaling." Finally, we provide some recommendations on using weak PDFs in high-energy lepton collisions in Sec. 6.1. App. A provides some instructions for reproducing our results and using (un)polarized EW boson PDFs in mg5amc. For completeness, a derivation of helicity-polarized EW boson PDFs at LO is given in App. B.

The Effective Vector Boson Approximation for µ + µ − collisions
In this section, we summarize the EVA, i.e., the framework in which we work, and its use in evaluating scattering cross sections in many-TeV µ + µ − collisions. While we focus on muons, the EVA is, in principle, applicable to any lepton-lepton, lepton-hadron, and hadron-hadron collider configuration. Extension to other colliders, however, may require substitutions of gauge coupling charges and/or convolutions with additional PDFs [5]. In Sec. 2.1, we state a scattering formula that will be the basis for all our numerical results and validation checks. In Sec. 2.2, we list the q 2 and p 2 T -evolved collinear PDFs that describe the density of EW bosons in muons at LO. Finally, we document for the completeness in Sec. 2.3 the PDFs for SM neutrinos from muons.

A scattering formula for µ + µ − collisions
To described the fully differential production of an n-body, final state F with momenta {p f } via the high-energy VBF process V λ A V λ B → F, where V λ A and V λ B are helicity-polarized EW gauge bosons, in µ + µ − collisions at a center-of-mass (c.m.) energy of √ s, we invoke the EVA. In practice, this means working from a scattering formula given by σ(µ + µ − → F + X) =f ⊗f ⊗σ + Power and Logarithmic Corrections (2.1) Here, σ is the muon-level (beam-level) inclusive cross section for the production of F in association with an arbitrary state X. Explicitly, X consists of at least two leptons l, where l = µ ± , ν µ , or ν µ , in addition to particles originating from radiative corrections. The summation runs over all polarized EW boson V λ ∈ {W ± λ , Z λ , γ λ }, with λ ∈ {0, ±1}. Formally speaking, when the collection of states {V λ } is extended to left-handed (LH) and right-handed (RH) states ν µL and (ν µ ) R , the beam remnant X includes weak bosons.
For beams k =1,2, the quantitiesf V λ /µ ± (ξ k , µ f ) are the bare PDFs that describe the likelihood that an unpolarized muon µ ± with energy E µ = √ s/2 and momentum p µ = E µ (1, 0, 0, ±1) contains a "parton" V with helicity λ, mass M V , energy E V = ξ k E µ , and no transverse momentum p T,V λ . Following Ref. [84], we adopt thef notation to stress that the PDFs in Eq. (2.2) are not resummed. Thef are related to resummed PDFs f by Generally, E V = E V in the frame of the (V V )-system since generally M V = M V . Iñ f V λ /µ ± , the quantity µ f is the collinear factorization scale and acts as the ultraviolet regulator of the bare PDF. Physically, µ f is the phase space upper bound on the norm of the space-like momentum transfer q = (p µ − p l ) carried by V λ (q); alternatively, µ f can be interpreted as the upper bound on the p T of lepton l in µ ± → V λ + l splitting. The (phase space) integrals over the momentum fractions ξ k are bounded by the (dimensionless) kinematic threshold variable τ 0 = min(M 2 V V /s) = min(M (F)/s). For M V V < √ τ 0 s, the (V V )-system has insufficient energy to produce the n-body state F. The separately Lorentz-invariant phase space measure dP S n is given by the usual expression In Eq. (2.2), dσ/dP S n is the totally differential, "parton-level" cross section for the hard-scattering process V λ A (p A ) V λ B (p B ) → F({p f }), which occurs at a hard scale Due to this equality, we use the terms "hard-scattering system" and "(V V )-system" interchangeably. The helicity-polarized cross sections can be computed from the formula [85] dσ Here, λ(x, y, z) = (x − y − z) 2 − 4yz, is the usual Källen kinematic function that accounts for the masses of initial-state particles. Unlike traditional leading-twist approximations that neglect masses of initial-state partons, the ME M(V λ A V λ B → F) is evaluated with nonzero M V , M V . (In none of our results are weak boson masses set to zero.) Moreover, unlike the scattering of unpolarized partons in unpolarized beams, no spin-averaging factor for initial-states V λ A V λ B is needed for helicity-polarized cross sections that are paired with PDFs for helicity-polarized partons. The summation in Eq. (2.6) runs over all discrete degrees of freedom (dof) related to F, e.g., electric charge and color helicity multiplicities. Importantly, if the summations V λ A ,V λ B and dof do not run over all helicity polarizations for V V and F, respectively, then the square of M is not Lorentz invariant. In such cases, an infrared-safe reference frame must be specified to define the helicities. For further details on evaluating helicity-polarized cross sections, particularly in relation to PDFs for polarized partons and polarized parton showers, see Ref. [85]. Implicit in Eq. (2.2) is a restriction on the phase space integration measure dP S n . The purpose of this restriction is regulate M and render it meaningful. For example: the ME for the process γγ → qq, where q is a massless quark, diverges without phase space cuts on t-channel momenta. Cuts should also ensure that s-, t-, and u-channel invariants in the V λ A V λ B → F hard process are comparable to one another and to the hard scale M V V . In principle, this means that logarithms (L) of ratios of these invariants, which appear in hard-scattering cross sections, are never numerically large. In practice, we vary phase space cuts to explore the growth of these logarithms in initial-state µ ± → V λ l splitting (see Sec. 4.5). While the MEs for all the processes that we investigate are regulated, we set looser phase space cuts on final-state kinematics to balance computational demands. For some processes, logarithm can grow as large as O(5 − 10), and therefore remain within perturbative limits in the sense α W × L 1. We have checked (see Sec. 4.5) that tighter phase space cuts do not qualitatively change our findings.
When deriving Eq. (2.2), a number of assumptions are made. Two important ones are both related to enforcing large separations of scales in V λ A V λ B scattering. The first is that weak bosons are massive but that the invariant mass of the , remain nonvanishing when contracted with µ → l currents. We reiterate that including initial-state parton masses here differs from typical treatments of QCD partons in hadron collisions, which are assumed massless in the absence of specialized schemes [12][13][14][15]. Outside this limit, Eq. (2.2) receives quasi-universal power corrections of the form δσ ∼ (M 2k V /M 2k+2 V V ) k for k > 1, the size of which are quantified in Sec. 4.3. The qualifier "quasi-universal" refers to the fact that such corrections originate from the derivation off V ± 0 PDFs, and therefore appear for any V λ A V λ V scattering process with at least one longitudinally polarized W ± 0 or Z 0 . (Specifically, they come from expanding the ME for µ → l + V 0 splitting.) It is worth noting that the Goldstone Equivalence Theorem requires that these terms be small [81,82]; for further insights on relationship, see Refs. [4,80,86,87].
The second important assumption is the stipulation that EW bosons are emitted at shallow angles in µ → l + V λ splittings, i.e., p T,l ∼ p T,V λ M V V . This is a standard but necessary condition for collinear factorization in gauge theories [70,88]. As in QCD computations, universal power corrections of the form δσ ∼ (p 2k T,l /M 2k+2 V V ) for k > 1 can be incorporated by higher-order perturbative computations, e.g., next-to-leading order (NLO) in α W or α, parton showers, or ME matching to higher leg multiplicities (see Sec. 4.6). To be explicit, "universal" here refers to the fact that such corrections originate from the derivation of bothf V ± 0 andf V ± T , meaning that they are present for any V λ A V λ V scattering process. (Specifically, they come from expanding the ME for µ → l + V λ splitting.) In its present form, Eq. (2.2) is subject to universal and quasi-universal logarithmic weak coupling constant and µ f has the physical interpretation as described above Eq. (2.4). Naïvely, one may argue that these corrections are sub-leading since they are coupling suppressed. However, g W is not a small number and collinear logarithms can compensate for this. For instance: taking µ f = 1−10 TeV implies corrections of δσ ∼ (g 2 W /4π) log(µ 2 f /M 2 W ) that are O(20% − 30%). While we ultimately report in Sec. 4 a prescription for obtaining agreement between full and EWA-based calculations, the uncertainties associated with choosing µ f reported there and in Sec. 5 undercut our findings.
Since Eq. (2.2) is only a LO expression, and therefore does not resum any logarithms, the only (quasi-)universal logarithms that we study are those coming from the µ ± → V T l splittings themselves. For precision computations, an RG-improved version EVA with renormalized PDFs f V ± λ , running couplings, and an EW Sudakov form factor are necessary. Equation (2.2) is written such that renormalized PDFs can be incorporated by the replacement: is acting as a phase space cutoff and the RGE scale.) Were we to replacef V ± λ with their renormalized versions, then the absence of a Sudakov factor still implies that the scattering formula is not scale invariant in an RG evolution sense. That is to say, the anomalous dimensions associated with f V λ A /µ ± and f V λ B /µ ∓ do not necessarily cancel those associated with a renormalized partonic cross sectionσ R (V λ A V λ B → F). Sudakov factors can be incorporated following the classic treatment of Ref. [69] or modern treatments like Ref. [89]. However, investigating and quantifying the impact of these improvements as well as those related to γ T /Z T mixing [4,26,56,[90][91][92][93], which we also neglect, is left to future work.

q 2 and p 2 T -evolved collinear EW PDFs
The expressions for EW boson PDFsf V λ A /µ ± depend strongly on their precise formulation; compare for example Refs. [5-7, 20, 94-96]. As discussed in Sec. 4.5, seemingly innocuous conceptual differences can lead to substantial numerical differences in real computations. Therefore, we now summarize the PDFs used in this study.
In PDFs for W λ , Z λ , γ λ bosons from high-energy muons, one has the freedom to parameterize the momentum transfer in µ(p µ ) → l(p l ) + V λ (q) splittings either by the squared virtuality q 2 = (p µ − p l ) 2 < 0 propagated by V λ , or by the squared transverse momentum p 2 T carried away by l. While the two quantities are related by q 2 (1 − ξ) = −p 2 T , where ξ = E V /E µ is the fraction of µ's energy held by V λ , the resulting PDF sets for transversely polarized V λ differ analytically. Consequentially, for fixed z, λ, and µ f , one can obtain large differences due a relative contribution that scale as δf V T /µ ± ∼ log(1 − ξ). This logarithm diverges in the large-ξ limit and corresponds to a nonzero q 2 but a vanishing p 2 T . Such differences have been sporadically discussed throughout the literature [18-20, 23, 24, 79, 95, 97] but not systematically compared. In light of this, we investigate both sets of PDFs.
For the couplings in Table 1, and assuming q 2 evolution, the LO PDFs for polarized and RH (f R ) fermions in the hard scattering frame arẽ Choosing instead to integrate over p 2 T leads analogously to the following PDFs for V λ : To obtain the LO PDF for γ λ from polarized muons in either evolution scheme, one must make the replacement M V → m µ in thef V T PDFs and neglect thef V 0 PDF. Given a scheme, we construct polarized EW boson PDFs for unpolarized muon beams, denoted bỹ f V λ /µ ± , from those PDFs for polarized muons, denote byf V λ /µ ± λ , through the relatioñ As a technical note, both schemes are available in mg5amc (see App. A for details) but stress that RG evolution of EW boson PDFs from leptons is not yet supported. Differences between the two sets of PDFs appear only in the collinear logarithms for transversely polarized V λ . In this sense, the impact of log(1 − ξ) corrections is process dependent and thus is labeled "quasi-universal." The absence of scale evolution in PDFs for longitudinally polarized V λ is well-known and implies that traditional means of estimating scale uncertainty in pQCD, e.g., three-point scale variation, are not applicable to longitudinally polarized weak boson PDFs. In principal, one can obtainf V T /f L/R fromh V T /f L/R , or vice versa, with appropriately chosen µ f . To further highlight the parallels with pQCD, we note that absorbing factors of (1 − ξ) into factorization scales is common practice in Soft-Collinear Effective Field Theory (SCET) [98,99]. We reiterate that the PDFs here are only accurate to LO. This means that charge-flipping splittings such as µ − → γ * → µ + and µ − → γ * → W + , which appear first at NLO, are neglected.

Collinear PDFs for SM neutrinos
We briefly note that the derivation of W ± λ PDFs in µ → l + V λ splitting also implies the existence of neutrino PDFs. As we are working in the SM, only massless, LH neutrinos (and RH antineutrinos) exist. Therefore, by probability conservation, the µ − L → ν µL PDF at leading order accuracy when evolved by q 2 and p 2 T arẽ As we are interested in VBF at µ + µ − colliders, we do consider further the role of neutrino PDFs from muon beams; for recent discussion on these PDFs, see Ref. [26]. Moreover, while we have also implemented these PDFs into the public release of mg5amc, access to them is temporarily restricted due to the unregulated divergence at ξ → 1. Likewise, throughout this study, we neglect the importance of µ → µ PDFs due to the complication of soft/collinear photon emissions, which necessitates resummation [100]; we refer readers to studies by Refs. [56,[100][101][102][103][104].

Computational Setup
In this section we summarize the computational framework used in this study. Here, we only document the Monte Carlo (MC) tool chain and its tuning. Details on the EVA itself and usage in mg5amc are documented in Sec. 2 and App. A.
To simulate high-p T muon collisions, we employ a development release of version 3.3.0 of MadGraph5_aMC@NLO [105,106]. In this software suite, fully differential events are obtained from tree-level ME that are constructed [107] and evaluated [108] using helicity amplitudes defined in the HELAS basis [109], with QCD color algebra decomposed according to color flow [110]. Helicity-polarized ME are obtained by truncating spin-averaging over initialstate states and/or spin-summing over final-state states [85]. Analysis of parton-level events is handled by MadAnalysis5 [111,112].

Standard Model Inputs
For all ME and PDFs, we take the following EW inputs and masses [113] M W = 80.419 GeV, M Z = 91.188 GeV, G F = 1.16639 × 10 −5 GeV 2 , This implies a QED coupling of α −1 QED (µ r = M Z ) ≈ 132.507. While we consistently modify EW couplings EW inputs are varied but we do not RG-evolve them. Importantly, we have structured mg5amc such that EW couplings and masses present in EW boson PDFs are set to those values stipulated in the param_card.dat configuration file. Changes to EW inputs in this file are automatically propagated into EW boson PDFs. We reiterate that all ME and PDFs assume non-zero W and Z boson masses. We use the light lepton masses m e = 510.9989461 × 10 −6 GeV and m µ = 105.6583745 × 10 −3 GeV for the collinear logarithms contained in the γ λ PDFs. These masses are hard-coded into the γ λ PDFs and are independent of param_card.dat. While it is technically possible use massless e ± /µ ± in MEs, in this paper we choose to use massive leptons.

EVA at high energies
A chief goal of factorization is to simplify in a systematic manner complicated, multiscale MEs that describe many-body processes into a set of simpler, 1-to-2-scale MEs. In practice, this divide-and-conquer approach improves the efficiency and stability of numerical computations. Importantly, the formal perturbative accuracy of factorized calculations can also be improved through quasi-universal RGE methods, e.g., Sudakov resummation and DGLAP evolution. For the specific case of VBF in multi-TeV µ + µ − collisions, factorizing collinear µ → V λ l splittings into weak boson PDFs enables one to reorganize computations of an inherent 3-scale, 2 → (n + 2) scattering process (the three scales being M V , p l T and M V V ) into the product of two 2-scale computations (M V with µ f ∼ p l T , and M V with M V V ) involving process-independent PDFs and process-dependent 2 → n MEs.
As described in Sec. 2.1, the EWA is accurate up to universal and quasi-universal power corrections of the order O(p l2 , which originate from expanding the ME for transversely and longitudinally polarized weak bosons in µ → V λ l splittings, as well as universal and quasi-universal logarithmic corrections of the order O[log(µ 2 f /M 2 V )], which stem from working at LO in the EW theory. In principle, both classes of corrections can be reduced via standard techniques, e.g., higher-order perturbative calculations and Sudakov resummation. In the absence of such improvements, however, there exist theoretical uncertainties in the formulation of weak PDFs that we now explore. In Sec. 4.1 we describe our common setup to study power-law and logarithmic corrections. In Sec. 4.2 we describe how we quantify uncertainties associated with the cutoff scale µ f . We then study e + µ − → HHν e ν µ 1.26 4.43 9.60 e + µ − → γγγν e ν µ 248 · 10 −3 558 · 10 −3 4.04 · 10 −1

Process choice and polarization decomposition
To quantify uncertainties that stem from factorizing polarized EW bosons from initial-state µ ± → V λ l emissions into PDFs, we chose the two benchmark processes e + µ − → HHν e ν µ and e + µ − → γγγν e ν µ . (4.1) Following Ref. [42], we work with e + µ − collisions in order to remove s-channel, µ + µ − annihilation diagrams in a gauge-invariant manner. As such channels have sizable contributions to inclusive cross sections, their removal helps isolate the VBF sub-processes. Under the EWA, these beam-level processes correspond to the partonic processes In practice, we restrict ourselves throughout this section to the EWA helicity configurations We consider these specific processes and configurations due to the high purity of helicity polarizations that drive them. For polarizations λ A and λ B defined in the (W + W − ) frame, we find by explicit calculation [85] that 97% − 99% of HH production in the EWA is dominated by longitudinally polarized W + W − scattering, i.e., (λ A , λ B ) = (0, 0), for √ s = 4 − 30 TeV. In contrast, γγγ production is driven at the 97% − 99% level, albeit with a large scale uncertainty, by transversely polarized W + W − scattering, i.e., (λ A , λ B ) = (T, T ), where T, T = ±1, when assuming the following fiducial phase space cuts on photons p γ T > 50 GeV, |η γ | < 3, and ∆R(γ, γ) > 0. 4. (4.4) In making this distinction between (0, 0) and (T, T ) configurations, we can showcase possible differences of the EWA as applied to longitudinal and transverse polarizations. Many other processes, such as heavy Higgs production and top quark pair production, receive comparable contributions from multiple polarization configurations, which we believe can lead to ambiguities in interpreting the following comparisons.

Defining scale uncertainties for unrenormalized W T /Z T PDFs
A key difference between the bare, LO PDFs in Sec. 2.2 and their renormalized variants is the definition of µ f . For renormalized PDFs, µ f is the RGE scale generated through dimensional regularization; varying µ f is a standard procedure for quantifying perturbative uncertainties in QCD predictions. In the present case of µ → V λ l splitting at LO, µ f is literally a boundary on a phase space integral over either the virtuality |q 2 | of V λ , if one uses Eq. (2.7), or the transverse momentum p l T of l, if one uses Eq. (2.8). For the PDFs of Sec. 2.2, setting µ f proportional to the (V V ) scattering scale M V V is a natural choice as this attempts to captures the whole phase space in V V scattering [5,94]. However, much smaller choices are also favored. Integrating up to p l T ∼ µ f ∝ M V V suggests a potential breakdown of the collinear approximation since one assumes p l T = p V T M V V . As discussed in Sec. 4.6, it is the wide-angle contribution of µ → V λ l splitting that coincides with the regime p l T = p V T ∼ M V V . Therefore, there is an ambiguity, or uncertainty, in the choice of µ f , and increasing or lowering µ f corresponds to conjecturing how much phase space is actually captured by collinear kinematics. It is not guaranteed that arguments used in hadron collisions to fix µ f , e.g., Refs. [114,115], are applicable here. Furthermore, this uncertainty is only one part of the possible uncertainties of the EWA, as illustrated in the factorization formula of Eq. (2.2). Exploring how such ambiguities relate to disagreements between the EWA and full matrix element computations is a reason for this study.
For the γγγ process, we focus on PDF evolution by virtuality (q 2 ) and set the baseline collinear factorization scale to be half the partonic c.m. energy, given by Three-point scale uncertainties for W ± T PDFs are obtained by varying ζ discretely over the range ζ ∈ {0.5, 1.0, 2.0}. While our inspiration to use this procedure draws from common practices in QCD, we reiterate that the physical interpretation is not the same as for renormalized PDFs in QCD. There are also alternative ways to quantify uncertainties in the EWA [79,80,83]. For representative collider energies √ s = 4, 14, and 30 TeV, the beam-level cross sections under the EWA (σ EWA ), scale uncertainties [%], and polarization fractions [%] for W + λ A W − λ B → HH and γγγ are summarized in the top two panels of Table 2. For comparison, we show in the lower panel of Table 2 the corresponding cross sections (σ Full ) using the full MEs, i.e., without the EWA. The sizable differences between σ EVA and σ Full , as well as the large scale uncertainties of the EWA result, will now be discussed.

Dependence on hard-scattering scale
We start our presentation on EWA uncertainties with what we find to be the most telling: that the accuracy of EWA cross sections for VBF depends crucially on the size of (M 2 Figure 2. power corrections. To show this, we plot in Fig. 2(a) the invariant mass distributions at √ s = 4 TeV of the (HH)-system using the full 2 → 4 ME (solid) and the EVA 2 → 2 ME (dashed). We assume two scenarios: one where the SM Higgs vev is its usual value √ 2 Φ = v SM ≈ 246 GeV (dark, lower curves), and a hypothetical situation where the vev is reduced by a factor 10 (light, upper curves), i.e., where √ 2 Φ = v SM /10 ≈ 24.6 GeV. In the small-vev scenario, we keep M H and all EW gauge couplings to be their SM values in the Thomson limit. This implies M W ≈ 8.04 GeV and M Z ≈ 9.14 GeV.
Focusing first on the SM case, we clearly see that the EWA and the full ME computations are in agreement for M (HH) 1 TeV. Below this threshold, the EWA curve significantly overestimate the full ME. In the lowest bins, the differences between the curves reach approximately factors of 3 − 5. This excess in the EVA prediction accounts for the differences in cross sections reported in Table 2. Differences between the full ME and f W ± 0 PDFs consist of corrections associated with expanding in powers of (M 2 W /M 2 (HH)) and (p ν2 T /M 2 (HH)) ∼ (M 2 W /M 2 (HH)). Importantly, we can rule out a meaningful dependence on µ f since W + λ A W − λ B → HH is driven almost exclusively by W + 0 W − 0 scattering. To check that these power corrections are driving the disparity between the full and approximated MEs, we turn to the reduced-vev case. Remarkably, if we reduce M W by a factor of 10, the disagreement between the EWA and the full ME disappears to within MC statistical uncertainties.
The same scenarios are presented for γγγ production in Fig. 2 (b). There, we plot the invariant mass distribution of the (γγγ)-system for the full (solid) and EWA MEs. As W + λ A W − λ B → γγγ is driven by W + T W − T scattering, there is an ambiguity associated with our choice for µ f in the W ± T PDFs. Therefore, we consider the envelop spanned by setting µ f = √ŝ (dash) and µ f = √ŝ /4 (dot). In the SM case (light curves), the scale uncertainty envelope spans a huge gap that sandwiches the full ME for M (γγγ) 750 GeV. This large scale variation can be understood by considering the logarithms in the W ± T PDFs themselves. For a fixed M (γγγ), the ratio of the two EVA distributions is given by (4.7) For the representative triphoton invariant masses M (γγγ) ∈ {0.5, 1, 2, 2.5} TeV, we obtain roughly the respective ratios {17, 4.9, 3.1, 2.8}, in agreement with the distribution. For M (γγγ) 750 GeV, the full ME curve sits just above the EWA envelope. This is in contrast to W + 0 W − 0 → HH, where the full ME distribution sits below the EWA rate. We attribute this to W + T W − T scattering having a weaker dependence on power corrections than W + 0 W − 0 scattering. Differences between the full ME and f W ± T PDFs are associated with expanding in powers of (p ν2 the distribution of the full ME approaches the EWA curve for µ f = √ŝ /4, suggesting a preferred choice for setting µ f . In the reduced-vev scenario (dark curves), we observe several noteworthy features. First is an improved agreement between the full and EWA distributions for M (γγγ) 250 GeV. For even lower invariant masses, the full ME is again higher than the EWA band. Second, we find that the full ME converges to the EWA curve for µ f = √ŝ /4 when M (γγγ) 750 TeV. Third is the appearance of a smaller scale uncertainty envelope, in accordance with Eq. (4.7). Numerically, this follows from the fact that for small variations of the argument x, the quantity log(x) varies less when x is large than when x is near unity. Physically, this means that in the reduced-vev case, typical M (γγγ) are further away from the W 's mass threshold, and therefore is less sensitive to O(log(M 2 (γγγ)/M 2 W )) variations. Despite being smaller in this scenario, we stress that the scale uncertainty band remains sizable. For instance: using Eq. (4.7) and our benchmark values for M (γγγ), we obtain the ratios {2.3, 2.0, 1.8, 1.7}. This indicates that for realistic EW boson masses, one must go to asymptotically large M (γγγ) in order to obtain O(10% − 20%) uncertainties. From an alternative perspective, the large µ f dependence is indicative of the need to extend the formula of Eq. (2.2) by an EW Sudakov form factor and/or RG evolution for weak boson PDFs, as studied in Refs. [4,26,56,89,90,93,[116][117][118].
From these distributions, we can conclude that the EWA is acutely sensitive to power corrections of the form (M 2 . This is particularly true when scattering longitudinally polarized weak bosons. Distributions also suggest a weaker dependence on power Table 3. Same as Table 2 but requiring M (HH) > 1 TeV and M (γγγ) > 1 TeV.
corrections when scattering transversely polarized weak bosons. We attribute this difference to the different expansions needed to derive longitudinal and transverse weak boson PDFs: W T PDFs require a single power expansion whereas W 0 PDFs require a double expansion. Altogether, this points to evidence of the EWA's success for both transverse and To further demonstrate this at the level of cross sections, we show in Table 3 the same quantities as in Table 2 but require also that M (HH) > 1 TeV and M (γγγ) > 1 TeV. The improved agreement between the full and EWA computations is due to the cuts on M (W W ).

Dependence on collider energy
In light of the above, we consider now the impact collider energy on the EWA's accuracy. Increasing √ s has two prominent effects on V λ A V λ B scattering: (i) For fixed momentum fractions ξ 1 and ξ 2 , more energetic µ + µ − collisions lead to more energetic (V V )-systems, with M 2 V V = ξ 1 ξ 2 s. The corresponding enhancement of collinear logarithms indicates an enlargement of collinear regions of phase space. (ii) For a fixed hard-scattering scale M V V , increasing the collider energy leads to probing smaller ξ 1 and ξ 2 . The corresponding enhancement of soft logarithms similarly indicates an enlargement of soft regions of phase space. (Soft logarithms appear after integratingf V λ (ξ i ) ∼ 1/ξ i over ξ i ; see, e.g., Ref. [42].) To explore these effects, we present in Fig. 3 the invariant mass distributions of (a) the (HH)-system and (b) the (γγγ)-system at √ s = 4 TeV (light), 14 TeV (darker), and 30 TeV (darkest), assuming the full MEs for the processes in Eq. (4.1) (solid), and the EWA MEs for the processes in Eq. (4.3) (dash). For the W + T W − T → γγγ process, we show the scale variation envelop obtained by setting µ f = √ŝ (dash) and µ f = √ŝ /4 (dot). Focusing on the M (HH) distribution in Fig. 3(a), several observations can be made. We start with the anticipated jump in cross section for increasing √ s. For both the EWA and the full MEs, we find that increasing the collider energy by a factor of 3.5 causes all total cross sections to increase by about a factor of 3 (see Table 2). Increases are much more dramatic at the differential level for M (HH) 1.5 TeV due to the significant opening of phase space. In this regime, we also find good agreement with the normalization and shape between the EWA and full MEs. At lower invariant masses, particularly for M (HH) 500 GeV, we find that the EWA overestimates the full ME in the same manner as observed in the previous section. In this regime, the EWA distributions increase more quickly with rising √ s than the full ME distributions: in the lowest M (HH) bins, the EWA ME overestimates the full ME by about a factor of 3 − 5 at √ s = 4 TeV and by about 3.4 − 4.7 at 30 TeV. As longitudinal weak boson PDFs do not contain collinear logarithms, the enhancements in Fig. 3(a) are driven exclusively by soft logarithms. This implies that the EWA favors the production of relatively softer W ± 0 , and hence lower M (HH), a phenomenon that is sometimes [4] described as "ultra collinear enhancements." Consequentially, increasing the collider energy reinforces the sensitivity to (M 2 W /M 2 (HH)) power corrections, which must be negative. Despite this, the distributions show that regardless of √ s the EWA converges to the full ME computation for M (HH) 1 TeV.
Turning to the M (γγγ) distribution in Fig. 3(b), we observe several of the same characteristics. Foremost we find that the full ME distribution consistently sits within the EVA scale uncertainty band for M (γγγ) 750 GeV − 1 TeV, for all √ s = 4 − 30 TeV. Though, for increasing √ s we find that the full ME expectation migrates away from the µ f = √ s/4 boundary and towards the envelope's center. For a fixed M (γγγ), we find that the thicknesses of the µ f uncertainty bands remain about the same for increasing √ s, with changes just outside MC statistical uncertainties. This is consistent with the ratio expression of Eq. (4.7), which does not obviously suggest an additional dependence on collider energy once M (γγγ) is fixed. An important difference with respect to the M (HH) case is the preference for larger values of M (γγγ) with increasing √ s. (For W + 0 W − 0 → HH, smaller invariant masses are preferred at increasing collider energy.) As W + T W − T → γγγ is devoid of possible ME-level enhancements from longitudinal polarizations, we attribute these behaviors to the collinear logarithm contained in the f W ± T PDF, which favors producing larger invariant masses. Notably, the collinear logarithms reinforce the accuracy of the EVA by favoring phase space regions where (M 2 W /M 2 (γγγ)) is small.

Dependence on evolution variable and phase space cuts
As shown in Sec. 2, EW boson PDFs can be constructed using either the virtuality q of V λ or transverse momentum p T of l as the evolution variable in collinear µ → V λ l splitting. Given that the expressions for the f W ± T PDFs differ, we investigate whether these schemes give appreciably different results. To explore this, we focus on the process PDFs are the same under both schemes.
In Fig. 4(a), we show the invariant mass distribution of the (γγγ)-system in W + T W − T → γγγ at √ s = 4 TeV using the full ME for the 2 → 5 scattering process (solid) and the analogous process under the EWA. For the EWA curves, we consider the scale variation envelop spanned by setting µ f = √ŝ (dash) and µ f = √ŝ /4 (dot) for evolution by q 2 (thick curves) and evolution by p 2 T (thin curves). We also impose the fiducial cuts in Eq. (4.4). As found in previous subsections, the distribution for the full ME sits in the EWA envelope for M (γγγ) 750 GeV, when assuming evolution by q 2 . In comparison, the evolution-by-p 2 T envelope is systematically shifted upward to a larger set of rates. Notably, this leads to the µ f = √ŝ /4 curve for evolution by p 2 T to overestimate the full ME when M (γγγ) 2 TeV. This shift is entirely due to the different arguments in the PDFs' collinear logarithms.
More specifically, for evolution by . This difference implies that for a fixed factorization scale, the logarithm in the evolution-by-p 2 T PDF is relatively enhanced by a factor of 1/(1 − ξ). This favors larger momentum fractions and therefore harder invariant mass distributions since M 2 (γγγ) = ξ 1 ξ 2 s. This enhancement and subsequent disagreement do not imply that evolution-by-p 2 T scheme itself is incorrect. It only indicates that setting µ f = √ŝ or µ f = √ŝ /4 are poor choices of factorization scale for this process and collider energy. Obviously, setting µ f = (1 − ξ)ŝ or µ f = (1 − ξ)ŝ/4, which is a common practice SCET [98,99], would recover the results for evolving by q 2 .
In Fig. 4(b), we show the invariant mass distribution but for √ s = 14 TeV. Remarkably, we find that two the envelopes converge, and differences between the two evolution schemes are within MC statistical uncertainties. Importantly, the full ME distribution remains inside both envelopes for M (γγγ) 750 GeV. To understand this improved agreement, recall that when the hard scattering scale M (γγγ) is fixed, one probes smaller momentum fractions for increasing collider energies. This implies that as √ s increases, the 1/(1 − ξ) enhancements cease and systematically approach unity, i.e., 1/(1 − ξ) → 1. In Fig. 4(c) and (d), we investigate the impact of tighter phase space cuts on the final-state system. In particular, we consider the cases when photons are (c) more central with |η γ | < 1.5, and (d) harder with p γ T > 150 GeV. Aside from an obvious reduction in cross section, a few qualitative observations are worth reporting. In figure (c), we observe a shift in the distribution to smaller M (γγγ), whereas in figure (d) the shift is to larger M (γγγ). In both cases, the change is kinematical. Consider the invariant mass of the entire three-photon system in terms of two-photon systems, i.e., Here, η i and φ i are the pseudorapidity and azimuth of γ i . Requiring a smaller |η γ | window leads to smaller diphoton masses, and hence smaller triphoton masses. Likewise, requiring larger p γ T leads to larger diphoton masses, and therefore larger triphoton masses. Despite these different shifts, we do not find clear improvement or worsening of the agreement between either evolution scheme and the full ME. In both (c) and (d), the full ME distribution sits inside both envelopes for M (γγγ) 750 GeV. However, as in the baseline case, the p 2 T -scheme again overestimates the full ME for M (γγγ) 2 TeV due to too large log(1 − ξ) enhancements. This suggests that the specific dynamics of the hard scattering process may not have an appreciable impact the validity of the EWA, as one would hope.

Matrix Element Matching with Collinear W T PDFs
As a final check of our implementation of EW boson PDFs into mg5amc and as a proofof-concept demonstration of the potential capabilities, we briefly explore matrix element matching (MEM) with transverse weak boson PDFs. The idea behind MEM is that one can divide computations for complicated, many-leg final states that are susceptible to numerical instabilities, e.g., µ + µ − → ν µ ν µ γγγ, into two easier, more stable parts: (i) a ME with a smaller final-state multiplicity that represents a particular region of phase space of the original process, e.g., the process µ + W − → ν µ γγγ with af W − /µ − PDF, which describes the collinear µ − → W − ν µ splitting; and (ii) the ME for the original process but where the phase space for (i) is excluded, e.g., µ + µ − → ν µ ν µ γγγ process with only wide-angle µ − → W − ν µ splittings. In principle, summing the two components should recover the full phase space for the original ME, up to power corrections that are formally small. The aim of this procedure is to efficiently describe regions of phase space that are otherwise difficult to model simultaneously due to instabilities associated with soft and/or collinear radiation. If MEM is successfully implemented, then the sum of (i) and (ii) should not only reproduce the original cross section, up to uncertainties, but also display an insensitivity to the (artificial) cutoff scale that divides the original process into regions (i) or (ii).
This subsection serves as a proof-of-concept check and exploration of the power-law-like corrections in the factorization formula of Eq. (2.2). The discussion here also touches upon whether there is a natural or preferred choice for µ f , which might resolve the collinear/wideangle ambiguity described in Sec. 4.2. (It is not obvious that arguments for setting µ f in hadron collisions, such as those given in Ref. [114,115], are applicable here.) As discussed below, there are technical nuances at play that merit comprehensive exploration. However, this is beyond the scope of our work. Future studies that expand on this section are therefore encouraged.
To sketch MEM conceptually for the case of matching collinear W T (or Z T ) PDFs with full ME, we focus on the scattering process e + µ − → ν e ν µ γγγ. One can schematically divide the cross section of the process (σ Tot. ) into three disjoint pieces that describe different modes of µ − → W − ν µ splitting: (i) a collinear piece (σ C ), (ii) a quasi-collinear piece (σ QC ), and (iii) a hard piece (σ H ). Taking p νµ T as the evolution variable, one can write: (4.10) Here, µ f M W is a factorization scale separating collinear and wide-angle µ − → W − ν µ splittings. The second cutoff Λ, which satisfies µ f Λ p max T , is some arbitrary scale such that collinear factorization remains a good approximation. (The introduction of Λ is simply for bookkeeping: it ensures σ QC can be written in a convenient manner.) Finally, p max T is upper bound on p ν T allowed by momentum conservation and phase space cuts. After integration over p ν T , both σ C and σ QC will scale like a collinear logarithm and power corrections, which we neglect (retain) in the (quasi-)collinear expression: Combining the quasi-collinear and hard terms, one obtains the cross section (σ W ≡ σ QC + σ H ) for wide-angle µ − → W − ν µ splitting that is independent of Λ (since all Λ-dependent terms are kept). Moreover, were one to combine the leading logarithmic term of σ C with σ W , then the logarithmic dependence on µ f would vanish identically. One would also recover the total cross section, up to the neglected power corrections. Explicitly, one finds This indicates that µ f can be interpreted in MEM also as the "matching scale" that matching collinear and wide-angle regions of phase space in an inclusive calculation. To demonstrate that MEM is possible with transverse weak boson PDFs, we focus on W + W − T → γγγ in e + µ − collisions and define the following disjoint regions of phase space: When mediated by the EWA, we can identify e + W − T scattering in Eq. (4.14) as the collinear component of the inclusive e + µ − → ν e ν µ γγγ process. Analogously, we can identify Eq. (4.15) as the wide-angle component of the inclusive e + µ − → ν e ν µ γγγ process, when the appropriate phase space cut is applied to p νµ T or the µ − → ν µ momentum transfer q µνµ .
To regulate poles associated with final-state photons, to avoid instabilities associated with collinear e + → W + ν e splittings, and to minimize the power corrections described in Sec. 4.3, we impose the following phase space restrictions on both Eqs. (4.14) and (4.15): We treat the initial-state W − T using the appropriate PDF and consider when the PDF is defined in terms of (i) q 2 as given in Eq. (2.7), and (ii) p 2 T as given in Eq. (2.8). Assuming thatf W − T has been evaluated at a factorization scale µ f , then for the case of evolution by q 2 , we remove collinear and shallow-angle splittings in Eq. (4.15) by requiring that the norm of the µ − → ν µ momentum transfer is larger than µ f . Symbolically, this is given by This cut is implemented into mg5amc through the dummy_cuts function (file dummy_fct.f).
For the case of evolution by p 2 T , we require that ν µ satisfies the following restriction In Fig. 5, we show as a function of the factorization scale at √ s = 10 TeV: (i) the total fiducial cross section for e + µ − → ν e ν µ γγγ without restrictions on the µ − → W − splitting (flat black curve labeled "Total"); (ii) the same but with restrictions on the µ − → W − splitting (light purple curve labeled "Full ME + q cut" or "Full ME + p T cut"); (iii) the fiducial cross section with its scale uncertainty band for e + W − T → ν e γγγ (light green band labeled "EWA"); and (iv) the sum of the restricted cross section and EWA with scale uncertainty (dark blue or red bands labeled "Sum"). More specifically, we show the dependence when µ f = |q 2 µνµ | in Fig. 5(a) and when µ f = p To establish a baseline, we start with Fig. 5(a) for evolution by q µνµ at M (γγγ) > 1 TeV. As one can anticipate, the fiducial cross sections with and without restrictions on µ − → W − splittings converge to the same rate, about σ ∼ 0.3 fb, when effectively no cut is placed on the momentum transfer variable |q 2 µνµ |, i.e., when µ f → 0 TeV. For increasing µ f , we observe a logarithmic-like dependence for both the full ME with a |q 2 µνµ | cut as well as the EWA result. In the former (latter), the rate decreases (increases) with increasing µ f . We observe that the restricted rate and the EWA rate are comparable over the approximate range µ f ∈ (200 GeV, 500 GeV). For µ f 1.5 TeV − 2 TeV, we find that the full ME with a |q 2 µνµ | cut becomes negligible, whereas the EWA result remains comparable to the total cross section. Over the range of µ f investigated, we find that the EWA uncertainty band spans from about ±75% at µ f ∼ 200 GeV to about ±20% at µ f 2.5 TeV; we caution, however, that the large change to the uncertainty band reflects the extreme range of µ f we   investigate ) Importantly, when adding the restricted and EWA cross sections, we find that the sum reproduces the total fiducial cross section to within uncertainties and shows a strong insensitivity to the matching scale for µ f < 2 TeV. This indicates that MEM was achieved. While not shown, we report that the central value for the summed result consistently overestimates the total result by about 3% − 7% for µ f = 200 GeV − 1 TeV, and up to 20%−30% for µ f = 2 TeV−3 TeV. For the case of M (γγγ) > 3 TeV, we observe much of the same qualitative and quantitative behavior. The two notable differences are: (i) the obvious reduction in cross sections due to a more restrictive phase space cut, and (ii) a more stable summed result that features a central value (curve not shown) sitting nearly uniformly at about 10% above the total result for µ f up to µ f = 1 TeV.
Much can be learned from this exercise. First is that the power corrections (∆σ PC ) that distinguish the total fiducial cross section without t-channel cuts (σ Total ) from the summed result (σ Sum ), and which scale as ∆σ PC = σ Total − σ Sum ∼ O(µ 2 f /M 2 W W ), are negative. (This is clear from the fact that the summed rate exceeds the total rate.) Second is that the W T PDFs, and by extension the Z T PDFs, work best when µ f is set to lower values, e.g., µ f 1 TeV when √ s = 10 TeV. Naïvely, requiring relatively small µ f may appear at odds with standard practices in pQCD, where µ f are typically very high. However, in pQCD, one nearly always uses QCD PDFs that have been DGLAP-evolved; this has the effect of reducing the size of PDFs due to the RG running of α s (µ r ), thereby compensating for large collinear logarithms. The EW PDFs throughout this work are not evolved with EW-DGLAP equations and necessitates smaller µ f . The third observation is that for sufficiently large values of µ f , the EWA rate begins to overestimate the total cross section. With little doubt, this can be attributed to a breakdown of the collinear approximation, which requires collinear initial-state splittings, i.e., |q 2 µνµ | ∼ p ν2 T M 2 W W . Fourth is that the relative independence of the summed result on a matching scale indicates that the logarithmic dependence on µ f in each component effectively cancel, in accordance with expectations. This serves as a highly non-trivial check of MEM with W T /Z T PDFs but also demonstrates the potential to support it in MC event generators. Finally, we note that the very similar and consistent size of the scale uncertainty bands across all channels can be attributed to the fact that we are working with fixed µ f . For example: in (a), the scale variations are given by the ratio σ(µ f = ζµ 0 )/σ(µ f = µ 0 ) = log(ζ 2 µ 2 0 /M 2 W )/ log(µ 2 0 /M 2 W ) with ζ ∈ {0.5, 1, 2}. Focusing now on Fig. 5(b), we show the same quantities as in Fig. 5(a) but for evolution by p νµ T . Qualitatively, we see many similarities to the previous case, including that the summed result reproduces the total fiducial cross section to within uncertainties for µ 1.5 TeV. For larger µ f , the difference between the summed and total results exceeds the uncertainty band of the summed result. We attribute this breakdown of MEM to a breakdown of the collinear approximation. The breakdown is more explicit in this case since one is varying p νµ T and is influenced by the large-ξ enhancement discussed in Sec. 4.5.
To further explore the dependence on the kinematics of the hard W + W − T → γγγ scattering process, we show in Figs. 5(c) and 5(d) the same quantities as in Figs. 5(a), but consider the more restrictive phase space cuts (c) |η| < 1.5 and (d) p γ T > 150 GeV. We again observe many of the same qualitative features, indicating some degree of independence from the hard scattering process. (Fewer changes would suggest more universal-like behavior.) One notable difference is that the tighter η γ restriction helps extend the agreement between the summed and total results out to about µ f ∼ 3 TeV. Slightly better agreement is also observed for the tighter p γ T requirement, but only until µ f ∼ 2.25 TeV. Finally in Fig. 6(a) and (b), we show the same quantities as in Fig. 5(a) and (b) but for √ s = 30 TeV. As the qualitative and quantitative findings are highly comparable, little needs to be said. One noteworthy difference, however, is that higher collider energy further alleviates differences between the total and summed results, and extends the agreement to µ f 3 TeV. For all cases, the "summed" rate remains bigger than the "total" rate.

Polarized Vector Boson Scattering at Muon Colliders
A scenario in which the EVA promises to be highly relevant is that of a multi-TeV µ + µ − collider. It is worth reiterating that when considering a generic multi-particle state F, the VBF mechanism becomes an increasingly important, if not dominant, production vehicle in lepton collisions as the energy increases [42]. Given this, it is natural to consider muon colliders as effective weak boson colliders and take full advantage of the EVA in order to systematically organize and simplify the precision of scattering computations.
In this section, we explore the production of SM states generically parameterized by where F contains up to n F = 4 unpolarized states from the collection {H, t, t, W + , W − , Z, γ}.
All helicity polarizations are defined in the hard-scattering frame, i.e., the rest frame of F. We consider collider configurations over the range √ s = 2 − 30 TeV, and require final-state particles to obey the following kinematic and fiducial cuts: The invariant mass cut of 1 TeV on the system F is needed to ensure that power corrections of the form (M 2 V /M 2 (F)) are negligible, in accordance with findings of Sec. 4. We require moderate rapidities y to avoid t-and u-channel singularities and instabilities associated with final-state particles, as advocated by Ref. [26]; in the massless limit, the pseudorapidity value of y → η = 3 corresponds to a polar angle of θ ≈ 5.7 • . For all calculations involving V T V T or V T V 0 scattering, we set the central collinear factorization scale according to Eq. (4.5) and display three-point scale uncertainties. Scale uncertainties are unavailable for V 0 V 0 scattering as thef V 0 PDFs do not depend on µ f . Our survey is organized in the following manner: We start in Sec. 5.1 with associated and multi-Higgs production. In Sec. 5.2, top quark and associated top quark production are discussed, followed by diboson and triboson production in Secs. 5.3 and 5.4, respectively.

Higgs production
We first consider Higgs production in µ + µ − collisions and focus on the channels At the unpolarized level in Fig. 7(a), we observe a strong hierarchy between associated and multi-Higgs production, with V H production rates being more than an order of magnitude larger than HH production. More specifically, the ZH and W + H rates span roughly σ ∼ 5 − 25 fb for √ s ∼ 5 − 30 TeV, while HH production reaches about σ ∼ 0.5 − 1 fb over the same range. By C-symmetry, the W − H production rate is the same as W + H, and therefore is not shown. The production of γH is universally smaller than ZH and W H by about a factor of 2 − 3. Since ZH and γH are both mediated by W + W − scattering, the difference can be attributed to the difference in W W Z and W W γ gauge couplings, where σ ZH /σ γH ∼ (g 2 ZW W /g 2 γW W W ) = (g cos θ W /g sin θ W ) 2 ∼ 3. As the HH channel is also driven by W + W − fusion (recall that the Zµµ coupling is smaller than the W µν µ coupling), it is tempting to also attribute the relative size of the V H and HH production rates to the coupling ratio (g/λ) 2 ∼ 27. However, as shown in Fig. 7(b) and discussed in the next paragraph, this is not actually the case. All V H and HH channels are about 2-to-3 orders of magnitude larger than HHH production, which is yet another 2-to-3 orders of magnitude larger than HHHH production. The (relatively) tiny triple and quadruple Higgs cross sections follow from the compound effect of a small Higgs self-coupling and phase space suppression. For √ s 5 TeV, all rates are sensitive to small changes in collider energy due to threshold effects; above this scale, the energy dependence becomes milder. As a function of polarization, one sees from Figs. 7(b-d) several notable characteristics. For instance: The V H channels are driven almost exclusively by V 0 V T scattering, with a subleading component of V 0 V 0 scattering. HH is dominated by V 0 V 0 scattering, which accounts for the smaller dependence on factorization scales when summing over all polarizations, but also contains a sub-leading V T V T contribution. Notably, ZH, W + H, and HH all have comparable V 0 V 0 scattering rates, which is in line with the Goldstone Equivalence Theorem. This indicates that the σ V H -σ HH hierarchy observed in Figs. 7(a) is actually due to the compound effect of logarithmic enhancements inf V T PDFs and helicity configurations allowed by angular momentum conservation, e.g., V T V 0 → HH is helicity suppressed. We report that HHH production has significant and comparable contributions from V 0 V 0 and V T V 0 scattering, but only a marginal contribution from V T V T . Interestingly, HHHH Table 4. Fiducial cross sections with scale uncertainties for associated Higgs and multi-Higgs production in µ + µ − collisions through V λ A V λ B scattering, where V λ ∈ {γ λ , Z λ , W ± λ }, under the EVA for polarization-summed and polarized initial states, at √ s = 3, 14, and 30 TeV. All helicity polarizations are defined in the hard-scattering frame. Phase space cuts are summarized in Eq. (5.2). Also shown is the mg5amc syntax for modeling the process assuming the multi-particle definitions "vxp=z,w+,a" and "vxm=z,w-,a." The configuration (0, T ) implies a sum over both (0, T ) and (T, 0), and uses the syntax generate vxp{0} vxm{T} > ...; add process vxp{T} vxm{0} > ... receives comparable contributions from all polarization configurations. As one can expect, production from V T V T scattering exhibits larger scale uncertainties than in V T V 0 scattering. For each unpolarized and polarized scattering configuration, we document in Table 4 the relevant mg5amc process syntax that enables our computation and the fiducial cross section [fb] with its scale uncertainty [%], for representative collider energies.

Top and associated top production
Next, we address tt pair and associated tt production, focusing on the channels V λ A V λ B → tt, ttZ, ttW + , ttγ, ttH, and tttt. For the tttt process, we consider both the production though QCD and EW couplings (tttt QCD+EW ) as well as purely through EW couplings (tttt EW ). In Fig. 8 we present Figure 8. Same as Fig. 7 but for tt pair production and tt + X associated production.
the unpolarized and helicity-polarized cross sections as a function of collider energy in the same manner as in Fig. 7. An immediate observation is that for all processes each of the polarization combinations of initial-state EW bosons V λ A V λ B give a comparable contribution to the unpolarized process. This is in contrast to associated and multi-Higgs production in Sec. 5.1, where typically one particular V λ A V λ B configuration drives the total process. For all polarized and unpolarized cases, tt production exhibits the largest cross sections, whereas pure EW production of tttt exhibits the lowest rates. The difference between the channels is about three orders of magnitude. Mixed QCD+EW production of tttt sits just above the pure EW rate; notably, the difference between the two tttt processes Table 5. Same as Table 4 but for tt pair production and tt associated production.
is larger than their uncertainty bands. All ttV and ttH cross sections are sandwiched between the two processes and exhibit a rate hierarchy that is in line with naïve EW coupling enhancement/suppression. While the hierarchy is independent of collider energy, some reordering can be observed when individual helicity configurations are considered.
Focusing first on V 0 V 0 scattering in Fig. 8(b), we find that the ttV and ttH channels all have highly comparable rates. This is essentially due to all these channels being driven by either W + 0 W − 0 or W + 0 Z 0 scattering. (There is no γ 0 PDF.) Hence, appreciable differences in rate are due to differences in coupling constants. Moreover, unlike other polarization configurations, there is no log(µ 2 f /m 2 µ ) enhancement since this is associated with γ T PDFs. More explicitly, for V T V 0 scattering in Fig. 8(c) and V T V T scattering in Fig. 8(d), we find appreciably larger cross sections than for pure V 0 V 0 scattering. Again, we attribute this to the opening of γ T V λ scattering and logarithmic enhancements in transverse PDFs. Interestingly, we find that the hierarchy of ttγ/Z/W + /H depends on the precise polarization configuration of initial-state EW bosons. For example: ttZ has a larger rate than ttW + in V T V 0 scattering, but the opposite is true in V T V T scattering. Both configurations lead to larger rates than for ttγ, which is not the case for V 0 V 0 scattering.  Figure 9. Same as Fig. 7 but for V V production.
We summarize these results for representative √ s in Table 5.

Diboson production
We now turn to diboson production in EVA. In Fig. 9, we plot again the (a) polarizationsummed and (b-d) polarized production cross section with their respective factorization scale uncertainties, as a function of collider energy for the following six processes V λ A V λ B → W + W − , ZZ, ZW + , γZ, γW + , and γγ. (5.5) As in previous cases, there are several global features that one can infer. Foremost is that unlike Higgs (Fig. 7) and top quark (Fig. 8) processes, we find a clear hierarchy among initial-state EW boson polarizations. More specifically, we find that V T V T scattering is categorically the dominant production vehicle of EW boson pairs. The role of V T V 0 scattering is about one-to-two orders of magnitude smaller than V T V T scattering, and rate of V 0 V 0 scattering is yet another decade smaller. As of V 0 V 0 scattering is negligible, unpolarized diboson cross sections in the EVA exhibit a relatively larger scale uncertainty than Higgs and top quark production, which feature a larger dependence on V 0 V 0 scattering. Another consequence of the strong polarization dependence is that the hierarchy of fiducial cross sections shown in the unpolarized case largely mirrors the hierarchy in V T V T scattering. That said, we find that this same hierarchy is mostly preserved with other helicity configurations. This suggests a larger dependence on available partonic channels and gauge couplings than logarithmic enhancements from soft and collinear regions of phase space. For example: the W + W − production cross section dominates for all helicity configurations but also can be produced via the most number of partonic configurations, i.e., The reverse can be said for γγ and γW + production. These processes exhibit the lowest diboson cross sections but can only proceed through one or two partonic channels, namely W + λ A W − λ B and W + λ A γ λ B . We also note that helicity suppression stemming from angular momentum conservation also plays a role in this hierarchy. For example: while W + 0 Z 0 → W + 0 γ T exhibits a larger longitudinal polarization enhancement than W + 0 Z 0 → W + T γ T , the former (later) is disfavored (favored) since it must proceed through a high-wave (low-wave) angular momentum configuration.
We summarize these results for representative √ s in Table 6.

Triboson production
In the final part of our survey, we show in Fig. 10 the fiducial cross sections for triboson production (a) after summing over all initial-state helicity polarizations and (b-d) for individual (λ A , λ B ) configurations. For conciseness, we focus on the representative channels An immediate observation we can make is the qualitative similarities between triboson and diboson production in Fig. 9. In particular, we find that triboson production is driven by largely V T V T scattering. The V 0 V T rate is an order of magnitude or two smaller, and the V 0 V 0 rate is smaller by about one or two additional decades. The precise polarization composition has a slight dependence on the underlying process: The production of ZW + W − categorically exhibits the largest cross section, with In other words, ZW + W − production is very roughly Table 6. Same as Table 4 but for V V production.
Other similarities between triboson and diboson production include the cross section hierarchy that one observes in V T V T scattering largely appears also in V 0 V T scattering and V 0 V 0 scattering. One exception is γγW + and γγγ production in V 0 V 0 scattering, where the ordering inverts. We attribute this to the opening of orbital angular momentum configurations in 2 → 3 scattering, which can spoil the cancellations described in Sec. 5.3. While it is beyond our immediate scope, one can, in principle, investigate the hierarchy of triboson processes by considering the relative important of specific partonic and helicity channels, e.g., γ − γ + scattering or Z + Z 0 scattering. One can also further decompose the final-state triboson system into its helicity components, thereby allowing one to investigate the relative importance of high-and low-wave angular momentum configurations. Finally, due to the similarities between diboson and triboson production, we conjecture that comparable behavior will be observed for the production for four or more EW vector bosons.
We summarize these results for representative √ s in Table 7. Figure 10. Same as Fig. 7 but for V V V production.

Discussion, Outlook, and Conclusions
As a weakly coupled, non-Abelian gauge theory, the weak sector of the SM naturally exhibits similarities to the electromagnetic and strong sectors, e.g., coupling universality. However, as the weak sector is also spontaneously broken, scattering and decay rates involving weak bosons feature a power-law dependence on the scale ratio (M 2 V /Q 2 ) k or (v 2 /Q 2 ) k , with k > 0. In the absence of such contributions, i.e., when momentum-transfer scales are much larger than the EW scale, meaning that O(M 2 V /Q 2 ) terms are negligible, scattering and decay rates involving weak gauge bosons resemble the analogous expressions for massless gauge Table 7. Same as Table 4 but for V V V production.
bosons in QED and pQCD. Importantly, this resemblance also holds for the factorization of weak-boson emission in the collinear and soft limits. Precisely how factorization (and resummation) operates in the weak sector, and in particular how it differs from QED and pQCD, is no longer an academic intrigue as multi-TeV µ + µ − colliders and 100 TeV pp colliders are being seriously discussed as eventual successors of the HL-LHC program [40,41]. At these colliders, typical parton collisions readily satisfy known criteria for collinear factorization of weak bosons. Even at the √ s = 13 TeV LHC, such hard-scattering scales have already been observed in measurements of VBF/S. Therefore, establishing a fuller picture of the colliders' physics potential requires a better understanding of how collinear and soft weak bosons behave in multi-TeV collisions. Motivated by the fact that multi-TeV µ + µ − colliders are effectively "high-luminosity weak boson colliders" [42], we have revisited the treatment of weak bosons as perturbative constituents of high-energy leptons, the so-called Effective W/Z Approximation [5,6]. To conduct this investigation, we have implemented PDFs at LO for helicity-polarized γ, W ± , and Z bosons from e ± and µ ± into the MC event generator MadGraph5_aMC@NLO. ‡ This ‡ These features were released publicly in version 3.3.0, and available from the mg5amc repository. allows for the fully differential simulation of scattering processes that are initiated by one or two initial-state EW bosons. Starting from a formula for scattering partons from a muon, which we state in Eq. (2.2), we systematically explored the limitations of the EWA in Sec. 4. Novelties of our comparative investigation are the focus on: (i) universal and quasi-universal power-law corrections, (ii) universal and quasi-universal logarithmic corrections, (iii) phase space dependence, (iv) helicity polarization, and (v) many-body processes with VBF. Past studies were mostly restricted to 2 → 1 and 2 → 2 scattering processes, which possess special kinematics, at a single choice of factorization scale. We stress that some important issues, such as Z T /γ T mixing and Sudakov resummation, were not investigated and are left to future work.
As documented in Sec. 4.3, a key conclusion is that the EWA is acutely sensitivity to power corrections of the form (p 2 These corrections originate in the derivation of weak boson PDFs and are related to the accuracy of collinear factorization and the Goldstone Equivalent Theorem. Our results suggest that using the EWA to describe VBF requires This is in addition to the typical assumptions needed to justify collinear factorization. When VBF systems carry larger invariant masses, we find that the difference between the full computation and the EWA are within factorization scale uncertainties. As documented in and around Eq. (4.7), these uncertainties can be very large and demonstrate a need for RG evolution in order to achieve precise results with weak boson PDFs. We caution that we restricted our attention to dynamic factorization scales that are proportional to M V V for consistency across the many processes that we surveyed. As in pQCD, more "optimal" choices probably exist but are also (probably) process dependent and should be investigated thoroughly. Our implementation of EW boson PDFs in mg5amc can facilitate such studies.
Importantly, we show that the size of non-universal power corrections and the size of factorization scale ambiguities in multi-TeV µ + µ − collisions are due to the largeness of the W and Z masses. At first, this may seem at odds with collinear factorization in pQCD, where above even a few GeV, PDFs can describe full matrix elements involving light quarks. However, using the operator product expansion, one can show that the phenomenon of "precocious scaling," i.e., the emergence of asymptotic freedom at moderate energies, is due to the smallness of parton masses [119][120][121]. Power corrections associated with quark masses m q are of the form (m 2 q /Q 2 ) k . For a scattering scale of Q ∼ 2 − 3 GeV, these reach at most O(m 2 q /Q 2 ) 10 −5 − 10 −2 . Likewise, heavy quark PDFs become adequate and reliable tools for O(m 2 q /Q 2 ) 0.01 [8][9][10][11]13]. Both are consistent with requiring that .01 for the EWA to adequately describe full matrix elements. For M V V O(1 TeV), we find contrasting behaviors between longitudinal and transverse weak boson PDFs: whereas longitudinal PDFs overestimate full scattering amplitudes, transverse PDFs underestimate them. As documented in Sec. 4.4, increasing the collider energy does not necessarily improve the accuracy of the EWA for M V V O(1 TeV). Whereas the presence of soft logarithms slightly improve the accuracy of transverse PDFs when √ s is increased, these same logarithms worsen the accuracy for longitudinal PDFs. It is clear that a matching, subtraction, or re-weighting scheme akin to those already available for pQCD and QED is needed to correct EWA matrix elements in this region. To facilitate such developments, we also report the availability in mg5amc of both q 2 -and p 2 T -dependent PDFs for transversely polarized EW bosons. As shown in Sec. 4.5, p 2 T -dependent PDFs consistently lead to larger cross sections in EWA, but converge to the q 2 -dependent results when √ s increases. To further strengthen the parallels with pQCD, we give a proof-ofprinciple demonstration in Sec. 4.6 of matrix-element matching of W T PDFs with the full matrix elements. Despite its formally large scale uncertainty band, the matched result shows significant independence on the matching scale. Broadly speaking, these capabilities provide a starting point for matching EWA matrix elements to EW parton showers and more sophisticated EW boson PDFs that involve RG evolution. Given these considerations, we cataloged in Sec. 5 a litany of processes of the form where F contains up to n F = 4 states from the collection {H, t, t, W + , W − , Z, γ}. In comparing polarized and polarization-summed cross sections, we find an intriguing interplay between helicity polarizations and hard scattering processes. For example: whereas multiboson and many-boson production is driven by the scattering of initial-states in the (λ A , λ B ) = (T, T ) configuration, top quark and associated top quark production features a large (0, T ) and (T, 0) component. In further contrast, multi-Higgs processes are dominated by (T, T ) and (0, 0) helicity configurations but receive only a marginal contribution from mixed configurations. It is worth noting that scale uncertainties at √ s = 3 − 30 TeV reach about δσ/σ ∼ 20% − 50% for many processes, which is beyond expectations based on coupling order. As we remained inclusive with respect to the helicities of final-state particles, the processes we surveyed and their differential behavior can all be further investigated.

Recommendations for using W/Z PDFs in high-energy lepton collisions
Finally, while much about PDFs for polarized and unpolarized weak bosons remains to be investigated, we believe this work helps clarify quantitatively when the W and Z bosons can be treated as partonic constituents of high-energy leptons. To further this prerogative, we provide a set of recommendations on using weak boson PDFs in many-TeV scattering calculations. These guidelines draw heavily from the findings in Sec.

A How to use EVA in MadGraph5_aMC@NLO
The simulation of polarized and unpolarized EW PDFs in high-energy charged lepton collisions, i.e., EVA, is possible using a series of commands inside the mg5amc interface. Instruction on how to run/setup MadGraph5_aMC@NLO for various configurations involving unpolarized matrix elements can be found in Ref. [106]; for the setup of polarized matrix elements, see Ref. [85]. In this appendix, we describe how to setup a EVA computation in mg5amc and particularly focus the new options and syntax introduced for this mode.
As a concrete example, we consider the hard scattering process W + W − → hh in µ + µ − collisions. A typical set of mg5amc commands to simulate a process like this is set group_subprocesses False generate w+ w-> h h output DIRECTORY_OUTPUT launch DIRECTORY_OUTPUT The command "set group_subprocesses False" is currently mandatory and deactivates some internal optimization mechanisms that are not (yet) compatible with EW boson PDFs as implemented into mg5amc. The second command corresponds to the hard process, and operates at the level of initial-state weak bosons. (In this sense, EW bosons are treated as partons of e ± and µ ± .) Note that EVA is only implemented here at LO in perturbation theory, without Z T /γ T mixing, and without EW-DGLAP evolution. Such corrections have a nontrivial impact on numerical results [4,26,56,[90][91][92][93]. The "output" command defines the directory where the code containing MEs and phase space integration routines, i.e., MadEvent [108], are physically written on disk. The "launch" command activates an interface to configure, compile, and execute this code. As mg5amc works by numerically evaluating helicity amplitudes, the command above syntax is equivalent [85]  In both cases, the particle species W ± λ is paired with the polarized PDF f W ± λ /µ ± (ξ, µ f ). When the user interface is initiated, i.e., just after the "launch" command, the user is prompted with the ability to edit multiple configuration files. To run in EVA mode, a user will need to edit the file DIRECTORY_OUTPUT/Cards/run_card.dat. This file contains all configuration details related to the beam, factorization scales, phase space restrictions (cuts), etc. The list of the important and new parameters are summarized in Table 8, and are described below. It is important to stress that in our implementation of EVA, initial-and final-state W and Z bosons retain their masses in all helicity amplitudes; nowhere do we set M W , M Z = 0 GeV. As a consequence, other mg5amc modules, such as MadSpin [122,123], can be employed in conjunction with EVA computations. This allows one to study, for example, the full process, W + T W − T → h(→ cc)h(→ bb). Investigating new physics remains possible through the interface [107] to Universal FeynRules Object (UFO) libraries [124]. We caution, however, that EW boson PDFs are hard-coded into files ElectroweakFlux.f and ElectroweakFluxDriver.f in the directory LO/Source/PDF/. This means that modifications to the − − γ/Z and − ν − W vertices introduced by a UFO will not propagate into the PDFs. We have designed and organized the calling of EW boson PDFs in mg5amc such that the W and Z boson masses as well as the EW couplings are automatically set to those values listed in the configuration file DIRECTORY_OUTPUT/Cards/param_card.dat. The values of the electron and muon masses are not read from the param_card.dat. Instead, the values listed in Eq. (3.3) are hard-coded into the file DIRECTORY_OUTPUT/Source/PDF/ElectroweakFlux.inc.
In order to initiate a computation with the EVA, the most important parameter that must be set in the file run_card.dat is the PDF set. Choosing EW boson PDFs for both beams can be done via the "pdlabel" parameter, which now accepts three additional modes: "eva" for the EW boson PDFs described in Sec. 2; "iww" for the so-called Improved Weizsäcker-Williams (IWW) γ PDF of Ref. [125]; and "mixed" for enabling different  PDF configurations for beams 1 and 2. A fourth option "none" deactivates PDFs for both beams. The γ PDF described in Sec. 2 is known historically as the Weizsäcker-Williams approximation [16,17,88], and is analogous to the gluon PDF in QCD at LO. Setting "pdlabel=iww" calls the γ PDF derived using the IWW approximation [125]; this PDF is sometimes mislabeled in the literature. Simply put, the IWW γ PDF augments the original Weizsäcker-Williams γ PDF by terms that correspond to operators in the operator product expansion with a twist larger than 2, i.e., are relatively suppressed by powers of (m 2 /Q 2 ). Just like for partons in hadron PDFs, the appropriate PDF and ME are paired automatically by the routines of pdg2pdf.f and pdg2pdf_lhapdf6.f in the directory LO/Source/PDF/.
Presently, it is possible in mg5amc to polarize electron and muon beams in leptonlepton and lepton-hadron collisions in the absence of lepton PDFs, i.e., "lppX=0" [106]. We have extended this capability and it is now also possible to polarize electron and muon beams when "lppX=±3,±4" and "pdlabel=eva" or "pdlabelX=eva." This is done via the run_card.dat parameters "polbeam1" and "polbeam2." A value of "polbeamX=0" (default) corresponds to an unpolarized beam, while "-100" and "+100" indicate, respectively, that § We report the correction of a long-standing labeling ambiguity of this PDF in the file LO/Source/PDF/ PhotonFlux.f, which contains the implementation of this PDF. We also reiterate that this PDF describes a photon from a proton in the elastic limit [126] and therefore does not include DGLAP evolution. For inelastic emissions of photons from protons, set "lppX=+1(-1)" and use a proton PDF set with QED evolution.
100% of beam X is polarized in the LH and RH helicity state. ¶ Note that the helicity polarization of an EW boson cannot be changed at run time via run_card.dat. It can only be specified when executing the "generate" command; see Ref. [85] for details.
To implement beam polarization with the EVA, we have augmented Eq. (2.9), which describes a polarized EW boson V λ from an unpolarized muon by the expressioñ Here, −1 ≤ β L ≤ 1 is a parameter describing the degree of LH polarization of the parent beam. For β L = 0, which corresponds to setting "polbeamX=0," one recovers Eq. (2.9). Likewise, setting β L = −1 (1), and implies that the muon beam itself is purely in the LH (RH) helicity state corresponds to setting "polbeamX=-100 (100)." The collinear factorization scale µ f that enters into EW boson PDFs in the EVA can be either dynamical, i.e., determined on an event-by-event basis, or fixed. This is chosen in run_card.dat via the Boolean parameter "fixed_fac_scale." Setting this parameter to "true (false)" activates a fixed (dynamical) µ f . Like "pdlabel," the dynamical/static scale scheme can be set separately for each beam using the two Boolean parameters "fixed_fac_scale1" and "fixed_fac_scale2." If a static µ f is selected, then its value is set in units of GeV by the parameters "dsqrt_q2fact1" and "dsqrt_q2fact2." For dynamical choices of µ f , a user can choose from predefined or user-defined scale schemes using the parameter "dynamical_scale_choice." For details on this, see Ref. [106]. It is possible to rescale µ f by the scale factor "scalefact". ‖ Importantly, automated scale variation of EW boson PDFs is possible using the "systematics" feature and setting "use_syst=True." Users are reminded that W/Z boson PDFs are not defined for values of µ f < M W/Z .

A.1 Example usage of EVA in MadGraph5_aMC@NLO
To reproduce the cross sections for the process V λ A V λ B → HH in Table 4, we use the following syntax to generate our (polarized) matrix elements and work environments: set group_subprocesses false set gauge Feynman define vxp = w+ z a define vxm = w-z a generate vxp vxm > h h output vxvx_hh generate vxp{T} vxm{T} > h h ¶ We note for clarity that setting "polbeamX=-70" indicates that 70% of beam X is polarized in the LH state while the remaining 30% is unpolarized. This implies that 85% (15%) of beam X consists of leptons in their LH (RH) helicity state.
‖ This new release of mg5amc is the first to allow both "scalefact" to be set different from unity and simultaneously allow "use_syst=True." However, one should note that the deprecated implementation of scale computation (SysCalc [128]) is not compatible with either the EVA or the new "scalefact" capability. The new and important commands for using the EVA in mg5amc are documented in App. A and summarized in Table 8.

B Effective W, Z, γ Approximation
In this appendix, we derive the (polarized) W/Z PDFs at LO and construct Eq. (2.2) as implemented in the MC event generator mg5amc. We include this for completeness and to make this work more self-contained for the non-expert reader. In making these mechanics explicit, we hope to minimize ambiguities in definitions of PDF that sometimes appear in the literature. Such differences are due to equally reasonable choices of, e.g., evolution variables, but can lead to sizable numerical differences at LO. We direct expert readers to derivations of the EVA in the axial gauge [4,7,79] and a fully covariant formalism [80] for considerations of more nuanced issues, such as Goldstone mixing and renormalization.
The construction proceeds in the following manner: after establishing notation and conventions in App. B.1, MEs for polarized collinear splitting functions are derived in App. B.2, with the corresponding phase space decomposition/factorization given in App. B.3. The unrenormalized, transverse-momentum-dependent distribution functions at LO are then derived in App. B.4, and finally the analogous collinear PDFs are given in App. B.5. While the construction follow closely that of Ref. [88] for parton splitting in QED, the identification of a factorization scale µ f as an ultraviolet regulator of a phase space integral is more closely aligned with Ref. [7]. This identification makes clearer that the procedures for implementing matrix-element matching (see Sec. 4.6) and (potential) RG evolution with renormalized PDFs are similar to those in pQCD.

B.1 Helicity Amplitude Notation and Conventions
To derive polarized, weak boson splitting functions and PDFs, we work at the helicity amplitude level, in the so-called HELAS basis [109]. We assume the spacetime metric signature g µν = diag(+, −, −, −), and work in the unitary gauge at dimension d = 4.

Spin-1/2 Particles
For a fermion with mass m = E 2 − | p| 2 , and a spin axis aligned with its momentum, its 4-momentum (p µ ) and transverse momentum (p T ) can be generically parameterized by p µ = (E, p x , p y , p z ) = (E, | p| sin θ cos φ, | p| sin θ sin φ, | p| cos θ), (B.1) The two-component helicity eigenstates with respect to the 3-momentum directionp are As such, the RH/LH chiral projection operators P R/L are respectively, Immediately one sees that the projection operators satisfy the identities 1 = P R + P L and γ 5 = P R − P L . (B.9)

Massive Spin-1 Particles
For a vector boson with mass M V = E 2 V − | q| 2 and a spin axis aligned with its momentum, its 4-momentum (q µ ) and transverse momentum (q T ) can be parameterized by q µ = (E V , q x , q y , q z ) = (E V , | q| sinθ cosφ, | q| sinθ sinφ, | q| cosθ), (B.10) In the Cartesian representation, its polarization vectors are given by , q x , q y , q z = γ(β, sinθ cosφ, sinθ sinφ, cosθ). (B.14) In the last line we used the Lorentz factor γ = E V /M V and the relationship | q| = βE V . In this representation, the above polarization vectors satisfy the orthogonality relationships: q µ ε µ (q, λ) = 0 and ε µ (q, λ)ε µ (q, λ ) = −δ λ,λ for λ = x, y, z. (B.15) Using the above, we can build the polarization vectors in the polar representation. For the RH (λ = +1), LH (λ = −1), and longitudinal (λ = 0) polarizations, these are given by For λ = ±1, the orthogonality relationships are modified such that In the EVA it is useful [5] to rewrite the λ = 0 vector in the following exact form: Explicit computation shows that the auxiliary vectorε obeys the following inner products: The purpose of this decomposition is two-fold: The first is to make manifest that the Goldstone contribution, i.e., the term that scales as ε µ (λ = 0) ∼ q µ /M V , vanishes when contracted with external parton currents via the Dirac equation. Ultimately, this is due to SU(2) L current conservation of massless leptons, meaning that massless leptons do not participate in helicity-inverting couplings of longitudinally weak boson. The second purpose is to make manifest that the non-vanishing part of ε µ (λ = 0) is formally a quasi-universal, beyond-twist-two term that scales as ε 0 ∼ M V /E. Massless vector bosons have identical transverse polarization vectors but do not possess a longitudinal polarization vector.

Setup and On-shell Decomposition
To build collinear splitting functions for weak bosons from high-energy leptons, we consider the deeply inelastic scattering (DIS) process shown in Fig, 11, given by The process above describes f (p A ) → f (p 1 ) + V (q) splitting with external fermions f, f , and an internal weak boson V with mass M V . Throughout this section, we remain agnostic to the composition of B. In App. B.5, we will extend our result and present a scattering formula involving a PDF for each lepton beam. As we are only working at LO in the EW theory for both the "hard" scattering (V B → X) and the splitting processes (f → f V ), we automatically work in the Single Boson Exchange Approximation. Furthermore, throughout this analysis, we neglect non-factorizable contributions. As we are working in the Unitary gauge (for clarity), the class of diagrams where this is an acceptable approximation is known to be limited in comparison to working in other gauges. For discussions on this, see Refs. [4,7,79,80,83,129,130]. For the process in Eq. (B.25), the associated ME is given by Here M ν (V * λ V B → X) describes the V * λ V B → X hard scattering process that occurs at the scale Q = (q + p B ) 2 ∼ E V M V ∼ |q 2 |. The propagator for V in the unitary gauge is denoted by ∆ V µν . Particle helicities are denoted by λ. For a V − f − f chiral interaction with LH (RH) coupling g f L (g f R ) and a universal strengthg, the f λ A → f λ 1 V * λ V fermion current can be written generically as (see Table 1 in Sec. 2.2 for explicit expressions) An important caveat to this expression is that it assumes we are working solely within the confines of the SM's SU(2) L gauge group (or some similarly broken gauge theory). It does not account for the fact that quarks and leptons in the SM both reside in the SU(3) c ⊗SU(2) L product group. While working in this larger structure has trivial implications for EW splitting functions themselves, it does impact d.o.f. counting and cancellations at the level of cross sections. More precisely, for the f λ A → f λ 1 V * λ V color-singlet splitting process, the V − f − f vertex is modified by a Kronecker δ-function that ensures color conservation. In equation (B.27), this causes the modification where the indices I and J denote the colors of f and f , respectively, and run over I, J = 1, . . . , N c = 3. For quarks, δ II = N c , whereas I, J = 1 and δ II = 1 for leptons.
The treatment of V 's propagator is subtle. In the language of Ref. [88], we are considering the kinematic regime where V is "almost real." For weak bosons, means we are assuming that the norm of V 's virtuality, its mass, and the difference of these two quantities, are all small compared to the hard scattering scale. We are also assuming that V propagates only a single helicity between the (f f ) and B systems since helicity inversion is suppressed by (M 2 V /Q 2 ) 1. Following Ref. [88], this entails making the replacement The idea of this (heuristic) replacement is that if V , which possesses an invariant lifetime τ V = /Γ V , is "almost" on-shell and if the hard scattering scale is sufficiently large such that Q Γ V , then V * is so long-lived compared to the hard process that it can be approximated as an asymptotic state with definite quantum numbers, e.g., with fixed helicity. Since V 's helicity λ V remains unchanged during propagation, the contribution from the unphysical / auxiliary polarization, which scales as ε µ (q, A) ∼ q µ q 2 − M 2 V , vanishes when contracted with J µ (f λ A → f λ 1 ). That is to say, since we assume f and f are massless, J(f λ A → f λ 1 ) · q = 0 by the Dirac equation. The net contribution of the propagator ∆ V µν (q) reduces to a coherent sum over physical polarizations λ V , i.e., a summation over λ V at the squared ME level. We note that this replacement can be made stronger when working in other gauges [4,7,79,80].
In terms of bras and kets, this decomposition is equivalent to writing the ME as which translates to a squared ME of Naïvely, this suggests a double summation over the helicities of V . However, when one is totally inclusive over the final-state X, then by unitarity 1 = X |X X|. This implies that the squared ME is only non-vanishing when the polarization of V in f → f V λ is the same as the polarization for V in V λ B → X. Other V λ |V λ combinations are orthogonal: The unitarity of final-state X acts as an effective Kronecker δ-function δ λλ . This reduce the double summation over helicities into a single summation over helicities, resulting in Returning to our derivation, after enjoining the outgoing (ε * µ ) and incoming (ε ν ) polarizations vectors with the f λ A → f λ 1 V λ V splitting amplitude and V λ V B → X scattering amplitude, respectively, one obtains the pair of MEs: For fixed polarizations of external particles f, f , and V , the ME for the process in Eq. (B.25) can be factorized into two spin-correlated sub-amplitudes: Immediately, one can write at the squared ME level: This indicates that if V goes on-shell and its mass is small compared to the scattering scale, then the probability density for f λ A B → f λ 1 X, when mediated by the space-like exchange of V λ A , can be approximated as the spin-correlated product of the f λ A → f λ 1 V λ V and V λ V B → X probability densities.

Helicity-dependent Collinear Splitting Amplitudes
To enact the on-shell condition such that Eqs. (B.37) and (B.38) are valid, we work in the kinematic configuration where V and f are emitted at shallow angles in f → f V splitting, i.e., the collinear limit. Formally, this involves working to leading order in the parameter λ ≡ p T /E A 1, where p T is the transverse momentum 2-vector of f (p 1 ) and E A is the energy of f (p A ). Defining z = E V /E A as the energy fraction carried by V (q) in the hard scattering frame, we can parameterize the momenta of the process in Eq. (B.25) by p µ A = E A (1, 0, 0, +1), p µ B = E A (1, 0, 0, −1), (B.39a) truncating wide-angle contributions at O(λ 2 ) results in the following approximations: Intuitively, this indicates that external particles are massless or approximately massless (to O(λ 2 )). The internal V , on the other hand, carries a virtuality −q 2 ∼ p T much smaller than the incoming energy E A . As V is also nearly on-shell, the scaling relationship M V ∼ −q 2 E A must consistently hold. Implicitly, we also work in the domain where the energy fraction z is far from its boundaries at z = 0 (the high-energy regime) and z = 1 (the threshold regime), which would otherwise necessitate resummation.
As f, f are massless, the only non-zero f λ A (p A ) → f λ 1 (p 1 ) currents (as defined in Eq. (B.27)) involving a vector emission are those that conserve helicity, namely (B.43) Employing the transverse polarization vector in Eq. (B.17) for V λ V =± , we obtain as the helicity amplitude (M) for the f L → f L V + splitting process In the second line of the above, we evaluated the ME exactly. In the third, we expanded the angular dependence to lowest order in the opening angles, and in the final line substituted θ V for p T . Repeating the same steps for the f L → f L V − splitting process, we obtain Note that the Dirac equation for massless particles implies the following orthogonality: J µ (p A , p 1 ) · (p A − p 1 ) µ = (−ig) u(p 1 , λ 1 )( p A − p 1 ) g f L P L + g f R P R u(p A , λ A ) = 0. (B.52) Employing this and the longitudinal polarization vector in Eq. (B.20) for V λ V =0 , then the helicity amplitude for the f L → f L V 0 splitting process is (B.56) Importantly, the amplitudes for a transversely polarized V exhibit a dependence on p T , whereas the dependence is on M V for a longitudinally polarized V . As is well-documented throughout the literature, this distinction leads to qualitative differences in scale evolution for λ V = ± and λ V = 0 states. Moreover, as M(f L → f L V 0 ) vanishes as (M V /E V ) → 0, one can interpret the forward emission of V 0 bosons as a "beyond-twist-2" phenomenon. By parity inversion, the helicity amplitudes for the RH f R → f R currents are: Explicit computation reveals that the relative minus sign in the transverse polarization cases can be traced to the relative phase difference that defines the λ V = ± polarization vectors.
(This is also evident by the relative positive sign in the longitudinal case.) At the level of squared MEs (|M| 2 ), the above splitting helicity amplitudes are summarized in Table 10.

B.3 Phase Space Decomposition
For the f (p A ) + B(p B ) → f (p 1 ) + X(p X ) process, where X is some n X -body final-state system, the phase space volume of the (n X + 1)-body system can be organized to isolate the one-body phase space for f . In doing so, we can factorize its contribution from the Table 10. The squared helicity amplitude for the f λ A → f λ1 V λ V process in the collinear limit for each non-vanishing helicity permutation, and assuming the coupling normalization of table 1.
phase space for the hard B + V → X sub-process and remain inclusive with respect to the kinematics of f . Noting the definition of V 's momentum q = p A − p 1 , decomposing the momentum of X over its n X constituents, and expressing the phase space volume for the (p A + p B )-system in terms of the (q + p B )-system, we obtain dP S n X +1 (p A + p B ; p 1 , {p k }) = (2π) 4 In the last line, we exploited the limits of Eq. (B.40), which allows the one-body phase space for f to be written in terms of either evolution variable p 2 T or q 2 . Explicitly, this is d 3 p 1 (2π) 3 2E 1 = (−1) dφ f dz dp 2 T 4(2π) 3 (1 − z) = dφ f dz dq 2 4(2π) 3 . (B.67)

B.4 Transverse Momentum-Dependent Distribution Functions for EW Bosons
To build a set of splitting functions that can be used with both spin-averaged and helicitypolarized MEs, we follow the formalism of Ref. [85] and consider the f + B → f + X scattering process when the external states f and f are in definite helicity states λ A and λ 1 but all other external states are unpolarized. (This implies that helicities are averaged and summed for B and X.) We work in the rest frame of X, which one can eventually identify as the hard-scattering frame. In the collinear approximation, this frame is related to the beam c.m. frame by a boost along the z axis. For massless objects, which travel on the light cone, such longitudinal boosts cannot change the helicity; for massive objects, e.g., EW bosons V , helicity inversions are suppressed by factors of (M V /E V ), or some power thereof, and are strongly suppressed in our kinematic limit. Using the results of App. B.3, the semi-polarized, 2 → n X + 1 cross section is subsequently given by where the color-averaged but helicity-dependent, totally differential cross section is In this expression, the summation runs over all extraneous dof, including color and possible multiplicities of X. N f c and N B c are the color factors for f and B. Working now in the collinear limit, fixing the polarization of V to be λ V , and using the factorized squared ME of Eq. (B.38), we rewrite the helicity-dependent, differential cross section as In reaching the second line we exploited several observations. First is that in the kinematic limit in which we are working, the invariant masses of the (V B)-system (Q) and (f B)-system are related by Q 2 (1 − M 2 V /Q 2 ) = λ 1/2 (Q 2 , M 2 V , 0) ≈ zs AB . Here, λ(x, y, z) is the Källen function defined just below Eq. (2.6). Second is that f → f V splitting is a color-singlet process and that we are free to introduce color-averaging factors of 1/N V c = 1 for the V + B → X sub-process. This also implies that the squared amplitude for f → f V passes through the summation after color indices have been counted (see Eq. (B.63)), i.e., (B.73) Consequentially, we can consistently define the color-averaged but spin-dependent, total and differential cross sections for the (partonic) V + B → X sub-process at Q 2 to bê , and (B.74) (B.75) For the high-energy process µ + µ − → F + X, where the production of F is mediated by high-energy V λ A V λ B scattering, and X is some arbitrary, inclusive final state, we are able to extend the above results and write the scattering formula σ(µ + µ − → F + X) =f ⊗f ⊗σ + Power and Logarithmic Corrections (B.89) (B.90)