Spin Correlations in Parton Shower Simulations

Spin correlations are an important, but often neglected, effect in modern Monte Carlo event generators. We show that they can be fully incorporated in Herwig 7 using the algorithm originally proposed by Collins and Knowles in all stages of the event generation process and between the different stages of the event generation. In this paper we present the final missing ingredient, correlations in both the angular-ordered and dipole shower algorithms and between the parton shower and hard production and decay processes.


Introduction
Monte Carlo event generators are essential in modern particle physics as they provide state-of-the-art theoretical calculations which can be directly compared to experimental results. 1 In recent years most of the developments in these programs have focused on improving the accuracy of the hard perturbative calculation by matching with higher order and higher multiplicity matrix elements. There has however been less work focused on improving the accuracy of the parton shower algorithm itself. While there is some older work [2][3][4][5] there has recently been a revival of interest in including higher order corrections in the parton shower evolution in antenna [6][7][8] and dipole [9][10][11] showers as well as work on amplitude-based evolution to treat subleading colour effects [12,13].
One important, but often overlooked, effect is the azimuthal correlation of emissions which occurs even at leading order in the parton shower evolution due to the polarization of the gluons. Spin correlations also give rise to correlations between the emissions in the parton shower and the hard process, and between the production and decay of heavy particles, such as the top quark. In order to fully include 1 See Ref. [1] for a recent review. a e-mail: Peter.Richardson@durham.ac.uk these correlations, such that the complexity of the calculation only grows linearly with the number of final-state particles, the algorithm of Refs. [14][15][16][17][18] can be used. This algorithm was used in Herwig6 [19] to implement the correlations in the parton shower, between the hard process and the parton shower emissions [15][16][17] and between the production and decay of heavy particles, in both the Standard Model, the MSSM and its R-parity violating extension. 2 The same algorithm is also used in the EvtGen package [20] for correlations in the decays of hadrons.
In Herwig++ [21] and Herwig7 [22,23] this algorithm was used to provide correlations between the production and decay of heavy particles, and in the decay of hadronic resonances, in a unified framework. However the correlations in the parton shower were not included. In this paper we will present the necessary calculations in order to implement these correlations in both parton shower algorithms available in Herwig7 and results from their implementation which will be available in Herwig7.2.
In the next section we will first recap the details of the algorithm of Refs. [14][15][16][17][18] and then present details of the calculations needed for its implementation in the Herwig7 parton showers. We then present some results from the implementation in Herwig7 followed by our conclusions.

Review of the algorithm
It is best to consider the details of the algorithm of Refs. [14][15][16][17][18] by considering an example. In order to consider the correlations in the parton shower we will consider the example of the decay of the Higgs boson to two gluons, followed by the branching of the gluons into quark-antiquark pairs, i.e. h 0 → gg → qqq q .
The algorithm proceeds by first generating the momenta of the gluons produced in the Higgs boson decay using the appropriate matrix element. Obviously in this case this is trivial, however the structure of the matrix element will give rise to the subsequent correlations. One of the outgoing gluons is then selected at random and a spin density matrix calculated for its subsequent evolution where M h 0 →gg is the helicity amplitude for the decay of the Higgs boson into two gluons, λ g 1 and λ g 1 are the polarizations of the first gluon and λ g 2 is the polarization of the second gluon. The Einstein convention where repeated indices are summed over is used. The normalization N is defined such that the trace of the spin density matrix Tr ρ g 1 = 1 and is therefore equal to the spin-averaged matrix element squared used to calculate the momenta of the gluons in the previous stage of the algorithm.
We can now proceed and generate the parton shower emissions from the gluon in the standard way because the probability of such an emission occurring is not affected by the spin density matrix. 3 However the azimuthal angle of the emission is affected by the spin density matrix and is generated using the distribution where M g→qq is the helicity amplitude for the splitting of a gluon into a quark-antiquark pair and λ q and λq are the helicities of the quark and antiquark, respectively. Here we have used the spin-density matrix to encode the correlations from the hard process, rather than just averaging over the polarizations of the branching gluon. This relies on the standard factorization of the matrix element in the collinear limit into the production of an on-shell gluon and its subsequent branching, which forms the basis of the parton shower approach. Similarly when generating the decay of a particle the amplitude factorizes in the narrow width limit into an amplitude for the production of the particle and one for its decay.
In the case that we have further emissions from the quarks produced in the gluon branching we would then calculate a spin density matrix for emission from the quarks in the same way as above, i.e. 3 While this is true for emissions in both QCD and QED if we were to consider the emission of electroweak W ± and Z 0 bosons then the emission probability will also depend on the spin density matrix as the diagonal elements give the probability of having a left-or right-handed particle radiating.
where again the normalization is such that Tr ρ q = 1 and is equal to the spin contracted amplitude used to calculate the distribution of the g → qq splitting, Eq. 2.
Once there is no further emission, or decay of the particle, we can compute a decay matrix for the final parton shower branching or decay. In our example if there is no further emission from the quark or antiquark then the decay matrix for the gluon is where again the normalization is such that Tr D g 1 = 1. If there had been emissions from either the quark or antiquark then instead of summing over their helicities we would contract with the appropriate decay matrix. These decay matrices encode the distributions generated in the branchings so that any subsequent emissions are correctly correlated with them. We now need to generate the branching of the second gluon from the hard process. We first calculate the spin density matrix, however now instead of summing over the helicities of the other gluon we contract with the decay matrix encoding how the first gluon branched, i.e.
where again the normalization is such that Trρ g 2 = 1.
As before this spin density matrix is then contracted with the matrix elements for the branching of the gluon to determine the azimuthal angle according to Eq. 2 with g 1 → g 2 .
It is worth noting that due to the choice of the normalization of the spin density matrices the normalization used in the next step of the calculation is always equal to the distribution used to generate the previous step. This ensures that the final distribution is equal to the full result, up to the approximation used to factorize the full matrix element into different components. In the example we have been considering the final result used to generate the momenta of the particles is i.e. the full matrix element for the process in the collinear limit.
The same arguments also apply for initial-state radiation with the rôles of the spin density and decay matrices exchanged, as we evolve backwards from the hard process [14][15][16][17][18]. A similar example, considering the production and decay of a top quark-antiquark pair can be found in Ref. [21].
It is worth noting that we could consider the spin sum used where a branching or decay has not been generated as contracting with a spin density or decay matrix proportional to the identity matrix. This can be generalised, i.e. if we have information on the production, decay or branching of a particle we should use the appropriate spin density or decay matrix. This allows us to change the order in which we calculate the decay or branchings if required. For example, we postpone the generation of any tau lepton decays until after the parton shower has been generated to avoid applying large boosts to the low invariant mass systems produced which leads to numerical issues with momentum conservation. This has the price of increasing the number of numerical evaluations as any spin density and decay matrices which are affected by the decay then need to be recalculated to ensure consistency, however this is still a linear operation. This is also required in the dipole shower where successive emissions can be from different particles.

Helicity amplitudes for parton branching
In order to implement the algorithm described above we need the helicity amplitudes for the branching processes in the same limit as that used in the parton shower. In modern parton shower algorithms this is the quasicollinear limit of Ref. [24]. In order to efficiently implement the algorithm we need compact analytic expressions for these branchings, with the correct approximations, rather than relying on numerical evaluations of the helicity amplitudes, which are used in Herwig7 for the amplitudes for production and decay processes.
We choose to calculate the branchings in a frame where the branching particle is along the z-axis, which is the choice made in the angular-ordered parton shower algorithm [25] in Herwig7. In this frame the momentum of the branching particle is where m 0 is the on-shell mass of the branching particle and p is the magnitude of its three-momentum. In the parametrization used in Herwig7 all momenta are decomposed in the Sudakov basis, i.e. the momentum of particle i can be written as where p is the momentum of the branching particle before any emission, n is a light-like reference vector, α i and β i are coefficients, and q ⊥i is the transverse momentum relative to the directions of p and n. The reference vectors and transverse momenta satisfy p 2 = m 2 0 , n 2 = 0, p·q ⊥i = n·q ⊥i = 0 and p · n > 0. As we need an explicit representation we take If we consider a branching 0 → 1, 2 then in this parameterization the momenta of the particles produced in the branching are: (10a) where z is the light-cone momentum fraction of the first parton and where m 1,2 are the on-shell masses of the particles produced, p ⊥ is the magnitude of the transverse momentum of the branching and φ is the azimuthal angle for the first particle produced in the branching.
In the Herwig7 angular-ordered parton shower we use the evolution variablẽ for final-state evolution. The momentum of the off-shell particle initiating the branching is where such that the virtuality of the branching parton is We need to evaluate the branchings in the quasi-collinear limit in which we take the masses and transverse momentum to zero while keeping the ratio of the masses to the transverse momentum fixed. Practically this is most easily achieved by rescaling the masses and transverse momentum by a parameter λ and expanding in λ [24], i.e. 123 In order to define our conventions for spinors and polarization vectors we present the calculation of the g → qq splitting here. The results for the remaining QCD branchings are presented in Appendix A.
In this case we start with an on-shell gluon such that after rotating so that the gluon is along the z-direction the momentum and polarization vectors of the gluon are: where λ 0 = ±1 is the helicity of the gluon. For the branching g → qq , m 1 = m 2 = m, therefore expanding in λ the momenta of the particles in the branching are: In general we define the spinors using the conventions in [26]. This allows us to write the spinors for the outgoing particles: We define the helicity amplitudes F λ 0 λ 1 λ 2 for the branching such that the spin averaged splitting function is For the g → qq branching There is no simple form for all the helicities and therefore the functions, F λ 0 λ 1 λ 2 , are given in Table 2. These amplitudes give the correct spin averaged massive splitting function and reduce to the splitting functions used in Herwig6 from [15] Table 1 in the massless limit. 4

Examples
It is instructive to calculate the correlations in some simple cases. This will both illustrate how the spin correlation algorithm works and provide analytical results that we later use to check the general implementation in Herwig7.
The simplest example is the correlation of the angle between the planes of two successive parton shower branchings, i.e. 0 → 12 followed by 2 → 34. The only non-zero correlation is due to the polarization of an intermediate gluon.
If we consider the branching q → qg the spin density matrix for the radiated gluon is where z 1 is the momentum fraction and φ 1 the azimuthal angle of the quark, respectively. We have neglected the mass of the radiating quark. Similarly for the branching g → gg the spin density matrix for the radiated gluon is If we contract these spin density matrices with the appropriate helicity amplitudes for the subsequent branching of the gluon then we obtain the distribution where Δφ = φ 2 − φ 1 is the angle between the planes of the two branchings and the coefficients A and B are given in Table 1.
The simplest example in which there are final-state correlations in the parton shower due to the production process is the decay of the Higgs boson to two gluons followed by the branching of each of the two gluons into a quark-antiquark pair, i.e. h 0 → gg → qqq q , which we have already considered.
We take the amplitude to be 5 where p i=1,2 and i=1,2 are the 4-momenta and polarization vectors of the outgoing gluons, respectively. We have neglected the overall normalization of the matrix element which affects the total rate but not the correlations we are studying.
The non-zero helicity amplitudes for h 0 → gg are: where φ is the azimuthal angle of the 1st gluon. We can use the helicity amplitudes in Table 2 to calculate the decay matrix for the first gluon branching where φ 1 is the azimuthal angle of the quark q in a frame in which the 1st gluon is along the z-axis and where z 1 is the energy fraction of the quark,q 1 the evolution variable for the branching and m q is the mass of the quark. We can then calculate the spin density matrix needed to calculate the second branching, This gives the distribution where z 2 is the energy fraction of the quark q andq 2 the evolution variable for the second branching. The azimuthal angle of the second branching, φ 2 , is measured in a frame where the 2nd gluon is along the z-direction. If we rotate the quark produced in the first branching into this frame its angle If we neglect the mass of the quark and multiply the distribution by the spin-averaged splitting function for the first branching we obtain the distribution Integrating over z 1,2 between 0 and 1 and normalizing the resulting distribution gives where Δφ = φ 2 − φ 1 is the azimuthal angle between the planes of the two branchings.

Implementation in the angular-ordered parton shower
The algorithm described in Sect. 2.1 is ideal for implementation in the angular-ordered parton shower in Herwig7. In the angular-ordered parton shower each parton from a hard process or decay is showered independently and the individual splittings are calculated in the same order as that described in Sect. 2.1.

(2020) 80:83
There is one complication. In the angular-ordered shower, each parton incoming to or outgoing from the hard process is showered in a specific frame. The transformation between the production frame and showering frame of a parton can introduce a rotation of the basis states of the parton which must be taken into account in the treatment of the spin density and decay matrices of the parton. This is done through the application of 'mappings' to the matrices and is described in detail in Appendix B.

Implementation in the dipole shower
Several additional considerations are necessary to implement the spin correlation algorithm in the dipole shower.
A dipole consists of two colour-connected external partons, an emitter and a spectator. In general the spectator defines the soft radiation pattern and, for several types of dipole, is used to absorb the recoil momentum in a splitting. In other cases a set of outgoing particles is used to absorb the recoil momentum in splittings. The algorithm described in this paper does not include any formal treatment for spectators or the recoil in splittings. In each splitting the momentum of one or several particles are therefore modified in some way that is not described by the helicity amplitudes. One consideration that we do make is to ensure that any change in the basis states of a spectator, or other particle, arising from the change in its momentum is accounted for. In Sect. 3 we investigate the effect of the treatment of recoils on predictions made using the dipole shower.
The objects that evolve in the dipole shower are series of colour-connected dipoles, referred to as 'dipole chains' [27]. At a given point in the shower evolution each dipole in the selected chain has a probability to radiate. This is in contrast to the treatment in the angular-ordered shower in which the shower evolution proceeds independently for each parton from a hard process or decay. Extra care must therefore be taken in the dipole shower to ensure that the spin density and decay matrices required for the computation of a given splitting are up-to-date. Furthermore, in the dipole shower, each splitting takes place in a unique frame defined by the momenta of the emitter and spectator. We must therefore compute and apply the mappings described in Appendix B for every splitting in the dipole shower.
The splitting probabilities in the dipole shower are given by the Catani-Seymour splitting kernels [28,29]. These splitting kernels, and the kinematics used to describe splittings in the dipole shower, are written in terms of splitting variables specific to this formulation. On the other hand the helicity amplitudes in Sect. 2.2 and Appendix A are derived using a decomposition of the splitting momenta into the Sudakov basis. In particular the helicity amplitudes are written in terms of the light-cone momentum fraction z. In practice in the dipole shower in Herwig7 we actually generate the variables z and p ⊥ , the transverse momentum of the emission, and map these to the variables required to express the splitting kernels. No additional work is therefore required to use the helicity amplitudes in Tables 2 and 3 in the dipole shower.

Correlations inside the parton shower
The analytic result for the distribution of the angular difference between the planes of successive parton shower branchings is given in Eq. (20). This expression can be expanded for each of the four sequences of branchings that give rise to correlations using the expressions in Table 1. In this section we present the predictions of these distributions obtained using the angular-ordered and dipole parton showers in Her-wig7. The angular difference between two successive parton shower branchings is measured in the splitting frame of the second branching. This test verifies the implementation of the helicity amplitudes in both parton showers. Furthermore, in the dipole shower, it also probes the implementation of the basis state mappings between splittings.
We test the cases of final-state radiation (FSR) and initialstate radiation (ISR) separately. In the ISR case the first splitting is identified as that closest to the incoming hadron and the intermediate parton between the splittings is space-like. In the dipole shower we can divide FSR and ISR further according to the type of dipole considered. Specifically we divide FSR, radiation from a final-state emitter, according to whether the spectator is final-state (final-final) or initial-state (final-initial) and we divide ISR, radiation from an initialstate emitter, according to whether the spectator is final-state (initial-final) or initial-state (initial-initial). We include a separate result for each of these four types of dipole. Note that we do not consider radiation from decay processes in this test.
For the purpose of these tests we implement an artificial restriction on the splittings allowed in the dipole shower. Following the first splitting we only allow subsequent splittings from dipoles in which the spectator is the spectator of the previous splitting and in which the emitter was produced as a new parton in the previous splitting. This restriction has two purposes. First, by forbidding subsequent emissions from different legs of the hard process, it allows us to probe only the correlations in the shower. Second, by using the same spectator in subsequent splittings the frame of the second splitting is a suitable frame in which to measure the angular difference between the splittings. (c) (d) Fig. 1 The analytic result for the difference in azimuthal angle between the branching planes of subsequent final-state a q → qg and g → gg, b q → qg and g → qq, c g → gg and g → gg, and d g → gg and g → qq splittings compared to the distributions predicted using the angular-ordered (QS) and dipole parton showers in Herwig7. The pre-dictions obtained using only final-final (DS-FF) and final-initial (DS-FI) dipoles in the dipole shower are shown separately. The prediction of the angular-ordered (CorrOff) shower without spin correlations is included for comparison. The momentum fractions in the first (z 1 ) and second branchings (z 2 ) lie in the ranges given The results are shown in Figs. 1 and 2, for FSR and ISR respectively. We have chosen to measure the azimuthaldifference for splittings in which the momentum fraction in the first and second branchings lies in the range 0.9 < z 1 < 1.0 and 0.4 < z 2 < 0.5 respectively as this is the configuration in which the correlation is strongest. All of the results shown are for the case of massless quarks.
Each plot shows the analytic result and the parton shower predictions. In each plot we have included the prediction obtained using the angular-ordered shower with spin correlations switched off. In each case this gives rise to a flat line at 1/2π and we have confirmed that the dipole shower also pre-dicts a flat line. All of the parton shower predictions display good agreement with the analytic result in each case.

Correlations with the hard process
In this section we consider results that probe the correlations between the parton shower and the hard process. These tests verify that correlations are passed correctly between the hard process and the parton shower. In addition these tests also probe the treatment of spectators and splitting recoils in the dipole shower, as discussed in Sect. 2.5.

123
(a) (b) (c) (d) Fig. 2 The analytic result for the difference in azimuthal angle between the branching planes of subsequent initial-state a q → qg and g → gg, b q → qg and g → qq, c g → gg and g → gg and d g → gg and g → qq splittings compared to the distributions predicted using the angular-ordered (QS) and dipole parton showers in Herwig7. The pre-dictions obtained using only initial-initial (DS-II) and initial-final (DS-IF) dipoles in the dipole shower are shown separately. The prediction of the angular-ordered (CorrOff) shower without spin correlations is included for comparison. The momentum fractions in the first (z 1 ) and second branchings (z 2 ) lie in the ranges given The analytic result for the distribution of the azimuthal angle between the planes of the g → qq branchings in h 0 → gg → qqq q is given in Eq. (28). This analytic result and the predictions of the angular-ordered and dipole parton showers are shown in Fig. 3. In addition we include the result obtained from a sample of leading-order (LO) events generated using MadGraph5_aMC@NLO [30]. All quarks are massless and our analysis requires two gluon splittings to different quark flavours to enable perfect identification of the quark pairs. The shower predictions both display good agreement with the analytic result and the LO prediction. The differences that remain are due to the cutoff on the transverse momentum used in both parton showers, whereas the analytic result has no cutoff and the LO result includes a cut on the invariant mass of the quark-antiquark pairs which does not affect the shape of the distribution. The transverse momentum cutoff removes more of the region z → 0, 1 where the correlation is smallest giving a slightly larger correlation effect overall.
The above result probes the the treatment of FSR. In order to test the correlations between ISR and the hard process we also consider the Higgs boson production process gg → h 0 followed by the backward splitting of each of the two Fig. 3 The analytic result for the difference in azimuthal angle between the planes of the two branchings in h 0 → gg → qqq q compared to the distributions predicted using the angular-ordered (QS) and dipole (DS) parton showers in Herwig7. The angular-ordered shower (QS-CorrOff) and dipole shower (DS-CorrOff) predictions without spin correlations are included for comparison. The result obtained from a sample of LO events generated using MadGraph5_aMC@NLO (LO) is also shown gluons into an incoming quark and an outgoing quark. In order to obtain a finite leading-order result we require that the minimum transverse momentum of the outgoing quarks is 20 GeV.
The predictions for the distribution of the difference in the azimuthal angle between the planes of the branchings predicted using the angular-ordered and dipole parton showers are shown in Figs. 4 and 5 respectively. In each plot we also include the result obtained from a sample of LO events generated using MadGraph5_aMC@NLO for comparison.
In Fig. 4 we include predictions obtained using the angular-ordered parton shower with and without spin correlations. When spin correlations are not included in the parton shower the predicted distribution is simply a flat line. We find that, with spin correlations included in the parton shower, the angular-ordered parton shower prediction is similar to the LO prediction with some differences in shape due to corrections beyond the collinear limit.
The predictions obtained using the dipole parton shower display more complex behaviour and we have included several results in Fig. 5. We first note that the prediction produced using the dipole parton shower without spin correlations is not flat. This is due to the treatment of splitting recoils. In initial-initial dipoles the recoil in splittings is distributed amongst all outgoing particles other than the new emission, and in initial-final dipoles the outgoing spectator gains a transverse component to its momentum. The momentum of the outgoing quark produced in the first split- Fig. 4 The difference in azimuthal angle between the planes of two initial-state g → qq branchings in gg → h 0 predicted using the angular-ordered (QS) parton shower in Herwig7. The angular-ordered parton shower (QS-CorrOff) prediction without spin correlations is also included. The result obtained from a sample of LO events generated using MadGraph5_aMC@NLO (LO) is shown for comparison ting is therefore changed in a non-trivial way in the second splitting and this gives rise to a directional preference of the second splitting relative to the first splitting. This behaviour necessarily affects the prediction when spin-correlations are included and gives rise to the corresponding distribution in Fig. 5.
In order to demonstrate that the effects seen in the dipole parton shower predictions are indeed due to the treatment of recoil momenta in splittings, we have also included results obtained using a modified version of the dipole parton shower. In this modified parton shower we only allow splittings off initial-initial dipoles and we have modified the behaviour of these splittings such that the splitting recoil is entirely absorbed by the outgoing Higgs boson in both of the splittings. With these modifications the direction of the quark produced in the first splitting is not modified in the second splitting and when spin correlations are not included the predicted distribution is a flat line. As such the prediction with spin correlations included displays better agreement with the angular-ordered parton shower and LO predictions. Again there are differences in shape between the dipole parton shower prediction and the LO prediction due to corrections beyond the collinear limit.
Similar problems with the default recoil scheme in dipole parton showers were recently observed in Ref. [31] where it was shown that the same change in the recoil strategy we have used resolved issues with the logarithmic accuracy of the parton shower.
123 Fig. 5 The difference in azimuthal angle between the planes of two initial-state g → qq branchings in gg → h 0 predicted using the dipole parton shower (DS) in Herwig7. The dipole parton shower (DS-CorrOff) prediction without spin correlations is also included. Predictions obtained using the dipole parton shower restricted to allow branchings from II dipoles only and with a modified handling of splitting recoils, as described in the text, are shown with (DS-II) and without (DS-II-CorrOff) spin correlations. The result obtained from a sample of LO events generated using MadGraph5_aMC@NLO (LO) is shown for comparison

Correlations in decay processes
The spin correlations in the hard process can also affect the distribution of the particles produced in the subsequent decay of unstable particles, such as the top quark, and also give correlations between the decay products of different unstable particles. In these cases the physical effects are dominated by the correlation between the hard process and the decay, which are present even at leading order, and we want to ensure that the parton shower correctly preserves the correlations which are already present, i.e. that the spin information is correctly passed through the parton shower phase. Technically in Herwig this is vital as the decays are generated as part of the parton shower. The simplest example is the production and decay of the top quark. Figure 6 shows the azimuthal separation of the charged leptons in dileptonic pp → tt events at a centre-of-collision energy of 8 TeV, measured by CMS. In addition we show the predictions of this distribution obtained using the angular-ordered and dipole parton showers in Her-wig7. The hard process is produced using a LO matrix element. In the angular-ordered shower the top-quark decays are corrected to NLO in QCD while in the dipole shower no such correction is applied to obtain these predictions. There is reasonable agreement between the experimental result and both parton shower algorithms showing that both approaches Fig. 6 The azimuthal separation of the charged leptons in 8 TeV dileptonic pp → tt events, as measured by CMS [32] and predicted using the angular-ordered (QS) and dipole (DS) parton showers in Herwig7. The predictions of the angular-ordered (QS-CorrOff) shower and the dipole shower (DS-CorrOff) without any spin correlations are also shown. This is dominated by the correlations between the production and decay of the top quark correctly transmit the spin correlations from the production process to the decays.

Conclusions
We have implemented the spin correlation algorithm of Refs. [14][15][16][17][18] in the angular-ordered and dipole parton showers in Herwig7. This feature will be available for public use in Herwig7.2. We have compared the predictions obtained using each of the parton showers in Herwig7 to analytic calculations or predictions obtained using a LO ME. Through these comparisons we have confirmed that the spin correlation algorithm is functioning correctly in both showers.
The handling of splitting recoils in the dipole shower is not formally included in the spin correlation algorithm. We have discussed these limitations and presented results that show where these effects are evident. Despite these limitations, we find that the dipole shower, and the angular-ordered shower, produce a fairly accurate prediction of a spin-correlation sensitive observable, namely the azimuthal separation of the leptons, in pp → tt events.
While spin correlation effects are often unobservable in average distributions, as we have seen there are cases where they are important. Their implementation in Herwig7 is therefore an important part of improving the accuracy of the simulation. Table 3 Spin unaveraged splitting functions for q → qg. In addition to the factors given above each term has a phase e iφ (λ0−λ1−λ2) , where λ 2 = ±1 for the outgoing gluon and λ 0,1 = ± 1 2 for quarks andq is the Herwig++ evolution variable λ 0 λ 1 F λ0,λ1,λ2 λ 2 = + λ 2 = − In this case the helicity amplitudes for the branching are These functions are given in Table 3. 6

A.3 Squark branching
For squarks there is only one branching to consider, i.e. q →qg. The kinematics are the same as those for gluon radiation from a quark, i.e. m 0 = m 1 = m and m 2 = 0. The polarization vectors of the gluons are also the same as for radiation from a quark, Eq. 33. The spin unaveraged splitting function is given by If we sum over the helicities of the gluon we obtain the correct massive splitting function.

B Basis rotation
In the case that we have a spin density matrix from the production process and then wish to generate the correlations in the shower for simplicity we will have to deal with the rotation of the basis states used. Consider a process with production matrix element A prod and a matrix element A branch giving the subsequent time-like branching of the intermediate particle with momentum p and mass m. 6 If we make the same redefinition as before these functions agree with those in [15] Table 1 up to the overall sign which is irrelevant.
In the case of a fermion we will replace the spin sum with the Dirac spin sum, i.e p +m = a Ψ aΨa , or similarly for vector bosons −g μν + p μ n ν + p ν n μ p · n = a μ a * ν a We write both of these as where the basis state for the production is χ † a and for the branching is χ a such that the matrix element for the full process is using the Einstein convention. The matrix element squared is then If we normalise by the spin summed matrix element for the production process then we obtain where M prod = χ † a A prod and In order to simplify the calculation we have made a specific choice of the basis of the wavefunctions χ a when we calculate the branching, which may not be the same as those used in the calculation of the production, χ a , after Lorentz transformation into the frame in which the branching is calculated.
We therefore need to calculate the rotation between these two choices of basis states