1 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.Footnote 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 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.Footnote 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.

2 Algorithm

2.1 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\rightarrow gg\rightarrow q\bar{q}q'\bar{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

$$\begin{aligned} \rho _{g_1}^{\lambda _{g_1}\lambda _{g_1}'} = \frac{1}{N}\mathcal {M}^{\lambda _{g_1}\lambda _{g_2}}_{h^0\rightarrow gg}{\mathcal {M}^*}^{\lambda '_{g_1}\lambda _{g_2}}_{h^0\rightarrow gg}, \end{aligned}$$
(1)

where \(\mathcal {M}_{h^0\rightarrow gg}\) is the helicity amplitude for the decay of the Higgs boson into two gluons, \(\lambda _{g_1}\) and \(\lambda '_{g_1}\) are the polarizations of the first gluon and \(\lambda _{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 \(\mathrm{Tr}\,\rho _{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.Footnote 3 However the azimuthal angle of the emission is affected by the spin density matrix and is generated using the distribution

$$\begin{aligned} \rho _{g_1}^{\lambda _{g_1}\lambda _{g_1}'} \mathcal {M}^{\lambda _{g_1}\lambda _q\lambda _{\bar{q}}}_{g\rightarrow q\bar{q}}\mathcal {M}^{*\lambda '_{g_1}\lambda _q\lambda _{\bar{q}}}_{g\rightarrow q\bar{q}}, \end{aligned}$$
(2)

where \(\mathcal {M}_{g\rightarrow q\bar{q}}\) is the helicity amplitude for the splitting of a gluon into a quark–antiquark pair and \(\lambda _q\) and \(\lambda _{\bar{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.

$$\begin{aligned} \rho _q^{\lambda _q\lambda '_q} = \frac{1}{N}\rho _{g_1}^{\lambda _{g_1}\lambda _{g_1}'} \mathcal {M}^{\lambda _{g_1}\lambda _q\lambda _{\bar{q}}}_{g\rightarrow q\bar{q}}\mathcal {M}^{*\lambda '_{g_1}\lambda '_q\lambda _{\bar{q}}}_{g\rightarrow q\bar{q}}, \end{aligned}$$
(3)

where again the normalization is such that \(\mathrm{Tr}\,\rho _q=1\) and is equal to the spin contracted amplitude used to calculate the distribution of the \(g\rightarrow q\bar{q}\) 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

$$\begin{aligned} D^{\lambda _{g_1}\lambda '_{g_1}}_{g_1} = \frac{1}{N} \mathcal {M}^{\lambda _{g_1}\lambda _q\lambda _{\bar{q}}}_{g\rightarrow q\bar{q}}\mathcal {M}^{*\lambda '_{g_1}\lambda _q\lambda _{\bar{q}}}_{g\rightarrow q\bar{q}}, \end{aligned}$$
(4)

where again the normalization is such that \(\mathrm{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.

$$\begin{aligned} \rho _{g_2}^{\lambda _{g_2}\lambda '_{g_2}} = \frac{1}{N}\mathcal {M}^{\lambda _{g_1}\lambda _{g_2}}_{h^0\rightarrow gg}{\mathcal {M}^*}^{\lambda '_{g_1}\lambda '_{g_2}}_{h^0\rightarrow gg} D_{g_1}^{\lambda _{g_1}\lambda '_{g_1}}, \end{aligned}$$
(5)

where again the normalization is such that \(\mathrm{Tr}\rho _{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\rightarrow 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

$$\begin{aligned} \mathcal {M}^{\lambda _{g_1}\lambda _{g_2}}_{h^0\rightarrow gg}{\mathcal {M}^*}^{\lambda '_{g_1}\lambda '_{g_2}}_{h^0\rightarrow gg} \mathcal {M}^{\lambda _{g_1}}_{g\rightarrow q\bar{q}}\mathcal {M}^{*\lambda '_{g_1}}_{g\rightarrow q\bar{q}} \mathcal {M}^{\lambda _{g_2}}_{g\rightarrow q'\bar{q'}}\mathcal {M}^{*\lambda '_{g_2}}_{g\rightarrow q'\bar{q'}}, \end{aligned}$$
(6)

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.

2.2 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 quasi-collinear 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

$$\begin{aligned} p = (\sqrt{\mathbf{p}^2+m_0^2}; 0,0,\mathbf{p}), \end{aligned}$$
(7)

where \(m_0\) is the on-shell mass of the branching particle and \(\mathbf{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

$$\begin{aligned} q_i = \alpha _ip+\beta _in +q_{\perp i}, \end{aligned}$$
(8)

where p is the momentum of the branching particle before any emission, n is a light-like reference vector, \(\alpha _i\) and \(\beta _i\) are coefficients, and \(q_{\perp 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\cdot q_{\perp i}=n\cdot q_{\perp i}=0\) and \(p\cdot n>0\). As we need an explicit representation we take

$$\begin{aligned} n = (1,0,0,-1). \end{aligned}$$
(9)

If we consider a branching \(0\rightarrow 1,2\) then in this parameterization the momenta of the particles produced in the branching are:

$$\begin{aligned} q_1= & {} zp+\beta _1n +q_\perp ; \end{aligned}$$
(10a)
$$\begin{aligned} q_2= & {} (1-z)p+\beta _2n -q_\perp ; \end{aligned}$$
(10b)

where z is the light-cone momentum fraction of the first parton and

$$\begin{aligned} q_\perp= & {} (0;p_\perp \cos \phi ,p_\perp \sin \phi ,0); \end{aligned}$$
(10c)
$$\begin{aligned} \beta _1= & {} \frac{1}{2zp\cdot n}\left( p_\perp ^2+m_1^2-z^2m_0^2\right) ; \end{aligned}$$
(10d)
$$\begin{aligned} \beta _2= & {} \frac{1}{2(1-z)p\cdot n}\left( p_\perp ^2+m_2^2-(1-z)^2m_0^2\right) ; \end{aligned}$$
(10e)

where \(m_{1,2}\) are the on-shell masses of the particles produced, \(p_\perp \) is the magnitude of the transverse momentum of the branching and \(\phi \) is the azimuthal angle for the first particle produced in the branching.

In the Herwig7 angular-ordered parton shower we use the evolution variable

$$\begin{aligned} \tilde{q}^2 = \frac{q^2_0-m^2_0}{z(1-z)}, \end{aligned}$$
(10f)

for final-state evolution.

The momentum of the off-shell particle initiating the branching is

$$\begin{aligned} q_0 = p +\beta _0n, \end{aligned}$$
(10g)

where

$$\begin{aligned} \beta _0 = \beta _1+\beta _2 = \frac{1}{2p\cdot n}\left( \frac{p_\perp ^2}{z(1-z)}+\frac{m_1^2}{z}+\frac{m_2^2}{1-z}-m_0^2\right) , \end{aligned}$$
(10h)

such that the virtuality of the branching parton is

$$\begin{aligned} q_0^2 = \frac{p_\perp ^2}{z(1-z)}+\frac{m_1^2}{z}+\frac{m_2^2}{1-z}. \end{aligned}$$
(10i)

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 \(\lambda \) and expanding in \(\lambda \) [24], i.e.

$$ \begin{aligned} m_i\rightarrow&\lambda m_i \& p_\perp \rightarrow \lambda p_\perp . \end{aligned}$$
(11)

In order to define our conventions for spinors and polarization vectors we present the calculation of the \(g\rightarrow q\bar{q}\) 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:

$$\begin{aligned} p= & {} \mathbf{p}(1; 0, 0, 1); \end{aligned}$$
(12)
$$\begin{aligned} \epsilon ^\mu _{\lambda _0}= & {} \frac{1}{\sqrt{2}}(0; -\lambda _0, -i, 0); \end{aligned}$$
(13)

where \(\lambda _0=\pm 1\) is the helicity of the gluon. For the branching \(g\rightarrow q \bar{q}\) , \(m_1=m_2=m\), therefore expanding in \(\lambda \) the momenta of the particles in the branching are:

$$\begin{aligned} q_0= & {} \left( \mathbf{p}+\frac{m^2+{p_\perp }^2}{4\mathbf{p}z\left( 1-z\right) }\lambda ^2 , 0,0, \mathbf{p}-\frac{m^2+p_\perp ^2}{4\mathbf{p}z\left( 1-z\right) }\lambda ^2\right) \nonumber \\&+\mathcal {O}(\lambda ^3); \end{aligned}$$
(14a)
$$\begin{aligned} q_1= & {} \left( z\mathbf{p}+ \frac{m^2+p_\perp ^2}{4z\mathbf{p}}\lambda ^2, \lambda p_\perp \cos \phi , \lambda p_\perp \sin \phi ,\right. \nonumber \\&\left. z\mathbf{p}-\frac{m^2+p_\perp ^2}{4z\mathbf{p}}\lambda ^2\right) +\mathcal {O}(\lambda ^3);\end{aligned}$$
(14b)
$$\begin{aligned} q_2= & {} \left( \left( 1-z\right) \mathbf{p}+\frac{m^2+p_\perp ^2}{4\left( 1-z\right) \mathbf{p}}\lambda ^2,-\lambda p_\perp \cos \phi , -\lambda p_\perp \sin \phi ,\right. \nonumber \\&\qquad \times \left. \left( 1-z\right) \mathbf{p}-\frac{m^2+p_\perp ^2}{4\left( 1-z\right) \mathbf{p}}\lambda ^2 \right) +\mathcal {O}(\lambda ^3). \end{aligned}$$
(14c)
Table 1 The coefficients A and B for the correlation between the azimuthal angles of the first and second branching. The momentum fraction in the first branching is \(z_1\) and \(z_2\) in the second branching
Full size table

In general we define the spinors using the conventions in [26]. This allows us to write the spinors for the outgoing particles:

$$\begin{aligned} \bar{u}_{\frac{1}{2}}(q_1)= & {} \left[ \sqrt{2\mathbf{p}z}, \frac{\mathrm{e}^{-i\phi }p_\perp }{\sqrt{2\mathbf{p}z}}\lambda , \frac{m}{\sqrt{2\mathbf{p}z}}\lambda , \frac{ p_\perp m\lambda ^2 \mathrm{e}^{-i\phi }}{\left[ 2z\mathbf{p}\right] ^{3/2}} \right] ;\nonumber \\ \bar{u}_{-\frac{1}{2}}(q_1)= & {} \left[ -\frac{p_\perp m \lambda ^2\mathrm{e}^{i\phi }}{\left[ 2z\mathbf{p}\right] ^{3/2}}, \frac{m}{\sqrt{2\mathbf{p}z}}\lambda , -\frac{\mathrm{e}^{i\phi }p_\perp }{\sqrt{2\mathbf{p}z}}\lambda , \sqrt{2\mathbf{p}z} \right] ;\nonumber \\ \end{aligned}$$
(15a)

and

$$\begin{aligned}&v_{\frac{1}{2}}(q_2) = \left[ \begin{array}{c}-\frac{p_\perp \mathrm{e}^{-i\phi } }{\sqrt{2\mathbf{p}(1-z)}}\lambda \\ -\sqrt{2\mathbf{p}(1-z)}\\ {\frac{mp_\perp \lambda ^2 \mathrm{e}^{-i\phi } }{\left[ (2\mathbf{p}(1-z)\right] ^{3/2}}}\\ \frac{m}{\sqrt{2\mathbf{p}(1-z)}}\lambda \end{array}\right] ; \nonumber \\&\quad v_{-\frac{1}{2}}(q_2) = \left[ \begin{array}{c}\frac{m}{\sqrt{2\mathbf{p}(1-z)}}\lambda \\ -\frac{mp_\perp \lambda ^2 \mathrm{e}^{i\phi } }{\left[ 2\mathbf{p}(1-z)\right] ^{3/2}}\\ -\sqrt{2\mathbf{p}(1-z)}\\ \frac{p_\perp \mathrm{e}^{i\phi } }{\sqrt{2\mathbf{p}(1-z)}}\lambda \end{array}\right] . \end{aligned}$$
(15b)

We define the helicity amplitudes \(F_{\lambda _0\lambda _1\lambda _2}\) for the branching such that the spin averaged splitting function is

$$\begin{aligned} P(z) = \frac{1}{2}\sum _{\lambda _0,\lambda _1,\lambda _2} F_{\lambda _0\lambda _1\lambda _2} F^*_{\lambda _0\lambda _1\lambda _2}. \end{aligned}$$
(16)

For the \(g\rightarrow q\bar{q}\) branching

(17)

There is no simple form for all the helicities and therefore the functions, \(F_{\lambda _0\lambda _1\lambda _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.Footnote 4

2.3 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\rightarrow 12\) followed by \(2\rightarrow 34\). The only non-zero correlation is due to the polarization of an intermediate gluon.

If we consider the branching \(q\rightarrow qg\) the spin density matrix for the radiated gluon is

$$\begin{aligned} \rho _g = \left( \begin{array}{cc}\ \frac{1}{2} &{} -\frac{z_1\mathrm{e}^{2i\phi _1}}{1+z_1^2}\\ -\frac{z_1\mathrm{e}^{-2i\phi _1}}{1+z_1^2} &{} \frac{1}{2}\end{array}\right) , \end{aligned}$$
(18)

where \(z_1\) is the momentum fraction and \(\phi _1\) the azimuthal angle of the quark, respectively. We have neglected the mass of the radiating quark.

Similarly for the branching \(g\rightarrow gg\) the spin density matrix for the radiated gluon is

$$\begin{aligned} \rho _g = \left( \begin{array}{cc}\ \frac{1}{2} &{}-\frac{z_1^2\mathrm{e}^{ 2i\phi _1}}{2\left( 1-z_1(1-z_1)\right) ^2}\\ -\frac{z_1^2\mathrm{e}^{-2i\phi _1}}{2\left( 1-z_1(1-z_1)\right) ^2} &{} \frac{1}{2}\end{array}\right) . \end{aligned}$$
(19)

If we contract these spin density matrices with the appropriate helicity amplitudes for the subsequent branching of the gluon then we obtain the distribution

$$\begin{aligned} \frac{1}{2\pi }\left[ 1+A\,B\cos 2\varDelta \phi \right] , \end{aligned}$$
(20)

where \(\varDelta \phi =\phi _2-\phi _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\rightarrow g g \rightarrow q \bar{q} q'\bar{q'}\), which we have already considered.

We take the amplitude to beFootnote 5

$$\begin{aligned} \mathcal {M}_{h^0\rightarrow gg} = -p_1\cdot p_2 \epsilon ^*_1\cdot \epsilon _2^* +p_2\cdot \epsilon _1^* p_1\cdot \epsilon _2^*, \end{aligned}$$
(21)

where \(p_{i=1,2}\) and \(\epsilon _{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\rightarrow g g\) are:

$$\begin{aligned} \mathcal {M}_{h^0\rightarrow gg}^{++}= & {} -\frac{m_h^2}{2}\mathrm{e}^{-2i\phi }; \end{aligned}$$
(22a)
$$\begin{aligned} \mathcal {M}_{h^0\rightarrow gg}^{--}= & {} -\frac{m_h^2}{2}\mathrm{e}^{ 2i\phi }; \end{aligned}$$
(22b)

where \(\phi \) 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

$$\begin{aligned} D = \left( \begin{array}{cc} \frac{1}{2} &{} a_1(z_1,\tilde{q}_1){\mathrm{e}^{2\,i\phi _{1}}}\\ a_1(z_1,\tilde{q}_1){\mathrm{e}^{-2\,i\phi _{1}}} &{} \frac{1}{2} \end{array}\right) , \end{aligned}$$
(23)

where \(\phi _1\) is the azimuthal angle of the quark q in a frame in which the 1st gluon is along the z-axis and

$$\begin{aligned} a_1(z_1,\tilde{q}_1)= \frac{ \left( z_1 \left( 1-z_1 \right) -{\frac{m_q^2}{z_1 \left( 1-z_1 \right) {\tilde{q}_1}^2}} \right) }{\left( 1- 2z_1 \left( 1-z_1 \right) +\frac{2m_q^2}{z_1 \left( 1-z_1 \right) {\tilde{q}_1}^2} \right) }, \end{aligned}$$
(24)

where \(z_1\) is the energy fraction of the quark, \(\tilde{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,

$$\begin{aligned} \rho =\left( \begin{array}{cc} \frac{1}{2} &{} a_1(z_1,\tilde{q}_1){\mathrm{e}^{2i\phi _{1}-4i\phi }}\\ a_1(z_1,\tilde{q}_1){\mathrm{e}^{-2i\phi _{1}+4i\phi }} &{} \frac{1}{2} \end{array}\right) . \end{aligned}$$
(25)

This gives the distribution

$$\begin{aligned} \left( 1-2z_2 \left( 1-z_2 \right) +\frac{2m^2_{q'}}{z_2 \left( 1-z_2 \right) {\tilde{q}_2}^2} \right) \nonumber \\ \times \left( 1+4a_1(z_1,\tilde{q}_2)a_{2}(z_2,\tilde{q}_2)\cos \left( 4\phi -2\phi _2-2\phi _1\right) \right) , \end{aligned}$$
(26)

where \(z_2\) is the energy fraction of the quark \(q'\) and \(\tilde{q}_2\) the evolution variable for the second branching. The azimuthal angle of the second branching, \(\phi _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 \(\phi '_1=2\phi -\phi _1+\pi \).

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

$$\begin{aligned} P_{g\rightarrow q\bar{q}}(z_1)P_{g\rightarrow q\bar{q}}(z_2) \left[ 1+4a_1(z_1)a_{2}(z_2)\cos \left( 2\left[ \phi _2-\phi '_1\right] \right) \right] . \end{aligned}$$
(27)

Integrating over \(z_{1,2}\) between 0 and 1 and normalizing the resulting distribution gives

$$\begin{aligned} \frac{1}{8\pi } \left[ 3+2\cos ^2\left( \varDelta \phi \right) \right] , \end{aligned}$$
(28)

where \(\varDelta \phi =\phi _2-\phi '_1\) is the azimuthal angle between the planes of the two branchings.

2.4 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.

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.

2.5 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_\perp \), 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.

Fig. 1
figure 1

The analytic result for the difference in azimuthal angle between the branching planes of subsequent final-state a \(q\rightarrow qg\) and \(g\rightarrow gg\), b \(q\rightarrow qg\) and \(g\rightarrow q\bar{q}\), c \(g\rightarrow gg\) and \(g\rightarrow gg\), and d \(g\rightarrow gg\) and \(g\rightarrow q\bar{q}\) splittings compared to the distributions predicted using the angular-ordered (QS) and dipole parton showers in Herwig7. The predictions 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

Full size image

3 Results

3.1 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 Herwig7. 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 initial-state 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 initial-state 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.

The results are shown in Figs. 1 and 2, for FSR and ISR respectively. We have chosen to measure the azimuthal-difference 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\pi \) and we have confirmed that the dipole shower also predicts a flat line. All of the parton shower predictions display good agreement with the analytic result in each case.

Fig. 2
figure 2

The analytic result for the difference in azimuthal angle between the branching planes of subsequent initial-state a \(q\rightarrow qg\) and \(g\rightarrow gg\), b \(q\rightarrow qg\) and \(g\rightarrow q\bar{q}\), c \(g\rightarrow gg\) and \(g\rightarrow gg\) and d \(g\rightarrow gg\) and \(g\rightarrow q\bar{q}\) splittings compared to the distributions predicted using the angular-ordered (QS) and dipole parton showers in Herwig7. The predictions 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

Full size image

3.2 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.

The analytic result for the distribution of the azimuthal angle between the planes of the \(g\rightarrow q\bar{q}\) branchings in \(h^0\rightarrow g g \rightarrow q \bar{q} q'\bar{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\rightarrow 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 \rightarrow h^0\) followed by the backward splitting of each of the two 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.

Fig. 3
figure 3

The analytic result for the difference in azimuthal angle between the planes of the two branchings in \(h^0\rightarrow g g \rightarrow q \bar{q} q'\bar{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

Full size image

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 partonshower 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 splitting 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 indipole 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.

Fig. 4
figure 4

The difference in azimuthal angle between the planes of two initial-state \(g \rightarrow q\bar{q}\) branchings in \(g g \rightarrow 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

Full size image
Fig. 5
figure 5

The difference in azimuthal angle between the planes of two initial-state \(g \rightarrow q\bar{q}\) branchings in \(g g \rightarrow 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

Full size image

3.3 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\rightarrow t\bar{t}\) 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 Herwig7. 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 correctly transmit the spin correlations from the production process to the decays.

Fig. 6
figure 6

The azimuthal separation of the charged leptons in 8 TeV dileptonic \(pp\rightarrow t\bar{t}\) 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

Full size image

4 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\rightarrow t\bar{t}\) 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.