Interleaved resonance decays and electroweak radiation in the Vincia parton shower

We propose a framework for high-energy interactions in which resonance decays and electroweak branching processes are interleaved with the QCD evolution in a single common sequence of decreasing resolution scales. The interleaved treatment of resonance decays allows for a new treatment of finite-width effects in parton showers. At scales above their offshellness (i.e., typically Q > Γ ), resonances participate explicitly as incoming and outgoing states in branching processes, while they are effectively “integrated out” of the description at lower scales. We implement this formalism, together with a full set of antenna functions for branching processes involving electroweak (W/Z/H) bosons in the Vincia shower module in Pythia 8.3, and study some of the consequences. Copyright H. Brooks et al. This work is licensed under the Creative Commons Attribution 4.0 International License. Published by the SciPost Foundation. Received 10-09-2021 Accepted 03-03-2022 Published 22-03-2022 Check for updates doi:10.21468/SciPostPhys.12.3.101


Introduction
In an interleaved evolution algorithm [1], two or more distinct types of shower-style resummations are performed jointly, by ordering them in a common measure of resolution scale and combining them into a single "interleaved" sequence of evolution steps. Examples are the interleaving of initial-state radiation (ISR) and multiple parton interactions (MPI) in Pythia 6 [1,2], the subsequent interleaving of final-state radiation (FSR) as well [3] in Pythia 8 [4], and more recently the interleaving of Pythia's MPI treatment with alternative FSR+ISR shower models such as Vincia [5] and Dire [6].
In such a framework, the full perturbative evolution probability can be written schematically as: dP dQ 2 as a function of a generic resolution variable, Q, for which both Pythia and Vincia use measures of transverse momentum, and where the sums run over (systems of) evolving QCD and QED charges respectively. The first line contains the unresummed probability densities (here including MPI as well as the QCD and QED branching kernels) and the second line represents the combined Sudakov-like factor that expresses the total no-evolution probability between a preceding scale Q i−1 and the current scale Q.
Interleaving provides a physically intuitive picture of the fine-graining of event structure with resolution scale. Moreover, to the extent that the different evolution steps can affect each other, e.g. by changing the available phase space and/or the number or character of the evolving entities, the joint resummation will have different properties than what a sequential or factorised approach would produce 1 . While we are not aware of a mathematical proof that interleaved evolution is more accurate than alternative approaches 2 , we believe it is a physically well-motivated paradigm which is worth extending to further physics aspects.
One such aspect which, to our knowledge, has not been considered in an interleaved context until now is that of decays of short-lived resonances (though a formalism with similar effective consequences was put forth by Khoze and Sjöstrand in [9]). In hard processes that involve such resonances, the standard approach in Monte Carlo (MC) event generators is to adopt a factorised approach based on the narrow-width limit. Modulo spin-correlation effects, this is a reasonable starting point, in particular for resonances in the Standard Model (SM), which all have Γ O(GeV) M . Thus, as a first step, one performs ISR+FSR on the Bornlevel hard process treating any outgoing resonances as stable. To account for finite-width effects, one could in principle impose a cutoff or dampening of radiation off resonances at scales of order their widths, to reflect that they cannot physically radiate coherently at wavelengths longer than ∼ ħ hc/Γ , but since this is not a big effect for SM resonances (see e.g., [9][10][11]), such effects are currently not accounted for in Pythia. In a second step, one iteratively treats each resonance decay, adding additional FSR showering to each decay process as needed. In these secondary showers care must be taken to preserve the invariant mass of each resonance-decay system, in order not to generate unphysical corrections to the Breit-Wigner shapes. Again, perturbative finite-width effects are typically regarded as small [9][10][11] and are not formally included in Pythia's current modelling.
In section 2, we extend the concept of interleaving to include decays of short-lived resonances. This opens for a new treatment of finite-width effects, in which resonances can act as physical emitters (and recoilers) during an initial stage of the evolution but effectively are "integrated out" of the physics description at scales below Q ∼ Γ . This automatically suppresses low-frequency radiation off the resonances and allows for low-scale interference effects to generate modifications to the Breit-Wigner resonance shapes. The model has qualitative similarities with (but is not identical to) the preferred scenario of [9]. Moreover we propose a recursive picture which allows for a natural treatment of sequential (cascade) decays.
A related aspect concerns how weak-scale resonances are created in the first place and, more generally, how electroweak (EW) corrections are treated. Pure QED showers aside, the dominant approach is to compute cross sections using fixed-order expansions in the relevant weak couplings. This is a reasonable starting point at relatively low ("electroweak-scale") energies where the integrated probabilities for multiple weak branching processes over phase space, and hence the corresponding electroweak logarithms (see, e.g., [12]), are small. An alternative approach which is better suited for (asymptotically) high energies is that of electroweak showers [13][14][15][16][17], by which electroweak logarithms can be resummed to all orders, at the expense of neglecting subleading-and process-dependent non-logarithmic terms in a similar compromise as that made in ordinary (QCD/QED) showers 3 . The current implementation of EW showers in Pythia [13] focuses on spin-averaged vector-boson (V ∈ [W, Z]) emissions in jets only. Notably, it does not include triple-boson vertices such as V → V V , V → V H, etc., nor does it account for the spin dependence that arises due to the chiral nature of the EW sector. Three alternative comprehensive formalisms have since been developed, in [15][16][17].
In sec. 3, we adapt the EW shower formalism of [16] to the p ⊥ -ordered Vincia shower model and discuss its implementation in terms of antenna functions, evolution variables, Sudakov factors, recoil strategies, and treatment of neutral-boson interference effects. The implementation includes a full set of explicit EW antenna functions with helicity dependence, that are recapitulated in app. A, with overestimating trial integrals for p ⊥ -evolution collected in app. B. The combination of a full-fledged EW shower with other parts of MC event generation raises a number of issues. In sec. 3.5, we discuss how we treat branching processes that are present both in the EW shower and as resonance decays, such as t → bW , Z → qq, etc. Furthermore, in sec. 3.6 we discuss how to avoid double-counting between EW and QCD showers from different Born-level hard processes, such as, e.g., in the case of pp → V V j which can be reached both via a Born-level V V event + a QCD emission, or via a Born-level V j event + an EW V emission.
Extending eq. (1) to include interleaved resonance decays and EW shower branchings, it becomes: where the sum over dP RES /dQ 2 in the second line (and the corresponding negative integral in the first line) expresses the interleaving of resonance decays with the rest of the evolution via probability densities for 1 → N "branchings" (decays), each normalised to integrate to unity, and the full set of EW antenna kernels, P EW , has replaced the corresponding QED ones in the first term. The precise interpretation of eq. (2) will be elaborated on in secs. 2 -3 below. Finally, in sec. 4 we present a set of validations and preliminary results, and use them to discuss the physical implications of our approach, before summarising and concluding in sec. 5. The complete set of EW antenna functions is collected in app. A. Integrals of trial-antenna functions are worked out in app. B. Methods for Breit-Wigner sampling and expressions for the partial widths of Higgs bosons, transversely and longitudinally polarised vector bosons, and top quarks are collected in app. C.

Interleaved Resonance Decays
The starting point for the modelling of resonance production and decays in event generators is the narrow-width limit, Γ /M → 0, also called the pole approximation. In this limit, an infinite timelike interval separates the production and decay of the resonance, and there is no interference between radiation emitted before and after the decay. Formally, the decay of a narrow resonance can be factorised from its production process.
If the resonance carries spin, the factorisation takes the form of a tensor in spin space [19,20]. This can be incorporated in Monte Carlo simulations through a clever linearly scaling algorithm [21][22][23][24][25], as is done for instance in Herwig [24,26], while Pythia employs spinaveraged expressions with correlations imposed on a case by case basis e.g. by using the full 4-fermion matrix element to generate correlated decay angles for the two W bosons in e + e − → W + W − → 4 fermions. In both cases, spin correlations (when they are included) manifest themselves as angular correlations between partons from different decays and/or between partons from a decay and partons in the hard process.
In addition, event generators typically make several finite-width improvements on the strict narrow-width assumption, notably Breit-Wigner distributions BW R (Q 2 ) for resonance invariant masses instead of δ functions, as well as options for allowing partial widths (and hence relative branching fractions) to vary with Q 2 . The latter allows, for example, to account for kinematic thresholds and effects of running couplings across a reasonable range around the pole, but cannot be pushed too far, especially in the electroweak sector where masses and couplings are not independent of each other.
However, the production and decay of the resonance is still treated separately, without accounting for any (perturbative) interference effects between them beyond colour conservation and, in some cases, spin correlations. We refer to such treatments as Breit-Wigner-improved pole approximations (BWPA).
For example, in both Pythia and Herwig, a hard process like g g → tt (with independently selected Breit-Wigner distributed masses for both tops) is first subjected to both initial-and final-state showers starting at the evolution-scale maximum defined by the hard process and ending at the infrared shower cutoff. After this, each of the top-decay processes, t → bW , are subjected to an internal "resonance shower". The latter is done in a way that preserves the invariant mass of the resonance-decay system so that the Breit-Wigner shape of the decaying top quark is preserved (i.e., there are no momentum exchanges with any partons outside of the top-decay system), again only stopping when the infrared shower cutoff is reached. Finally the W decay systems are showered similarly.
The implicit assumption is that interference between radiation emitted in each of these stages (top production, top decay, and W decay) is negligible. The fundamental reason why this is a good assumption, at least for perturbative QCD radiation off SM particles, is that none of the SM resonances (top, Higgs, W , and Z bosons) have widths that are much larger than the shower cutoff for QCD radiation, Q cut ∼ 1 GeV, hence the region of the phase space for perturbative QCD shower evolution over which interference effects could be relevant is very small. The strong suppression of such interference effects have also been verified by explicit theoretical and phenomenological studies e.g. of e + e − → W + W − [27] and e + e − → tt [9,10,28,29].
Nevertheless, the experimentally achievable statistical precision on top-quark mass measurements at hadron colliders has now reached the order of a few hundred MeV [30][31][32][33], making it important to evaluate (and preferably control) QCD uncertainties at that level or better. This has catalysed a reassessment of possible non-perturbative uncertainties such as colour reconnections [34][35][36], and also of the effects of soft perturbative radiation [37,38] and of finite-width effects in fixed-order matrix elements matched to parton showers [39,40]. So far, the latter efforts have focused mainly on improvements to the treatment of finite-width effects on the fixed-order side, and on how to match these consistently with showers, without substantial modifications to the showers themselves.
Here, we note that the BWPA is, strictly speaking, not quite consistent with the strongordering condition in parton showers. Strong ordering expresses the simple fact that the leading singularity structures of QCD (and QED) amplitudes correspond to Feynman diagrams in which each successive propagator has a much smaller virtuality than the preceding one (or next one, for initial-state legs). Physically, this reflects a formation-time principle; short-lived fluctuations do not have time to emit low-frequency radiation. However, for unstable particles in the BWPA, one can have precisely the situation that a particle which has nominally been assigned an invariant mass quite different from the pole value does emit low-frequency radiation. In the corresponding Feynman amplitudes, there are then two (or more) off-shell propagators, which ought to be suppressed relative to amplitudes in which the low-frequency radiation is emitted after the decay. This leads us to consider an interleaved paradigm for showers off resonance-production + decay processes, in which resonance decays are inserted in the overall event evolution when the perturbative evolution scale reaches a value of order the width of the resonance.
The desire to connect with the strong-ordering criterion in the rest of the perturbative evolution, as the principle that should dictate the leading amplitude structures, leads us to prefer a dynamical scale choice for resonance decays, whereby resonances that are highly off shell will persist over shorter intervals in the evolution than ones that are almost on shell. We note that this has the consequence that the on-shell tail will be resolvable by soft photons or gluons, albeit suppressed by the survival fraction. To illustrate this, fig. 1 shows the survivial fractions (denoted ∆ R ) as functions of evolution scale, for t, Z, and W resonances, for three different options for dynamical scale choices, all of which are roughly motivated by the propagator structure: where m 0 is the pole mass and m its BW-distributed counterpart. Near resonance, options i) and ii), illustrated in the left and middle panes of fig. 1, are functionally almost equivalent, differing mainly just by an overall factor 2, while for option iii), illustrated in the rightmost pane, m = m 0 ± Γ /2 translates to Q 2 ∼ m 0 Γ , so that option is primarily intended to give an upper bound on the effect that interleaving could have. Alternatively, our model also allows for using a fixed scale, Q RES ≡ Γ , irrespective of offshellness. In that case, the resonance will not be resolved at all by any photons or gluons with scales Q < Γ . We regard this as a good starting point for the width dependence but have not selected it as our default since the fixed-scale choice by itself does not automatically extend strong ordering to the resonance propagators; this can only be achieved by allowing the choice to be dynamical. Our default choice, eq. (4), is constructed to have a median scale of 〈Q RES 〉 = Γ , while simultaneously respecting strong ordering event by event. This implies that soft quanta will be able to resolve the resonance with a suppressed magnitude ∝ ∆ R , which acts as a form factor.
Finally, interleaved resonance decays, and especially the dynamical scale choices, have the additional benefit of enabling a particularly simple way to match resonance decays and EW branching processes. This is necessary since processes such as t → bW occur not only as resonance decays but also as EW shower branchings (and for sufficiently high energies can occur multiple times); from the point of view of the EW shower, the "decay" process is merely the last such branching, which happens with unit probability. Since the EW shower evolution is ordered in a measure of propagator offshellness, a simple and elegant transition between the EW shower evolution which dominates at very high offshellnesses and the unit-probability Breit-Wigner-distributed resonance decays which dominate near the pole, can be achieved by also ordering the latter in offshellness, viz. by interleaved resonance decays. We return to this in sec. 3.5.

The Model
The change to an interleaved treatment of resonance decays is achieved by introducing a probability density for 1 → n decay(s) in the perturbative evolution, as follows: 1. We assume that we start from a hard process (or a previous decay or EW shower branching process) which has produced a resonance according to the BWPA, i.e. a distribution of the form where m 0 is the pole mass and Γ is the width of the resonance 4 . The overall normalisation factor N ensures that the probability density integrates to unity. In this work, the above formula is also used for resonances produced by the EW shower (see below), with explicit expressions given in app. C.
2. By default, we define the evolution scale associated with the decay of the resonance to be given by eq. (4). As discussed above, this implies, e.g., that the decay of a resonance which has m = m 0 ± Γ /2 will be performed when the interleaved perturbative event evolution reaches a scale Q ∼ Γ , cf. the middle pane of fig. 1. The left-and righthand panes illustrate two alternative dynamical choices, with lower and higher averages, respectively.
3. Having determined a scale Q RES for each resonance in the hard process, the interlaved perturbative event evolution commences, at the shower/MPI starting scale, Q 0 , determined by the hard process.
4. If any Q RES values are greater than Q 0 , the corresponding resonance decays are performed immediately; i.e., those resonances are replaced by their decay products. Any other resonances are kept stable (for now).
5. The interleaved event evolution continues (with shower emissions and multi-parton interactions generated as usual), until either the shower cutoff or a Q RES scale is reached.
When a Q RES scale is reached, the following happens: Figure 2: Recursive interleaving of resonance decays with QED and QCD radiation, for a hard event with t → bW followed by W → qq . The vertical axis represents the evolution scale, from high (top) to low (bottom). The scales m t and m W represent the (BW-distributed) invariant masses of the respective particles, setting upper kinematic limits for each resonance shower. The new aspect is the introduction of the scales Q t→bW and Q W →qq (typically of order Γ t and Γ W respectively), below which each of the resonance-decay systems are merged into their respective production system(s).
6a. The resonance in question is replaced by its decay products. (From the point of view of the rest of the perturbative evolution, this looks like a 1 → n branching.) To give an example, the diagram in fig. 2 illustrates a hard process with an outgoing top quark. Initially, the top quark is treated as a stable outgoing particle which is subject to showering as usual, but when the evolution scale reaches the scale associated to the top-quark decay, here denoted Q t→bW , the evolution of the hard event is put on hold, and the top-decay treatment is started. This step is represented by the leftmost diagonal black arrow.
6b. A "resonance shower" begins. As described in detail for non-interleaved resonance showers in [38], this stage only involves the decaying resonance and its decay products, with the four-momentum and quantum numbers of the initial resonance preserved as boundary conditions. There are no momentum (or quantum number) transfers to/from any partons outside of the resonance-decay system. The starting scale Q 0 for the resonance shower is set to the phase-space maximum for FSR in the given decay process, exactly as in the non-interleaved case. In fig. 2, the change from Q RES = Q t→bW to Q 0 = m t is illustrated by red-dotted lines.
6c. Any new resonances produced by the decay process are associated with their own Q RES scales, which can be higher or lower than Q RES . In the example of fig. 2, the top decay produces a W boson, so a scale for the W boson decay, Q W →qq , is determined.
6d. The resonance shower continues, with FSR shower branchings, until the evolution reaches either the next Q RES , or the original scale Q RES . Fig. 2 exemplifies the former case, with Q W →qq > Q t→bW . (Since Γ W > Γ t in the SM, this is the typical case; it also allows us to illustrate the recursive nature of our approach.) Thus, the scale associated to the W decay, Q W →qq , will be reached before the end of the top-decay resonance-shower stage.
6e. When a nested Q RES is reached, the current resonance shower is placed on hold, while the decay + shower of the corresponding nested resonance decay occurs, recursively, starting back from step 6a. In the example of fig. 2, this is illustrated by the rightmost diagonal black arrow and the subsequent shower evolution in the W → qq system.
6f. When the resonance-shower evolution again reaches Q RES , the resonance-shower stage ends, and the resonance-decay system (including any radiation emitted as well as any nested systems treated so far) is merged with the parton system that produced it. This is illustrated in fig. 2 by the bottom edges of the two resonance-shower systems, at which the showered partons are collected into an effective vertex which is inserted in the parent process. The evolution of the merged system then recommences at step 5, starting from Q RES . Momentum exchanges between partons that originated from the decayed resonance and partons from outside that system may now in principle occur, to the extent prescribed by colour and/or charge flow.
When Q cut is reached, the following happens: 7. As a last step, any resonances with Q RES below Q cut are decayed in no particular order, just as in the conventional (non-interleaved) treatment of resonance decays. However, any nested decays that have Q RES > Q cut will still be performed recursively, as above. For example, the decay of a top quark with Q RES < Q cut will take place after the shower evolution of its production process has finished, but if the W produced in its decay has Q RES > Q cut , then that decay will still be interleaved within the top-decay shower evolution.
Extending eq. (1) to include interleaved resonance decays, it becomes: where it is understood that the ISR+FSR term includes a sum over QED and QCD radiators, and similarly the RES term includes a sum over decayers. Different from conventional interleaved parton-shower and MPI kernels, the term dP RES /dQ 2 is not exponentiated in the Sudakov factor. This is because the probability density expressed by the Breit-Wigner distribution is already unitary and contains its own infiniteorder resummation. In other words: if a resonance is produced, its decay happens once, and once only; there is no need for a Sudakov-style resummation of it. Due to the interleaving with in particular the EW shower, there is, however, a finite probability (given by the EW Sudakov factor) that the resonance will undergo one or more EW branchings before it gets a chance to decay. We return to this in sec. 3.
We emphasise that, at the current stage, this proposal can only be considered a heuristically motivated paradigm. Applying the strong-ordering principle to finite-width propagators produces a kind of forced marriage between two different all-orders summations, the self-energy Breit-Wigner one, and the LL bremsstrahlung one. It captures the basic feature that radiation at frequencies below the resonance width should be suppressed, and we therefore consider it of phenomenological interest to explore its consequences. Should it become relevant to the community, a more formal mathematical investigation would be welcome.
Note also that the systematic inclusion of non-resonant effects would require future extensions of matching strategies, beyond the scope of this paper to explore. A final point left for possible future investigations is that resonances with low off-shellnesses can in principle persist to arbitrarily low scales. This raises the question whether, e.g., top quarks that are assigned off-shellness values less than the infrared shower cutoff (or less than Λ QCD ) should be allowed to hadronise.

Summary of Consequences
To summarise, the main consequences of the interleaving of resonance decays with the rest of the perturbative evolution are: • Due to the interleaving, unstable resonances effectively disappear from the evolution at an average scale Q ∼ Γ . They will therefore not be able to act as emitters or recoilers for radiation below that scale; only their decay products can do that.
• After the resonance has disappeared, recoils to partons originating outside of the decay system are in principle allowed, and may distort the Breit-Wigner shape. In practice, such recoil effects are still expected to be relatively small, for several reasons. Firstly, the fact that the interleaving only "kicks in" below the offshellness scale limits any outof-resonance recoil effects (e.g., in terms of p ⊥ kicks) to be smaller than that scale. Secondly, in decays of QCD colour singlets, such as Z and W bosons, there are no leading-colour (LC) dipoles to any partons outside of the decay system and hence no out-of-resonance QCD recoils at all. Even top-quark decays only involve a single such connection, corresponding to the colour flowing through the decay. Analogous arguments also apply to QED radiation, with α s → α EM and the colour of the resonance replaced by its overall electric charge.
• With the dynamical choice of decay scale, highly off-shell particles disappear from the evolution at higher evolution scales than ones nearer the pole mass value, translating to an increasing distortion of the resonance shape further away from the pole. This roughly corresponds to the notion of strong ordering in the rest of the evolution.
At the technical level, the interleaved handling of resonance showers does not dramatically change the amount of time it takes to generate events. To illustrate this, we consider fully showered and hadronised e + e − → tt events at s = 500 GeV, taking the case of sequential decays with Γ t = 1.5 GeV as reference. Fig. 3 shows the relative change in run time for the fixed-scale and default dynamical-scale choices, as functions of the top quark width, for the case of stable W bosons (left plot) and with W decays included (right plot). Note that there is an overall uncertainty on the time measurements of a few per-cent.
Without W decays (left plot), the total run-time is nearly the same with and without interleaving for the SM value of the Γ t , while interleaving actually decreases the run-time for large widths. We interpret this as due to the non-interleaved treatment having to perform two evolutions in the range Γ t > Q > Q cut , while the interleaved treatment only goes through that region once, removing the double-counting. In the right plot, with W decays included, the interleaved treatment is 10% -20% slower than the sequential one. (Note: the W width was not varied, only the top-quark one.) For pp → tt, interleaving does not measurably change the total run time at all, which is dominated by the generation of multi-parton interactions, colour reconnections, and hadronization.   Figure 3: Relative time consumption to generate e + e − → tt events at s = 500 GeV, as functions of the top quark width, for sequential (blue) and interleaved resonance decays with fixed (red) and the default dynamical (yellow) scale choices. Left: with stable W bosons. Right: with inclusively decaying W bosons. In both cases, sequential decays with Γ t = 1.5 GeV was used as the reference case. The run-time for the reference sample without W decays was 1.7 ms/event, increasing to 3.5 ms/event with W decays included. Tests performed on a single 2.7 GHz Intel Core i7 CPU, with clang version 13.0.0 with -O2 optimisation.

Electroweak Showers
In this section, we discuss the implementation of electroweak radiation in the Vincia parton shower. The realisation in Vincia draws heavily from the formalism set out in [16]. We provide a brief summary of the common points here and discuss the adjustments that have been made to assimilate it with the Vincia QCD shower. A comprehensive description of the QCD shower, including details like its antenna functions, exact phase-space factorisation and kinematic maps may be found in [5,41,42].
Vincia's QCD shower is based on the antenna subtraction formalism [43,44] and allows for the evolution of states with definite helicity [42,45,46]. This property is especially important in the electroweak sector due to its chiral nature [15,16]. However, it does not equate to a complete treatment of spin interference effects. Their inclusion in parton showers has received some attention recently [24,25]. However, an extension of the algorithms described there to the EW sector would not be straightforward. Spin does not affect emission rates in QCD, but only causes modulation in azimuthal distributions. Practically, this means that one does not need to modify the Sudakov form factor when applying corrections due to spin correlations. This is no longer the case in the EW sector, where the vector boson couplings depend on the spin of the branching particle. The intermediate solution of evolving definite helicities is able to describe such effects, while not requiring an increase in the algorithmic complexity of the shower.

Antenna Functions
To facilitate evolution with polarised particles, branching kernels for every helicity configuration are required. The kernels described in [16] are computed with the spinor-helicity formalism [47], similar to those used in [48,49] to compute QCD antenna functions. While these branching kernels capture the appropriate quasi-collinear singular structure of electroweak branchings, they suffer from numerical issues away from the singular limits due to a mismatch with the Vincia shower variables. In the new implementation, they have thus been replaced by antenna functions, bringing the formalism further in line with the QCD sector. While the electroweak antenna functions are defined in terms of the same phase-space invariants as those of the QCD sector, we currently make no attempt to enforce soft coherence in the electroweak sector. In [50], the complete multipole structure of photon emissions was included in a dedicated QED shower without spin-dependence. An extension to the EW sector requires a generalisation of that algorithm to include all EW couplings and helicity dependence, which we reserve for future work.
The calculation of the antenna functions for all final-state EW branchings, as well as vector boson emissions off fermions in the initial state is detailed in A. Crucially, the formalism ensures that its polarisations match up with those used by the HELAS [51] routines implemented in MadGraph5 [52] which Vincia uses to polarise the hard scattering. Special care has to be taken to avoid gauge-dependent relics associated with the Goldstone component of the longitudinal polarisation, which would cause unphysical energy scaling that would violate unitarity if not treated correctly [15][16][17]. The antenna functions calculated here have been validated to produce the same collinear limits as those of [16].

Evolution Variables and Overestimate Determination
In correspondence with Vincia's QCD shower, pre-branching momenta are denoted by capital letters and post-branching momenta by lowercase letters. The EW shower includes emissions from two final-state momenta, I, K → i, j, k, and from two initial-state momenta, A, B → a, j, b. Other types of emissions, such as those off one initial-state momentum and one final-state momentum, are not included. The current implementation does not attempt to correctly model coherent weak gauge boson emissions, and as such all relevant singularities can be incorporated with the above two types of branchings. The shower formalism is set up in terms of Lorentz-invariant quantities, e.g. s i j = 2p i ·p j . It is convenient to define The ordering scales are then given by where in practice all initial state momenta are treated as massless. The Sudakov form factor, i.e. the probability that no EW branching occurs between scales Q 2 n > Q 2 n+1 , may then be written as where Here, α(Q 2 ) is the running electroweak coupling evaluated at first order and a i are the antenna functions defined in app. A, written in terms of the ordering scale Q 2 and an auxiliary variable ζ. The PDF ratio R f is given by Finally, dΦ ant is the parton-shower component of the antenna phase-space factorisation as defined in app. B. The shower distributes its EW branchings according to the probability distribution Practically, this is sampled using the Sudakov veto algorithm [53][54][55], where the individual components inside the evolution integral eq. (11) are overestimated by strictly larger and simpler expressions, and branchings get accepted with probability whereα is a constant overestimate of the coupling constant,R f is a constant overestimate of the PDF ratio and a i,trial (Q 2 , ζ) is a simple overestimate of the antenna functions a i (Q 2 , ζ). The coupling constant α(Q 2 ) < α(Q 2 n ), the latter of which functions as the constant overestimate. Like the full antenna functions, the trial antennae are defined in terms of the Lorentz-invariant dot products. Due to the vast number of possible branchings in the electroweak sector that all come with varying masses, couplings and helicity dependence, the process of determining trial antennae is automated. To that end, parameterised overestimates are defined as where x i and x j are collinear momentum fractions defined in app. A. Note that in the finalfinal case, a term proportional with the mass m I appears, but terms proportional with m i and m j are absent. The latter terms are not required because the associated mass corrections in the antenna functions tend to be negative, while those that scale with m I tend to be positive and become dominant when Q 2 ≈ m 2 I . The trial antenna for initial-initial branchings consists of a single term because it only has to model vector boson emissions. While appropriate values for c II are easily computed for all flavor and helicity configurations, values of the coefficients c FF 1 through c FF 4 are determined automatically. A large number of branchings is sampled from randomly sampled antenna masses, and the corresponding antenna function is determined for each of them. A suitable trial antenna may then be found by minimising the average distance between the trial-and real antennae while ensuring that for no point in the sample a FF i,trial < a FF i . Such optimisation procedures are instances of linear programming, for which many external solvers are available. We make use of the Python package PuLP [56]. The trial evolution integrals eq. (11) are worked out in app. B.

Recoiler Selection
While in Vincia's QCD shower antenna are naturally spanned between colour-connected partons, no such guiding principle exists in the electroweak sector. In fact, since the electroweak shower makes no attempt to incorporate soft interference effects, the second particle only acts as a recoiler to the collinear emitter. In [16] it was shown that through a probabilistic choice of recoiler selection the effects of recoil of previous branchings may be partially mitigated. In [15], a similar modification was included as a multiplicative factor on the branching kernel. We adopt the strategy of probabilistic choice here, where the probability to select a final-state recoiler K for final-state particle I may be written as where the functions A(p i , p j ) are the collinear helicity-dependent 1 → 2 branching amplitudes computed in [16]. The sum over K in the denominator runs over all recoiler candidates, and the sums over X and X run over all electroweak clusterings of particles I and K or K respectively. The probability given by eq. (16) gives preference to recoilers K that were most likely to have emitted I. Because the final-final kinematic mapping preserves the total antenna momentum, the propagator of the corresponding particle X then remains preserved. Initial state branchings are always set to recoil against the other initial state.

Bosonic Interference
The electroweak sector has the unique feature of introducing interference effects between neutral vector bosons [15,16,57]. In [15] a physically accurate but computationally expensive treatment was outlined using evolution of mixed states in density matrices instead of a single definite event. We instead opt for a much simpler treatment, assigning an event weight to the event every time a neutral boson disappears through splitting to fermions or W bosons. For post-branching momenta p i and p j , the weight associated with the interference between neutral bosons V 1 and V 2 is given by where p i j = p i + p j . The pair V 1 and V 2 corresponds with either a photon and a transversely polarised Z boson, or a Higgs and a longitudinally polarised Z boson. Interference between spin states is not included in correspondence with the rest of the shower formalism. The sum over X runs over all particles that may have emitted the neutral vector bosons through the collinear branchings X → X V 1 and X → X V 2 .

Resonance Matching
The electroweak shower involves a number of branchings that would normally be associated with resonance decays, such as t → bW , Z → qq, etc. In such cases, the EW shower correctly describes the collinear dynamics of these branchings at off-shellness scales that are much larger than the EW scale Q 2 EW . At such scales, the resonance masses are subleading corrections to the regular collinear factorisation dynamics. However, at scales close to the resonance width, resonance decays are more correctly described by a Breit-Wigner distribution. Pythia's default method of treating resonance decays is to sample the resonance mass from a Breit-Wigner distribution upon production. Vincia follows this procedure, making use of helicity-and massdependent decay widths that contain O(α s ) corrections. These widths, as well as the sampling procedure, is detailed in app. C.
Resonances may be produced as part of the hard scattering, or by the EW shower. In either case, the EW shower may sample resonance-like branchings at scales sufficiently above the EW scale. The Breit-Wigner-sampled masses on average differ from the on-shell mass only by O(Γ ), and thus do not affect the EW shower significantly. The resonance decay-like branchings are suppressed by a factor Q 4 /(Q 2 + Q 2 EW ) 2 , such that the shower branching probability smoothly vanishes at scales close to Q 2 EW . Then, when the shower evolution reaches the scale associated with the offshellness of the resonance, given by eq. (4), without branching it, the resonance is decayed according to the procedure outlined in app. C. This ensures a smooth matching between the two descriptions, as is illustrated in fig. 4 for the case of a left-handed top decay. The value of Q 2 EW is implemented as a tunable parameter in this context. All high-energy SM resonances produced by the EW shower show the behaviour illustrated in fig. 4, where the EW shower dominates the Breit-Wigner shape at high invariant mass. This means that the exact treatment of running width effects far away from the peak, which are prone to gauge violations, are of little practical importance. The interleaving procedure described in the previous section is available for all resonance decays, whether they are produced or decayed by the EW shower or otherwise.

Overlap Veto
When the parton shower generates both QCD and EW branchings, care has to be taken to avoid double-counting the approximated matrix element of certain topologies. This doublecounting may occur when different Born-level events can populate the same phase-space point after shower branchings. For example, a final state V V j may be reached by starting from either a V V born-level state and adding a QCD initial-state branching, or by starting from a V j bornlevel state, and adding an EW branching. A solution to this issue was described in [13] and is included in the Pythia EW shower in the case of a single vector boson. Here, we generalise this procedure to be universally applicable.
Double-counting can be avoided by ensuring that the QCD and the EW showers only populate the regions of phase space that they are most accurate in. For instance, in the previously mentioned case of V V j where the two vector bosons are collinear, the path of a V j + an EW emission should be preferred as the EW shower models collinear vector boson emissions accurately. On the other hand, if the jet is soft, the preferred path is V V + a QCD emission off the initial state. The above procedure can be implemented generally by introducing a resolution measure that indicates which shower most accurately describes a particular phase space point. Note that this choice of resolution measure is not unique, and any choice that correctly separates the singular regions describes by both showers will suffice. Different choices will lead to different behaviour in intermediate regions where neither shower provides an accurate description.
Here, we choose to follow [13] and make use of a resolution measure inspired by the one used by the k t jet algorithm [58], generalised to account for particle masses, defined as where p 2 t = k 2 t + |m 2 i + m 2 j − m 2 I |, k t is the transverse momentum, ∆ i j is the angular distance between i and j. The parameter R determines the relative weight associated with initial-state and final-state emissions, again affecting behaviour in the non-singular intermediate regions.
It is left as a tunable parameter with default value equal to the implementation of [13]. The mass corrections are included to more accurately reflect Vincia's ordering scale. For instance, in a branching like W ± → W ± Z, the appropriate measure is p 2 t = k 2 t + m 2 z . The absolute value is added to prevent the measure from dropping below zero. It is only relevant in resonancelike branchings where m I > m i , m j , in which case no overlap between the QCD and EW sector exists anyway.
The resolution measure d i j is small in all singular regions of phase space for both QCD and EW emissions, and can thus be used to define a veto procedure. More specifically, if, for instance, the shower generates a QCD emission, eq. (18) is evaluated for that branching, as well as for all possible 2 → 1 EW clusterings in the post-branching event. Then, the QCD branching is accepted if d , and vetoed otherwise. The reverse procedure is also applied if an EW branching is generated. This procedure effectively sectorises the QCD and EW showers, ensuring that every phase-space point is only populated by either a QCD or EW shower history. Note that such a sectorization is also relevant for the construction of shower histories in the context of matrix element merging. Vincia's QCD shower is itself sectorized [42], ensuring only a single QCD shower history is associated with every phase space point. In a future version, the EW shower may also be sectorized, leading to a single shower history through both showers in combination with the overlap veto.

Validations, Results and Discussion
In this section, we present several validations and results concerning both the implementation of interleaved resonance decays, as well as the EW shower.

Interleaved Resonance Decays in ee → tt
We take the Vincia sector-antenna shower [42] implemented in Pythia 8.306 [4] as our baseline for illustrating the effects of interleaving. (For the non-interleaved case, a detailed study of Vincia's treatment of radiation in top-quark decay can be found in [38].) In all of the figures presented in this section, we take the conventional (sequential) treatment of resonance decays as our baseline (blue), and compare with the new interleaved method with a fixed scale equal to the width Γ 0 (red), or a dynamic scale choice, P(m), given by the inverse-propagator distribution, eq. (4) (yellow). We note that the latter treatment is now the default in Vincia since Pythia version 8.304 (while Pythia's simple showers retain the sequential treatment as default).  As a first simple test case, we consider e + e − → tt at s = 500 GeV, for the same mass values as were used in [35], m t = 173.3 GeV and m W = 80.385 GeV. To focus on the radiation emitted from the top quarks (and their decays), we keep the W bosons stable, and QED bremsstrahlung from the incoming beams is switched off.
In fig. 5, we show the spectra of FSR branchings as a function of ξ = ln p 2 ⊥ , where p ⊥ is Vincia's shower evolution scale. For reference, ξ(m t ) = 10.3, ξ(20 GeV) = 6, ξ(2Γ 0 ) = 2.2, and ξ(Γ 0 ) = 0.8, with the top-quark width Γ 0 = 1.5 GeV. The upper pane includes all branchings, regardless of where in the shower they occur, while the three lower panes separate the contributions from (left) emissions off top quarks, (middle) branchings during the top-decay resonance showers, and (right) all other branchings. Starting with the lower left-hand pane of fig. 5, the rate of emissions coming directly from t and/ort quarks is unchanged for large ξ (where this is the dominant component of the spectrum). At low ξ, the fixed scale choice produces a sharp cutoff at Γ 0 , while the dynamic choice produces a tapering off of the emission rate as the top quarks gradually disappear from the evolution. In the middle pane, we see similar effects inside the top-quark resonance showers. In the sequential case, these continue all the way to the QCD shower cutoff, while they terminate at a higher (fixed or dynamic) scale in the interleaved cases. In other words, the first two lower panes of fig. 5 illustrate that, in the interleaved treatment the top quarks themselves (and the component of their decay systems that is still treated as factorised from the rest of the event) are only relevant at evolution scales greater than their offshellnesses.
The rightmost lower pane of fig. 5 shows that this difference is almost perfectly made up for by an increase in radiation not associated directly with the top quarks (nor with their resonance showers). In the interleaved case, this is the component that dominates the evolution at low ξ, while in the sequential case those branchings still take place within the top resonance showers. Interestingly, for the dynamic scale choice, the increase in this component persists up to quite high scales. For reference, ξ = 5.4 corresponds to a scale of order ten times the width. This is due to the constraint of preserving the invariant mass of the top-decay system being lifted quite early on in the evolution for highly off-shell tops and the opening up of a larger phase space, especially relevant on the low-mass tail.
Although the relative contributions of each of the three classes of branchings shown in the left-hand pane thus change by substantial amounts, it is reassuring to note that, at the summed level, the total emission rate illustrated in the main (upper) pane of fig. 5, only changes by small amounts. We note that the total emission rate of very soft emissions, for scales below ∼ 1 GeV (ξ ∼ 0) at the extreme left-hand edge of the plot, drops a little. We regard this as a reasonable physical consequence of preventing the top-production and top-decay systems from both radiating at scales below the width; now only the decay products can do that. With the dynamic scale choice, there is moreover also a ∼ 10% enhancement of radiation at scales up to an order of magnitude larger than the width, due to the aforementioned opening-up of the phase space after exiting the resonance shower. We believe this is an interesting effect which exemplifies a non-trivial interplay between the two resummations.
The evolution-scale spectrum itself is of course model dependent and not physically observable. We show it mainly to illustrate the scale progression in the underlying algorithm. with |m − m 0 | < Γ 0 /2 (corresponding to about half of all top quarks). The emission rates are all normalised to the total cross section in their respective regions. The dominant feature is that interleaving suppresses radiation at energies at and just below the width. For the dynamic choice of interleaving scale, this is more pronounced for off-shell tops (left and right panes) than for on-shell ones (middle pane), while for the fixed scale choice it is independent of the offshellness, as expected. Interestingly, the fixed scale choice in particular also results in a small increase in radiation at very low scales (ξ ∼ −2 corresponding to E ∼ m π ) though it is doubtful whether this would have any observable consequences.
Turning to hadron-level quantities, we first consider two observables proposed in [9], defined by clustering all particles in the event (omitting any that originate from the W decays) into exactly two jets (using the e + e − k t algorithm [59] via FastJet [60]), which are stand-ins for the Born-level b quarks. This is most appropriate near the tt threshold where there is no radiation from the tt pair before decay. In that case, neglecting spin correlations between the two decays, the two jets should be distributed independently and isotropically, i.e. with a uniform distribution in cos θ J J . In our case, we are considering a situation somewhat above threshold, with a CM energy of 500 GeV corresponding to a Lorentz boost factor for each of the top quarks of γ = 1.45. The cos θ J J distributions shown in the left-hand pane of fig. 7 are therefore peaked towards -1, but there is still an interesting shift towards less strong anticollimation when interleaving is enabled, which is especially pronounced in the (suppressed) region where the two jets end up with a relative angle of less than 60 degrees (cos θ J J > 0.5). This agrees qualitatively with the conclusions of [9].
In the right-hand pane of fig. 7, we show how the average charged multiplicity depends on cos θ J J . Not surprisingly, the largest multiplicities occur when the two jets are approximately back to back, while fewer particles are produced when the two jets are closer to each other. This trend becomes stronger when interleaving is enabled, again in line with what was found in [9].
We round off our discussion of e + e − → tt by considering the hadron-level energy distribution as a function of the angle from the b quark (which we take as a proxy for the b-jet). This is somewhat motivated by the studies in [10], which however used the top direction as the   Figure 9: Hadron-level energy distribution in leptonic e + e − → tt events at s = 500 GeV, as a function of the angle from the b jet, for (upper pane) events in which the b quark from the top-quark decay was emitted "forwards" (with cos θ t b > 0.9) and (lower panes) at right angles (with sin θ t b > 0.9), with the latter divided into two hemispheres, cf. fig. 8. reference. We use the b-quark direction instead since we are still relatively close to threshold and this places the peak of the main observed jet at a well-defined location in the plots. To focus on the top-associated particle production, we consider fully leptonic events, and remove the leptons from the W decays from the energy distributions. The coordinate system is illustrated in fig. 8. We first consider the case when the angle between the t and b directions, θ t b , is small, cos θ t b > 0.9. Then, the t and b directions are approximately degenerate, and thet is at 180 • . The energy distribution for these events, as a function of the angle to the b quark, is shown in the left-hand pane of fig. 9. In the interleaved cases, there is a slight (∼ 5%) increase in the energy deposited in the outer regions of the b jet, relative to the sequential decay.
In the middle and right-hand panes of fig. 9, we consider events with a large angle between the t and b directions, sin θ t b > 0.9. We divide these events into two hemispheres, as illustrated in fig. 8. In the middle pane, we show the angular distribution in thet hemisphere, which contains the flight direction of thet. This hemisphere thus tends to contain theb jet, with is smeared over a large angular region due to varying event kinematics and the fact that thet is not highly boosted at s = 500 GeV. No significant effects of interleaving are evident in this hemisphere. The t hemisphere, on the other hand, is relatively free from contamination of theb jet, and here we see a small suppression of the energy density at basically all angles.
Thus, we find that there can be kinematics-dependent modifications of up to ∼ 5% to the energy flow in the event, here illustrated taking the b-quark direction as reference axis.

Interleaved Resonance Decays in pp → tt
As our final example of the effects of interleaved resonance decays, we take pp → tt at s = 8 TeV, repeating the analysis of semileptonic events that was used to estimate colourreconnection effects on the hadronically reconstructed top-quark mass in [35].
In the left-hand pane of fig. 10, we show the multiplicity of central charged particles, with |η| < 2.5 and p ⊥ > 0.5 GeV. We see that both of the interleaved options produce a slight reduction in the multiplicity of charged tracks relative to the sequential decay treatment. The multiplicity of anti-k t jets [61] with R = 0.5 and p ⊥ > 20 GeV, shown in the right-hand pane, is not affected.
The result of the primitive hadronic W and t mass reconstruction of [35] is shown in fig. 11. The anti-k t jets found above are combined to form jet pairs which are accepted as W candidates if they have a combined invariant mass in the vicinity of the W mass (|m j j − m pole W | < 5 GeV). The best such W candidate is then combined with a third jet and accepted as a top candidate if |m j j j − m pole t | < 20 GeV. The difference between the invariant mass of the best W candidates in each event and the W pole mass value is shown in the left-hand pane of fig. 11, with interleaving producing no significant differences in the distribution. The slight increase at higher masses translates to an increase in the hadronically reconstructed W mass of 35 MeV (with a statistical MC uncertainty of 6 MeV) for both of the interleaving options.
The right-hand pane of fig. 11 shows difference between the hadronically reconstructed top-quark mass and the pole-mass value. Again, reassuringly, there is essentially no difference visible at the resolution scale of the plot. Nevertheless, the directly calculated mean of the ∆m t distribution decreases by about 140 MeV (with a statistical MC uncertainty of 25 MeV) when interleaving is switched on (contrary to the W mass which changed in the opposite direction). This change is somewhat larger than the 30-MeV figure obtained for e + e − → tt in [9], and may have relevance for high-precision top-quark mass measurements at LHC. A more elaborate assessments with full-fledged top-mass extraction methods, including the effects of in-situ calibration, is beyond the scope of this work but would be well-motivated.

Validation of the EW Overlap Veto
In this section, we validate the implementation of the overlap veto procedure, as well as the correct implementation of triple-vector boson interactions in the EW shower. To that end, we compare the direct matrix element for pp → V V j at 14 TeV with the shower prediction. Here, the vector bosons include Z, W + and W − and the jet is a gluon or a light quark. Fig. 12 shows the comparison of the leading-order matrix element with several shower histories, as a function of the angular separation between the vector bosons ∆R V V . This observable serves as a natural separator between the regions of phase space where the individual shower contributions should approximate the matrix element well. For small angular separations, the vector bosons are collinear and the EW shower should perform best, while at large angular separation the vector bosons are back-to-back, and the QCD shower should be preferred. To allow for sufficient phase space for the EW shower to radiate in, the hard scattering is ensured to be highly energetic by requiring 0.5 TeV < p ⊥,jet < 1 TeV. The dashed lines show the individual contributions of pp → V V with a QCD emission and pp → V j with an EW emission, and the solid blue line shows their sum. When the overlap veto is not enabled, the sum of contributions clearly overshoots the exact matrix element significantly. On the other hand, when the overlap veto is enabled, the summed Vincia shower prediction matches the matrix element very well and the showers correctly cover their associated regions of phase space.

Fragmentation in Heavy Dark-Matter Decays
As another test of the EW shower implementation, we consider the computation of the prompt decay spectra of heavy dark matter that decays to Standard Model particles. If such dark matter particles are indeed heavy enough, an EW shower develops as it decays, and the shower products may appear on earth in the form of cosmic rays. Recently, these decay spectra were computed by Bauer, Rodd and Webber in [62], expanding upon previous work [63,64]. Their methods involve a numerical evolution of DGLAP equations above the EW scale in the unbroken phase of the Standard Model. The result is then matched to the broken phase at the EW scale, after which the rest of the evolution is performed with Pythia, until only stable particles S = {γ, e ± , p ± , ν e,µ,τ ,ν e,µ,τ } remain. With the inclusion of all possible collinear branchings in Vincia's EW shower, the same physics content should be included. However, the methods of [62] allow for evolution up to dark matter masses up to the Planck scale, while Pythia and Vincia can go up to O(100 TeV) before issues related to numerical precision start to appear. Fig. 13 shows a comparison between their results and the Vincia and Pythia predictions for the energy spectrum of a number of stable particles in the decay of a 200 TeV dark matter particle to electron-neutrinos. Already at this energy, the difference between the Pythia prediction and the two other results is apparent. It is again caused by the missing physics content of Pythia's EW shower, which includes the absence of triple-vector boson interactions and a treatment of spin. The difference is particularly striking in the hard photon spectrum, where the Pythia prediction drops off rapidly while the other two lines show a characteristic bump. Another striking feature is the fact that Vincia shows consistently relatively decreasing soft spectra, while Pythia appears to match up with the results of [62] up to a vertical shift. This is likely caused by the fact that [62] use Pythia to do the low-scale evolution, while in the former case the complete evolution is performed by Vincia. Furthermore, there are significant differences between the treatment of [62] and Vincia's EW shower. These include treatment of spin interference, a different treatment of soft interference and the matching procedure at the EW scale, which is not required in Vincia as it performs all evolution in the broken description of the Standard Model.

Electroweak Sudakov Logarithms
The previous two sections have illustrated the ability of the EW shower to correctly incorporate real EW corrections. In this section, we instead focus on the associated virtual corrections. Virtual EW Sudakov logarithms are known to have large effect in the hard tails of observables already at the LHC, and in particular at future colliders. While such corrections are regulated by the EW scale, they are physically relevant without the inclusion of the associated real corrections, because those lead to different experimental signatures. Much work has already been done on the analytic calculation and resummation of these corrections [65][66][67][68][69][70][71][72][73]. Here, we illustrate that they may also be calculated with Vincia's EW shower. Fig. 14 illustrates the size of the negative virtual EW corrections in the process pp → Z Z → e + e − µ + µ − at 14 TeV and 100 TeV as a function of p ⊥,Z . The EW shower incorporates these virtual corrections through its unitary nature. When it is enabled, some Z bosons will, for instance, branch to W + W − , or another Z may be radiated from the initial state. As a result, events without such corrections effectively get weighted by the EW nobranching probability where the virtual corrections are exponentiated. Fig. 14 shows these virtual corrections at the Monte Carlo level, as well as at the fiducial level, where the cuts 65 GeV < m Z < 115 GeV, p ⊥,l > 25 GeV and |η l | < 3.5 (19) are applied, to emphasise that they arise as a result of the different experimental signatures of the real corrections. At 14 TeV, the corrections reach −30% towards the highest end of the spectrum, while they go down to −70% at 100 TeV. Fig. 15 compares the run time penalty as a function of p ⊥,Z for activating the EW shower in pp → Z Z → e + e − µ + µ − at 100 TeV. As might be expected, this penalty increases with p ⊥,Z up to a factor of two as a larger phase space opens up for EW radiation in the final state, the products of which might in turn induce further activity.

Summary and Conclusions
We have presented two interrelated extensions of the Vincia shower framework in Pythia 8, introducing interleaved resonance decays and electroweak shower branchings, respectively, in the perturbative evolution. The latter is based on the formalism presented by one of us in [16], while the former is a new proposal that shares some features with earlier work by Khoze and Sjöstrand [9].   Figure 15: Relative time consumption to generate pp → Z Z → e + e − µ + µ − at 100 TeV as a function of p ⊥,Z , without EW shower (blue) and with EW shower (red). Tests were performed on a single M1 Pro core with clang version 13.0.0 with -O2 optimisation.
When decays are interleaved with the final-state shower evolution, each unstable resonance is assigned an offshellness scale, which, in our implementation, can either be chosen to be static (e.g., the width, Γ 0 ) or dynamic (e.g., a Breit-Wigner distributed offshellness). When the overall evolution reaches this scale, the resonance is removed from the event and replaced by its decay products; these are then subjected to a "resonance shower" which conserves the total 4-momentum of the decay system and thereby preserves the shape of the original (Breit-Wigner) resonance-mass distribution. This shower starts at the mass of the decaying resonance and ends when the evolution again reaches the offshellness scale, after which the resonancedecay plus shower products are merged with the production process and the evolution is continued starting from the offshellness scale. In this last step, any resonance-final antennae [38] are replaced by corresponding (non-resonant) final-final ones, which allows for the possibility for further radiation to generate out-of-resonance recoils which can distort the Breit-Wigner shape. The procedure also results in a suppression of low-frequency radiation from short-lived particles. We regard these features as making good physical sense and in line with conclusions of previous studies [9][10][11]27,29]. Further properties of our model include that cascade decays can be nested recursively when successive offshellness scales are increasing, and the virtue of enabling a robust and relatively simple matching between EW shower branchings and EW resonance decays.
As a complement to the more sophisticated treatment of interleaved resonance decays, the implementation of a full-fledged EW shower is now part of Vincia. It relies heavily on the helicity shower formalism set out in [45,46] to model the chiral nature of the EW sector. The EW shower includes all possible EW collinear splittings, including those originating from triplevector boson interactions or coupling to the Higgs. These were notably previously missing from Pythia's default EW shower [13]. A number of interesting features appear in the EW sector, including neutral boson interference effects and the matching of resonance-like branching to a Breit-Wigner distribution. Our implementation ensures that such branchings that occur at large offshellness are modelled by the EW shower, while mirroring Pythia's Breit-Wignerbased treatment at small offshellness. Note that the interleaving of resonance decays becomes especially relevant at the large offshellness values probed by the EW shower. Finally, when both the QCD and EW showers are enabled, a risk of double-counting appears due to the possibility of reaching identical states from different Born configurations. To avoid this, a generallyapplicable veto procedure was introduced that avoids any overlap and ensures every phase space point is populated by the most accurate shower.
The effects of interleaved resonance decays are particularly relevant for coloured resonances, such as top quarks, but also apply in the QED sector, to electrically charged resonances such as W bosons. In the latter context we note that Vincia has a fully coherent treatment of QED radiation, including multipole interference effects [50,74]. It is, however, not yet possible to enable the new EW shower at the same time as the fully coherent QED module, as it does not support treatment of helicity-dependent photon emissions. One currently has to choose which functionality is most important for the study at hand, with the coherent QED module being the one that is enabled by default.
Case studies in the context of tt production at e + e − and pp colliders were presented in sec. 4. For the most part, these validated the expectation of relatively modest effects, due to the comparative smallness of the SM top-quark width. Nevertheless, we found a shift of of about 140 MeV on a (primitive) reconstruction of the invariant mass of hadronically decaying tops in pp collisions, which may be worth investigating further, with more full-fledged analysis techniques.
Sec. 4 also includes a validation of the overlap veto of the EW shower, where the V V j matrix element is compared to the shower approximations V V + QCD and V j + EW. The results show that the EW shower correctly includes the triple-vector boson interactions, and that the overlap veto is necessary to avoid overshooting the three-body matrix element. Furthermore, a comparison to the work of [62] was shown, where the energy spectra of the decay high-mass dark matter to SM particles was computed analytically. We found that, while there are still major differences in the modelling of the physics, Vincia's EW shower closely matches their results. Finally, the size of virtual EW Sudakov logarithms in pp → Z Z were investigated at the LHC, as well as at future collider energies. The results show that, as has long been known [65][66][67][68][69][70][71][72][73], these effects are especially large at future collider energies, and an EW shower includes them systematically.
Potential avenues for future research in the area of interleaved resonance decays include studying the formal interplay between the Breit-Wigner and parton-shower resummations and the modelling non-perturbative effects in the tail of long-lived coloured resonances. On the EW shower side, possible ways forward are the enabling of coherent QED radiation in conjunction with the EW shower and especially assessing its logarithmic accuracy with the modern techniques explored in [75][76][77] and comparing to the work of e.g. [78,79]. The coupling constants are defined in Table 1. There, v and a are the vectorial and axial couplings of ff V vertices, while g v and g h are the couplings of V V V and V V H vertices. The calculation of the antenna functions then proceeds by computing the relevant 1 → 2 collinear branching in the collinear limit, taking care to remove the gauge-dependent spurious terms associated with longitudinal polarisations as described in [15][16][17]. After taking the quasicollinear limit, all reference vectors defined in eq. (22) become identical, and the branching processes can be expressed in Lorentz-invariant inner products that involve only the branching momenta and a common reference vector r. These are finally replaced by momentum-fraction variables that are defined in terms of the Vincia phase-space parameterisation.

A.1 Final-Final Antennae
are made, where m 2 I K = (p I + p K ) 2 and x i and x j are collinear momentum fractions defined in terms of the Vincia phase-space variables. To further shorten notation, it is convenient to definem Wherever appropriate, quark-flavour off-diagonal branchings are also included, weighted with the appropriate CKM matrix entry.

A.1.1 Vector Boson Emission off (Anti)fermions
The antenna for vector boson emission off antifermions are identical with the exception of the exchange (v ± λa) ↔ (v ∓ λa).  x

A.1.6 Transverse Vector Boson Splitting to Two Vector Bosons
x i x j .

A.1.7 Longitudinal Vector Boson Splitting to Two Vector Bosons
in the limit m i = m j = 0, and where we have used that for all electroweak branchings m k = m K . Since the massive phase space is contained in the massless one, the requirement of the positivity of the three-body Gram determinant may be incorporated with a veto. Transforming to transverse momentum and the above definitions of the auxiliary variable, we find where the sum runs over all decay channels. Because the electroweak shower produces resonances with definite helicity states, the partial widths should also be computed as such. We make use of the variables The partial widths are given by where N c = 3 for decays to quarks and N c = 1 for decays to leptons. These widths include full mass corrections, as well as O(α s ) corrections in accordance with Pythia [2]. Note the appearance of y 0 in the top width. It appears due to the cancellation of gauge-dependent terms associated with the scalar component of the longitudinal polarisation of the W boson. The corresponding Goldstone boson couples to the top quark through the Yukawa coupling, which corresponds with the on-shell mass. On the other hand, the kinematic mass m may be off-shell and dictates the running of the width. We point out that here, and in the electroweak antenna functions, one could in principle account for the running of the Yukawa couplings, but such effects are currently neglected. Technically, sampling of the Breit-Wigner distribution with running width is achieved through rejection sampling using an overestimate distribution of the form where the second term is required to match the behaviour of the running width at high masses. The values of parameters b 1 through b 4 are optimised using a simple Monte Carlo procedure. If the resonance survives the EW shower long enough to reach its offshellness scale, it is decayed by selecting a channel with relative probability Due to its definite helicity, the angular distribution of the decay products is not uniform. A polar angle in the centre-of-mass frame is thus selected with probability proportional to the 1 → 2 matrix element summed over the final-state spins. Then, a helicity configuration is selected according to the individual spin channels, after which the decayed state may be constructed, boosted to the event frame, and inserted in the event.