Helicity Antenna Showers for Hadron Colliders

We present a complete set of helicity-dependent 2->3 antenna functions for QCD initial- and final- state radiation. The functions are implemented in the Vincia shower Monte Carlo framework and are used to generate showers for hadron-collider processes in which helicities are explicitly sampled (and conserved) at each step of the evolution. Although not capturing the full effects of spin correlations, the explicit helicity sampling does permit a significantly faster evaluation of fixed-order matrix-element corrections. A further speed increase is achieved via the implementation of a new fast library of analytical MHV amplitudes, while matrix elements from Madgraph are used for non-MHV configurations. A few examples of applications to QCD 2->2 processes are given, comparing the newly released Vincia 2.200 to Pythia 8.226.


Introduction
The description of bremsstrahlung processes in parton-shower event generators typically starts from the probability density for unpolarised partons to emit unpolarised radiation, i.e., DGLAP kernels or dipole/antenna functions summed over outgoing and averaged over incoming polarisations/helicities. One way of incorporating nontrivial polarisation effects, used in Pythia [1], is to correlate the plane in which a gluon is produced, with the plane in which it subsequently branches, taking linear-polarisation effects into account on the intermediate propagator, and casting the result in terms of a non-uniform selection of the azimuthal ϕ angle around the direction of the branching gluon, see, e.g., [2]. A more complete, but also more cumbersome, alternative, used in Herwig [3], is to keep track of spin correlations explicitly, using a spin-density matrix formalism [4][5][6]. In both cases, the nontrivial angular correlations ultimately arise from dot products between reference vectors expressing linear polarisations.
By contrast, a helicity basis does not rely on any external reference vectors, and hence helicity-dependence in and of itself does not generate any nontrivial angular correlations. Nonetheless, helicity-dependent radiation functions, as used for final-state radiation in Vincia for a few years [7], do have some advantages: helicity conservation can be made explicit, allowing to trace helicities through the shower; unphysical helicity configurations are prevented from contributing to sums and averages; and the explicit helicity assignments allow faster evaluations of matrix-element correction (MEC) factors, since only a single (or a few) helicity amplitudes need to be evaluated for each ME-corrected parton state [7].
The concept of ME corrections was first developed to improve the description of radiation in Pythia (then called Jetset) outside the collinear region to agree with first-order matrix elements for e + e − → 3 jets [8,9], and was since extended to correct the first emission in a wide range of resonance-decay processes and some (colour-singlet) production processes [10,11]. It was also used as a component of the first ME correction strategies in Herwig [12,13], and it forms the basis of the treatment of real corrections within the Powheg

Helicity-Dependent Showers
A helicity-dependent antenna shower for final-state radiation has already been introduced in [7]. The extension to hadronic initial states is straightforward. We start with a brief review of how emissions are generated and helicities selected. In cases where an event with unpolarised partons is showered by Vincia, a polariser function is first called, which uses helicity matrix elements to assign explicit helicities to all partons. Since the events are also assigned colour flows, we first define the joint probability to select a parton configuration with a colour flow i and a set of helicities h, where the full-colour (FC) and leading-colour (LC) matrix elements squared are defined by with M i the amplitude for colour-ordering i. We also make use of the notation for the fraction of the full-colour helicity matrix element squared that is projected onto LC colour flow i. As written here, the easiest would be to start by generating a helicity configuration, using the first factor in eq. (2) and then subsequently generate a colour flow using the second factor. For events which already have colour-flow assignments, the conditional probability for choosing helicity configuration h is simplest to define in terms of the redefined LC matrix elements, (If the corresponding matrix elements do not exist in Vincia, the event will remain unpolarised and showered using helicity-averaged and -summed antenna functions.) For events with explicit helicities, trial branchings are generated just as in the helicity-independent shower, i.e., using unpolarised trial-antenna function overestimates. After generating the post-branching kinematics (see, e.g., [17,18]), the total probability for accepting a branching (denoting pre-branching partons by AB and post-branching ones by ijk) 2 is: for fixed helicities h A,B of the parent partons. The sum over daughter helicities, h i,j,k , in the physical antenna function, A phys , runs over all possible (physical) helicities for the ijk partons, with each term, , being a helicity-dependent antenna function. To avoid clutter, and for ease of reference, we collect the precise forms for these functions in the appendix. We note that some of the functions differ (by nonsingular terms) from those used in previous versions of Vincia, in particular those in [7,18]. We also note that the accept probability defined by eq. (7) is in general identical to the unpolarised one (i.e., where one averages over h A and h B as well), up to nonsingular terms. In case of initial-state radiation, eq. (7) will be multiplied with the accept probability for the PDF ratios, just as in the unpolarised case [18]. Explicit helicities are then selected for the daughters according to the relative probabilities given by the antenna functions, where the denominator is equal to the numerator in eq. (7). Helicities are assigned to initial-state partons as well, using the same formalism. With the assumption that positive-helicity partons appear equally often as negative-helicity ones in the (anti)proton, the algorithm does not require any modifications when considering initial-state partons. Helicity conservation implies that, for gluon emission off (massless) quarks or final-state gluons, the parent partons do not change their helicities. A subtlety arises, however, for emissions off initial-state gluons. In the perspective of forwards evolution, such a branching looks like g I i → g I A g F j , where superscript I (F ) denotes an initial-state (final-state) parton; clearly, the helicity of parton i can be inherited by either parton j or parton A without violating helicity conservation. Hence the reader should not be confused by the appearance of physical initial-state antenna functions for which h A = h i in apps. A.3 and A.4, with corresponding DGLAP limits given in app. A.6.
LO antenna functions such as the ones discussed here are of course only accurate in the single-unresolved soft and collinear limits. To estimate the amount of uncertainty caused by physical shower emissions being away from these limits, we use a two-pronged approach based on reweighting [16]: 1) variation of the nonsingular terms of the antenna functions to estimate how close a given branching is to the logarithmicallydominated region, and 2) variations of the antenna-function renormalisation scale to estimate the potential impact of subleading-logarithmic terms. We emphasise that both types of variations are performed so that they preserve the total cross section (i.e., the variations appear with equal and opposite signs in real and virtual corrections, respectively [22]). The technical implementation in Vincia is quite similar to that in Pythia 8; see the respective HTML User Manuals and app. B. The variation of the renormalisation scale in a helicity-dependent shower is performed just as for an unpolarised shower, Optionally, an NLO-level compensating term can be introduced for gluon emission, which forces the variation to be equal to the result for k = 1 through order α 2 s : where b 0 = (33 − 2n F )/6π, n F is the number of active flavours at the scale t and m ant the mass of the parent antenna. The prefactor z is s ik /m 2 ant for final-final and m 2 ant /max(s ik , m 2 ant + s jk ) for initial-initial and initial-final branchings, with the post-branching invariants s ij , s jk , and s ik , The variation of the antenna functions by nonsingular terms, is performed such that the additional nonsingular term C NS /m 2 ant is distributed evenly amongst all helicity configurations for a specific antenna function, i.e. all helicity-dependent antenna functions obtain the same fraction of the nonsingular term. Note also that the nonsingular-term variations are cancelled by ME corrections (up to the corrected order) and therefore only need to be carried out for uncorrected orders.
For any given (bin of a) physical observable, a large dependence on C NS indicates that corrections from hard matrix elements with higher numbers of legs are needed, while a significant dependence on the renormalisation scale indicates a need for further corrections at the loop level.
Finally, it is worth emphasising that the statistical fluctuations of the uncertainty variations are generally larger than for the central (non-varied) predictions. This is due to the central prediction being unweighted (in our setup) and the the variations being computed by reweighting. See [23] for an example of how weighting ("biasing") the central distribution can improve the relative statistical precision of the uncertainty bands.

Matrix-Element Corrections
The GKS formalism for iterated matrix-element corrections [16] was originally based on so-called smoothly ordered showers, with a Markovian (history-independent) choice of restart scale after each branching. This allows the shower algorithm to generate phase-space points that violate the nominal ordering condition of the shower, at a suppressed but still non-zero rate, thus filling previously inaccessible regions of phase space; the correct (tree-level) emission rates can then be obtained via matrix-element corrections just as in the ordered part of phase space. However, general arguments indicate that the effective Sudakov factors for the nonordered histories, are probably not correct [18,24,25]. Recent efforts [19,25] have therefore shifted focus back to filling the phase space for multiple hard emissions while remaining within the paradigm of strong ordering. In particular, we take the strongly-ordered iterated-MEC formalism presented in [19] as our starting point, and adapt it to include explicit helicities.
The question of Markovian vs non-Markovian behaviour comes about since the value of the shower evolution parameter in conventional strongly-ordered showers depends on which parton was the last one to be emitted. This cannot be uniquely determined merely by considering a given parton configuration; the value is a function of what shower history (or path) led to the configuration in question; a non-Markovian aspect. In the context of iterated ME corrections, non-Markovianity implies that the MEC factors contain nested sums over shower histories involving clusterings all the way back to the Born configuration (while a Markovian algorithm only requires a single level of clusterings [16]).
Within the formalism presented in [19], the splitting kernels are redefined by multiplying them with the correction factor |M(Φ n+1 )| 2 denotes the matrix element squared of the Φ n+1 state and A (Φ n+1 /Φ n ) the antenna function, associated with the clustering Φ n+1 → Φ n . The denominator sums over all possible ways the shower could have produced the n+1-particle state Φ n+1 from a given Born state Φ 0 , including the correct weights of every shower step on the way. This yields the recursive structure of eq. (12) and the dependence on the correction factors of the previous orders. In addition the (process-dependent) scale t(Φ 0 ), at which the shower starts the evolution off the Born state is taken into account.
For a helicity-dependent correction, we modify eq. (12) such that, for a given polarised Φ n state, the sums over the intermediate states Φ n−1 . . . Φ 0 are extended to include all possible helicity configurations. As an example, consider a possible clustering of a final-state qq pair to a gluon. In the unpolarised case, one term corresponding to the clustering qq → g contributes with the respective unpolarised antenna function and matrix element (which both implicitly involve helicity sums of course). For a polarised q +q− pair, two different clustered helicity states are possible, q +q− → g + and q +q− → g − , each contributing according to their antenna function and matrix element. The evolution variable, however, is the same as in the unpolarised case. This concludes our discussion of helicity-dependent matrix element corrections.

MHV Basics
For fast evaluation of certain types of helicity configurations Vincia uses Maximally Helicity Violating (MHV) amplitudes. MHV amplitudes have the advantage of being compact analytical expressions which are independent of Feynman diagrams; see [26,27] for reviews. In this section, we briefly introduce the concepts and notation relevant to understanding the conventions and properties of the small library of MHV amplitudes implemented in Vincia.
In the following we consider all particles to be outgoing and massless. We recall that in this limit a particle's helicity corresponds to its chirality, and define our spinors in the helicity basis: The notation ij and [ij] is used for inner products of such spinors: in terms of the (light-cone) momentum For more details about spinor inner products and their properties see e.g. [26,27]. Note that in recent literature one often finds the convention [ij] = ij * , which is different to above (see e.g. [28]).
In the all-outgoing convention, helicity conservation implies that at least two pairs of opposite-helicity partons must exist for an n-parton amplitude to be nonzero 3 . If the remaining n − 4 partons are all chosen to be of the same helicity (+ or −), the amplitude is called maximally helicity violating (MHV), and has a remarkably simple structure. The first MHV amplitude to be discovered was the all-gluon Parke-Taylor amplitude [29]. In the following years this was extended to include one [30,31] and two [32][33][34] quark pairs, as well as to the case of a quark pair and a massive vector boson which decays leptonically [35,36].
All-Gluon Amplitudes: To use these amplitudes we first note that the colour information can be factorised from the kinematics. In the n-point all-gluon case we use: where g s is the strong coupling (g 2 s = 4πα s ), the normalisation convention is t a = λ a √ 2 with λ a being the generators of SU (3), p i is the gluon momentum, h i the gluon helicity, Tr(t a σ(1) . . . t a σ(1) ) the colour factor and A n (σ(p h1 1 ), . . . , σ(p hn n )) the kinematic part of the amplitude. The sum is over all non-cyclic permutations σ of the particles. The Parke-Taylor amplitude then describes the kinematic part of eq. (16) and is given by: where gluons i and j have negative helicity, and all other particles have positive helicity.
One Quark Pair: If we add a qq pair we require that the quark and antiquark have opposite helicities (consistent with the gluon having spin 1), and use the following colour basis: where q, h q , and i (q, hq, and j) are respectively the quark (anti-quark) momentum, helicity, and colour index; and the sum is over all permutations of the gluons. If the quark and gluon i each have negative helicity and all other particles positive helicity, then the kinematic amplitude is the given by: where the numbers refer to the (colour-ordered) gluons. If we exchange the helicities on the quarks, it is sufficient to exchange the exponents in the numerator: Two Quark Pairs: The four-quark, n − 4 gluon colour structure is given by: where {ij} = ij for positive-helicity gluons and {ij} = [ji] for negative-helicity gluons; q and Q label the two quark lines; A 0 (h q , h Q , h g ) is a kinematic function which depends on the helicities of the two quarks and the gluons, with opposite-helicity cases obtained using parity transformation ij ↔ [ji]; and the two functions A (0) n and A (1) n are kinematic amplitudes, for which we have used the short-hand notation q ≡ q hq , i ≡ σ(p hi i ) etc.: We must sum over all possible partitions of gluons between the two quark colour lines, and also over all possible permutations of gluons within those partitions. If there are no gluons propagating off a particular colour line, then that colour line is described by a Kronecker delta. Note that this decomposition only works for the MHV configuration.
Drell-Yan, DIS, and hadronic Z decays: To create MHV amplitudes with a single quark pair, a single lepton pair, and an arbitrary number of gluons, the four-quark amplitude can be recycled with all gluons coming from a single quark line. The second quark line is now equivalent to a ll pair up to couplings and an overall propagator factor. The amplitude then has the form where the sum is again over all gluon permutations. The kinematic amplitude is given by where the braces have the same meaning as in eq. (24), and the function M l V is given by is the coupling of lepton l (quark q) with helicity h l (h q ) to vector V , and M V and Γ V are the mass and width of the vector boson respectively.
Finally, we remark that in all of the above expressions, flipping the helicity of every particle is equivalent to exchanging each ij ↔ [ji]. This concludes our brief recapitulation of the basics of the MHV formalism and convention choices.

MHV within Vincia
The MHV amplitudes that are made available in standalone Vincia are summarised in tab. 1. Note that these amplitudes are so far only used for QCD 2 → n matrix-element corrections, and that the second quark pair must have a different flavour to the first.
The colour-summed squared matrix element is calculated using the following matrix equation: where FC stands for the full colour-summed matrix element squared as in eq. (3), C ij is a colour matrix obtained by multiplying the colour factor from permutation σ i onto the conjugate colour factor from σ j , and the sum is over all colour orders. We optimise the all-gluon amplitudes by diagonalising C ij for the 4 and 5-gluon matrix elements, and partially diagonalising C ij for the 6-gluon matrix element as done in [26]. By default, Vincia uses MHV amplitudes wherever possible to compute its matrix-element correction factors, thus ensuring the fastest possible run time. However, this can be turned off (e.g., for cross checks with amplitudes from Madgraph) using the flag vincia:useMHVamplitudes. To calculate an MHV ME correction, Vincia actively crosses the initial-state partons into the final state, rearranges the partons into the correct colour order, calculates all of the explicit spinor products needed, and then calculates the matrix element squared.
The calculation of ME corrections for MHV configurations exhibits the nice feature that all clustered states in eq. (12) are MHV configurations as well. Helicity conservation does not allow ++ → − nor −− → + clusterings (in the all-outgoing convention). This results in clustered states being either MHV configurations themselves or unphysical states with a vanishing amplitude. Consider n positive-and 2 negative-helicity outgoing partons as an example. Here clustered states contain either n − 1 positive-and 2 negative-helicity partons (MHV) or n positive-and 1 negative-helicity partons (unphysical).
For instructions on how to use Vincia for calculating spinor products or MHV amplitude in a standalone context, see the online user guide [37].

Polarising events with MHV
The fact that Vincia assigns helicities to unpolarised events, with relative probabilities according to the corresponding helicity matrix elements squared, was briefly discussed in sec. 2. An interesting simplification occurs when all of the contributing amplitudes are of the MHV kind, as is, e.g., the case for all QCD 2 → 2 and 2 → 3 processes. The simplification follows by noting that the full-colour (FC) MHV matrix elements squared all have the following form (so long as there is at most one quark pair): where h is a label denoting the helicity assignments, M h n ≡ |A h n (1, . . . , n)| 2 is some function of the helicities and momenta, σ is the relevant set of permutations, CF is the relevant colour factor at the amplitude level, and | σ F (σ)| 2 is the square of the sum over colour permutations. For example, in the all-gluon amplitude A h n (1, . . . , n) could be ij 4 . We have therefore factored out the helicity information M h n from the colour information. This also works for the LC matrix elements LC h i which are given by eq. (29) above without the sum of permutations. That is, Recall that the conditional probability defined in (6) used to pick helicities for configurations that already have colour assignments has the form: We can use eq. (29) to simplify this: while the mostly plus factors are given by the usual parity relation. Note that this also holds for the full-colour amplitudes used for selecting helicities at the colour-summed level, cf. eq. (2), The preceding argument also works for 4-quark MHV amplitudes with distinct quark pairs provided one changes eq. (29) to include the second colour connection. However, this doesn't work for all 4-quark MHV amplitudes because there is an extra colour-connection when two identical quarks have the same helicity. Hence the colour factor depends on the helicity and cannot be factorised.

Speed Comparisons
At the level of a pure shower (before ME corrections are imposed), the change from helicity-summed to helicity-sampling radiation functions requires the generation of one more random number per n → n + 1 branching, to select the helicity of the emitted parton. This comes in addition to at least three random numbers for the one-particle phase space. All else being equal, a helicity-sampling shower should therefore not be more than a factor 4/3 slower than a helicity-summed one. (Similar arguments hold for the initial polarisation step for the hard process). However, since there are many common components which must be computed regardless of the choice of helicity treatment, one expects the effective slowdown of the full shower algorithm to be milder than this upper limit. This is also borne out by explicit tests with Vincia, which exhibit slowdowns of less than 10% when switching on helicity-sampling. (See also the first bin of fig. 1 below.) As a measure of the relative speed of helicity-dependent vs helicity-summed ME corrections, and the difference between using MHV matrix elements or Madgraph 4 ones, we consider the following specific (but fairly representative) benchmark case: qg → qg Born-level processes, with a minimump ⊥ of 100 GeV, in pp collisions with E cm = 10 TeV. A technical point is that, for this comparison, we switch g → qq branchings off  Figure 1: Speed comparison for helicity-independent ("Non-hel") and helicity-dependent ("Hel") showers as a function of the number of ME-corrected legs, for qg → qg + gluons withp ⊥min = 100 GeV, for pp collisions at E cm = 10 TeV. The dashed horisontal line indicates the time it takes to generate MPI and hadronisation for the same events. Results were obtained from 10,000 events generated for each run, on a single 2.9 GHz Intel Core i7 processor, using the clang compiler (v3.9), with -O2 optimisation.
in the shower, so that the generated shower configurations are all of the simple qg → qg + gluons type. This allows us to illustrate speeds of ME corrections with up to three additional legs while, if g → qq branchings had been switched on, the current version of Vincia is restricted to ME corrections with up to two additional legs. (This restriction will be lifted in a future update.) Fig. 1 illustrates the number of milliseconds it takes to generate one shower, as a function of the number of legs that are requested to be ME-corrected. The solid (red) line without symbols uses helicity-summed showers and matrix elements, while the two blue curves (with symbols) show the dependence of the helicitydependent formalism, with or without enabling the library of MHV matrix elements, respectively. For reference, the thick dashed horisontal line shows the time it takes to generate multi-parton interactions (MPI) and hadronisation for the same events 4 . For 0 or 1 corrected emissions, the helicity-summed shower is actually slightly faster, since the Born-level polariser and the helicity selection in the shower take a little extra time and the first-order ME corrections are very quick to evaluate even when summing over helicities. At two legs, however, the helicity-dependent formalism is up to 30% quicker (with the MHV library switched on) than the helicity-summed one. At three legs, the difference is a factor 4, with the MHV library allowing to shave an extra ∼ 15% off the shower-generation time relative to using only MG4 matrix elements.
One also notices that by two corrected legs, the showering time is becoming comparable to the time it takes to generate MPI and hadronisation for the events, hence this is the point at which the showering speed would start to be felt in the context of generating full events. By three corrected legs, the ME corrections dominate the event-generation time. The default in the current version of Vincia is that ME corrections are enabled for QCD 2 → 2 processes up to two additional legs; the event-generation time should therefore stay within roughly a factor 2 of that of the uncorrected algorithm. The complete set of matrix elements required for 3rd-order corrections will be provided in a future update. For hadronic Z, W , and H production or decay, the full set of 3rd-order matrix elements are already available in the current version. (We note that the implementation of the iterated-MEC algorithm itself is general and could in principle handle any number of legs, if provided with the required matrix elements.)

Example Application
To illustrate the properties of the ME-corrected algorithm (and uncertainty variations) in the context of a realistic application, we consider showers off gg → gg Born-level events and compare Pythia 8.226 and Vincia 2.200 on three observables sensitive to different aspects of the evolution: early branchings, late branchings, and a polarisation effect, respectively: 1. Early branchings: the 3-jet resolution scale, d 23 , using the longitudinally invariant k ⊥ -jet algorithm with R = 0.4. The analysis is adapted from the code used in [39], originally written by S. Höche in the Rivet [?] analysis framework.
2. Late branchings: the 6-jet k ⊥ resolution scale, d 56 , with the same jet algorithm and analysis as above.
3. Gluon polarisation: the angle between the event plane (characteristic of the original gg → gg Bornlevel event) and the plane of a subsequent g → bb splitting. Here, the anti-k ⊥ jet algorithm with R = 0.2 is used (so that the b jets can be resolved down to small separations), and we impose a minimum jet p ⊥ of 50 GeV. The analysis is performed in the Rivet + Fastjet [40] framework. (For further ideas on how to exploit heavy-flavour tags to probe g → qq splittings at colliders, see e.g. [42,43].) The basic 2 → 2 QCD process is sampled with the cutp ⊥ ≥ 500 GeV on the final-state partons. For consistency with the shower α s parameters, Vincia's default tune uses two-loop running for the strong coupling with α s (m 2 Z ) = 0.118 for the hard process. To compare predictions on an equal footing we apply the same settings for the underlying Born process in Pythia. To focus on the showering off the hard process all comparisons are done with multiparton interactions switched off.
To obtain dimensionless variables, the jet resolution measures d 23 and d 56 are normalised by a factor 1/d 12 , i.e., they are effectively measured relative to a scale representing thep 2 ⊥ scale of the underlying Born process 5 . The resulting quantities exhibit a fixed-order behaviour for large values and a Sudakov suppression for low values. Especially for well-resolved radiation, we therefore expect these observables to be sensitive to low-order ME corrections, and hence the uncertainty associated with nonsingular-term variations should be reduced when Vincia's ME corrections are switched on. (Note: Pythia does not incorporate ME corrections for QCD 2 → 2 processes.) Parton-level results for showered gg → gg events are presented in fig. 2 with uncertainty bands.
The ME corrections in strongly-ordered events exhibit a modest effect of up to 20% for large values of d 23 /d 12 and d 56 /d 12 , with the ME-corrected rate being larger than that of the pure Vincia shower. Shape differences between the predictions of Pythia and Vincia are visible throughout most of the distributions, with the uncorrected Vincia shower generating a somewhat harder d 23 /d 12 spectrum than Pythia. ME corrections increase the rate for large d 56 /d 12 values, bringing the predictions of Vincia closer to that of Pythia. Given the different choices of shower α s parameters, evolution variable, and radiation functions, we do not consider this level of disagreement between the two models surprising. The evolution of the hard process starts at the factorisation scale for both showers. However, depending on the form of evolution variable, the hardest possible scales correspond to different values of d 23 .
All predictions exhibit some rather large fluctuations in the uncertainty bands. The dijet system with the cutp ⊥ ≥ 500 GeV as underlying hard process is typically accompanied by a large number of additional jets. Given the nature of the reweighting algorithm of [22] (and similarly for [23,44]) this may easily result in fluctuating weights. In addition we expect larger fluctuations in the nonsingular-term variations for the helicity shower, compared to the helicity-independent one. As discussed in sec. 2, the additional nonsingular terms are distributed evenly between all helicity configurations. This results in a larger spread of weights, when considering helicity configurations that constitute either a large or a small fraction of the helicity-summed antenna functions. To mitigate the effects of weight fluctuations, we conclude that further development of these reweighting methods would be useful, in particular for large phase spaces (long shower chains). E.g., the authors in [23] have demonstrated that combining biasing with reweighting can improve the relative  statistical precision of the uncertainty variations, at the price of generating some reasonably well-behaved weights for the central (non-varied) event sample.
The variation of the nonsingular terms (hashed bands) results in a larger band around small |d 23 /d 12 | and |d 56 /d 12 | for Vincia without ME corrections, compared to Pythia. The ME corrections cancel the effect of varying the nonsingular terms in the radiation functions. Consequently, the respective uncertainty band for Vincia with ME corrections is very narrow, especially for d 23 . The renormalisation-scale variations (shaded bands) are quite similar in size for all predictions. They show the largest effect for small jet separation scales, where soft emissions and the Sudakov factor contribute to the distribution.
We now turn to an observable where polarisation effects are expected to contribute. In events with two b-jets a plane is defined by the two jets. A second plane is defined by the gluon-jet (the sum of the two b-jets) and the beam axis. In fig. 3 the angle between the two planes is shown. A flat distribution is obtained with Pythia without gluon polarisation effects in the final-state shower and Vincia without ME corrections. However, Vincia produces an around 15% higher total rate, compared to Pythia. We note that both codes generate a similar total rate of g → bb splittings in the shower, where the gluon splittings occur "later" in the evolution in Pythia (i.e., preceded by a larger number of other branchings). The b-quarks are therefore more likely to obtain a smaller invariant mass and might be clustered within the same jet. Together with the p ⊥ and invariant mass cuts on the jets, this may cause a smaller rate of events with two b-jets. The polarisation effects in Pythia leave the total rate unchanged, but increase the amount of events where the angle is close to π/2. The ME corrections in Vincia change the total rate by decreasing the number of events with splitting angles near 90 degrees. The qualitative effect is therefore the opposite of that in Pythia, where the total shower rate is preserved, but the region around 90 degrees is enhanced by the polarisation effect. We conclude that a measurement of this observable, and the development of alternative strategies for corrections beyond fixed order (e.g., along the lines proposed in [25]), would be desirable.

Conclusions
We have presented a helicity-dependent antenna shower for QCD initial-and final-state radiation, implemented in the Vincia shower model. The iterated ME correction formalism of [7,16,18,19] has been extended to cope with helicity-dependent clusterings and splitting kernels involving initial-state legs, and in this work has been applied to strongly ordered showers in a direct extension of the formalism presented in [19]. We further reported on new, user-specifiable uncertainty variations in Vincia, including renormalisation-scale and splitting-kernel variations.
The new approach and a library for tree-level MHV amplitudes enable a faster evaluation of MEC factors, as illustrated explicitly for the process qg → qg + gluons. While the pure shower is slightly slower due to the additional step of helicity selection, the evaluation of ME corrections can be done significantly faster when only a single or a few helicity matrix elements need to be evaluated per trial branching, relative to when helicity-summed matrix elements are used.
To illustrate the effect of the iterated ME corrections and uncertainty variations within the helicitydependent shower, we considered a few representative observables, based on showered gg → gg Born-level events. As expected, ME corrections reduce the overall amount of variation considerably in regions of relatively hard emissions, where process-dependent nonsingular terms (captured by the matrix elements) dominate over the universal logarithmic terms (captured by the showers). In regions of large scale hierarchies, the uncertainty due to renormalisation-scale variations dominates and remains uncompensated by tree-level ME corrections.
We also showed a more complex example, the angle between a Born-level gg → gg event plane and the plane of a subsequent g → bb splitting. In Pythia, a general implementation of gluon polarisation effects implies an enhancement of such splittings at 90 degrees to the original event plane (while the total shower rate of g → bb splittings is preserved); while in Vincia, ME corrections dominantly act to suppress the overall rate of g → bb splittings. Moreover, the suppression is most active for the most well-resolved branchings (at 90 degrees), leading to an opposite-sign effect than the one in Pythia. We conclude that there is a complex interplay between the rate and the angular dependence of these branchings, and intend to investigate this further in future studies.
Acknowledgments: AL and PS acknowledge support from the Monash-Warwick Alliance Development Fund Project "Collider Physics". PS is the recipient of an Australian Research Council Future Fellowship, FT130100744.

A.1 Notation and Conventions
We use capital letters to denote partons in the pre-branching n-parton configuration and lower-case letters to denote partons in the post-branching (n + 1)−parton configuration. Incoming partons are denoted a, b, while final-state partons are denoted i, j, k. Thus, for example, an initial-final antenna branching is written AK → ajk.
The scaled branching invariants for final-final antenna functions are and the energy fractions Note that, for gluon-emission antennae involving massive parent quarks, a helicity-independent negative correction to the eikonal is added, with helicity-summed average: For gluon-splitting antennae (Xg → Xq j q k ), the mass correction is positive:

A.2 QQ parents: Gluon Emission
The helicity averages for qq → qgq antennae are II : a(q A q B →q a g j q b ) = 1 s AB 2y AB y aj y jb + y jb y aj + y aj y jb IF : a(q A q K → q a g j q k ) = 1 s AK (1 − y aj ) 2 + (1 − y jk ) 2 y aj y jk where the slightly different nonsingular terms chosen for the IF case ensure positivity of in particular the (++ → + − +) antenna function over all of the IF phase space, while the nonsingular terms for the FF and II cases result from averaging over the corresponding helicity matrix elements for Z and H decays. The individual helicity contributions are: FF : II : IF :
The helicity averages for qX → gqX antennae (quark backwards evolving to a gluon) are IF : a(q A X K → g aqj X k ) = 1 s AK The individual helicity contributions are: IF : a(−X → − + X) = 1 s AK a(−X → + + X) = 1 s AK (1 − y AK ) 2 y aj .
The helicity averages for gX → qqX antennae (gluon backwards evolving to a quark) are II : a(g A X B → q a q j X b ) = 1 s AB 1 + (1 − y AB ) 2 2y aj (1 − y jb ) , IF : a(g A X K → q a q j X k ) = 1 s AK 1 + (1 − y AK ) 2 2y aj (y AK + y aj ) .

B Details of VINCIA Implementation
For completeness, we also report on the following changes in Vincia with respect to [18]: • The so-called "Ariadne factor" [45] for gluon splitting antennae has been removed completely, as it has only been applied to 4-jet events in hadronic Z decays and its influence cancels once ME corrections are used in the evolution.
• The CMW-rescaling of α s [46] is no longer applied to the soft-eikonal terms of the antenna functions, but rather as a global rescaling of Λ QCD , independent of the type of branching.
• By default the power shower approach [10,47] is used for hard process without QCD partons in the final state. This obviates the need for a separate event sample containing jets associated with scales larger than the factorisation scale, which has been introduced in [18]. For QCD-type processes the shower starts the evolution at the factorisation scale.
• The so-called "smooth ordering" [16], which allows the shower to populate phase-space regions beyond the reach of traditional ordered showers, is no longer used. Consequently, the MECs formalism so far used in Vincia is no longer applicable and the MECs method for ordered showers of [19] is applied. See sec. 3 for a brief review of the formalism.
• The CKKW-L merging implementation in Pythia 8 [48] is now also available in Vincia, making use of the parameters in Pythia 8. This allows to supplement the MECs method for ordered showers with non-shower-like events, as discussed in [19]. Note however, that it is not possible to combine the merging procedure with the helicity-dependent shower.
• Automated uncertainty variations are now user specifiable in the same manner as in Pythia 8 [22], with the following keywords controlling the type and size of variations in Vincia, for final-final (FF), initial-final (IF), and initial-initial (II) antennae respectively: • Renormalisation-scale variations (applied to all antenna functions): ff:muRfac ; if:muRfac ; ii:muRfac .
• Optionally, antenna-specific variations can be specified, which then supersede the antenna-independent variations. The full set of possible keywords is listed in Vincia's HTML User Reference.